## 1. textbook 12.8
library('lavaan')
lower <- "
0,
9, 0,
3, 7, 0,
6, 5, 9, 0,
11, 10, 2, 8, 0,"
D = as.dist(getCov(lower, names=c("1","2","3","4","5")))
hc.avg <- hclust(D, method="average")
hc.s <- hclust(D, method="single")
par(mfrow=c(1,2))
plot(hc.avg, hang=-1, main="Cluster Dendrogram (Average linkage)")
plot(hc.s, hang=-1, main="Cluster Dendrogram (Single linkage)")
## 2. textbook 12.9
library('lavaan')
lower <- "
0,
.76, 0,
2.97, .8, 0,
4.88, 4.17, .21, 0,
3.86, 1.92, 1.51, .51, 0,"
D = as.dist(getCov(lower, names=c("1","2","3","4","5")))
par(mfrow=c(2,2))
plot(hc.s <-hclust(D, method="single"),main="Single linkage") # single linkage
plot(hc.c <-hclust(D, method="complete"),main="Complete linkage") # complete linkage
plot(hc.avg <-hclust(D, method="average"),main="Average linkage") # average linkage
## 3.
x <- data.frame(X1=c(5,-1,1,-3),X2=c(3,1,-2,-2),row.names=c(LETTERS[1:4]))
k.cl <- kmeans(x, 2)
k.cl$cluster
## 4.
FOODP <- read.table("FOODP.DAT", skip=2, header = TRUE)
prin.cr <- prcomp( ~Bread+Hamburger+Butter+Apples+Tomatoes, data=FOODP, scale=TRUE)
pc.score <- prin.cr$x[,1:2] # get pc scores of the first 2 dim.
pc.score <- -pc.score # 調整正負號以便和價格高低相符
hc <- hclust(dist(pc.score),method="average")
groups <- cutree(hc,k=1:10)
# (a)
#---------------------------------------------------------
R.square <- function(x, method="average", k){
### R-square for hierarchical clustering given number of cutting groups
### x : original data (col=variables, row=observations)
### method: hierarchical clustering linkage method
### k : number of cutting groups
hc <- hclust(dist(x), method=method)
group <- as.matrix(cutree(hc,k)) # get the group vectors
## sum of squares for all the elements
ss <- function(x) sum(scale(x, scale = FALSE)^2)
## calculate BetweenSS, totalWithinSS, and R-squares
betweenss <- NULL # reserve space
tot.withinss <- NULL
r.square <- NULL
for(i in 1:length(k)){
## calculate cluster center "fitted" to each obs.:
fitted.x <- NULL # reserve space
for(n in 1:nrow(x)){
g <- group[n,i] # get the group for each obs.
fitted.x <- rbind(fitted.x, colMeans(x[group[,i]==g, , drop = FALSE]))
}
resid.x <- x - fitted.x
betweenss <- rbind(betweenss, ss(fitted.x))
tot.withinss <- rbind(tot.withinss, ss(resid.x))
r.square <- rbind(r.square, ss(fitted.x)/ss(x))
}
return(data.frame(K.groups = k,
betweenss = betweenss,
tot.withinss = tot.withinss,
totss = rep(ss(x),length(k)),
r.square = r.square))
}
#-------------------------------------------------------
RR <- data.frame(
Single = R.square(pc.score, "single", k=1:10)$r.square,
Complete = R.square(pc.score, "complete", k=1:10)$r.square,
Average = R.square(pc.score, "average", k=1:10)$r.square,
Ward = R.square(pc.score, "ward", k=1:10)$r.square)
RR
R2 <- read.table("R-squares.txt", header = TRUE)
plot(1:10,RR[,1], pch=15, type="b", xaxp=c(0,nrow(R2),nrow(R2)),
xlab="Number of clusters", ylab="R-squares",panel.first = grid())
points(1:10,RR[,2], pch=16, col=2, type="b")
points(1:10,RR[,3], pch=17, col=3, type="b")
points(1:10,RR[,4], pch=18, col=4, type="b")
legend(x="bottomright",legend=names(RR),
lty=rep(1,4), col=c(1:4), pch=c(15:18),bg="white")
# (b)
cl.center <- function(x, method, k){
hc <- hclust(dist(x), method=method)
memb <- cutree(hc, k) # get the group vectors
cent <- NULL # reserve space
for(j in 1:k){
cent <- rbind(cent, colMeans(x[memb == j, , drop = FALSE]))
}
return(cent)
}
centroid <- cl.center(pc.score,"single",k=4)
k.cl <- kmeans(pc.score, centers=centroid)
k.cl$cluster
table(k.cl$cluster) # 各群的人數
## 微調 k.cl$cluster[16] <- 4
# 先利用plot( , type="n" )將資料的範圍先畫出來，
# 再利用text( )把每一個資料的名稱都點出
plot(pc.score, type="n")
text(pc.score, col=k.cl$cluster, labels=row.names(pc.score))
abline(h=0, v=0, lty="dashed")
![多變量分析 [R]：學習筆記](https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhmm_jAdWd28bPWK4_1pigsXQglB2L4IhMz5xzAzM0pdsJ_KA5r3w6c2rPxAhpjCzxym9SJiSKNSCi1p7_pYGh7lj8g9Yop7jhOK92hVBS0ZLNPVcgqe1gd7wB3uiJhr_UoR_azj42J0GY/s1600/R+blogger+3.png)
