############################################################
#Borgen, F. H., & Barnett, D. C. (1987). 
#Applying cluster analysis in counseling psychology research.
#Journal of Counseling Psychology, 34, 456-468.
############################################################
#----------------------#
#Step1: Table1 (pp.456-7)
#----------------------#

##(a) データの読み込み
data1 <- read.csv("http://blue.zero.jp/yokumura/R/cluster/Borgen1987.csv",header=T, row.names="Student")
data1			#Table1



#----------------------#
#Step2: Figure 1 (pp.456-7)
#----------------------#
plot(Academic.Comfort~Introversion.Extraversion,data=data1, type="n", 
	xlim=c(30,80), ylim=c(10,80))
text(rownames(data1), y=data1[,1], x=data1[,2])
graphics.off()


#----------------------#
#Step3: squared Euclidean distance between students 1 and 2 (p.456)
#----------------------#
data1[1,]
data1[2,]
data1[1,] - data1[2,]
(data1[1,] - data1[2,])^2
sum((data1[1,] - data1[2,])^2)



#----------------------#
#Step4: 距離 Table2 (pp.456-8)
#----------------------#

##(a) ユークリッド距離の2乗をループで計算
proximity <- matrix(NA,ncol=nrow(data1), nrow=nrow(data1))

for(i in 1:nrow(data1)){
	for(j in 1:nrow(data1)){
		if(i != j){
			proximity[i,j] <- sum((data1[i,] - data1[j,])^2)
		}
	}
}


##(b) 関数使用
euc1 <- dist(data1, method="euclidean")
euc1
sqrt(proximity)			#平方根を取ることで一致する



#----------------------#
#Step5: Ward's method (pp.456-8)
#----------------------#

##(a) 小さな実例で確認
ward <- proximity			#wardという行列を作る
ward[upper.tri(ward)] <- NA		#上三角をNAとする

##ここから
(minimum <- min(ward, na.rm=T))				#最小値を求める
(minimum.number <- which(ward==minimum, arr.ind=T))	#最小値となる行列番号を得る
ward[which(ward==minimum)] <- NA			#最小値を取った行列をNAとする
##ここまで，何度も，繰り返すとWard's Methodが実感できる


##(b) 関数使用
clust1 <- hclust(euc1, method="ward")
clust1
plot(clust1,xlab="Student",sub=NA, main=NA)

graphics.off()


#----------------------#
#Step6: Table 3 Figure 4 (p.459)
#----------------------#
clus.group <- cutree(clust1, k = 4)
table(clus.group)


attach(data1)
   #subgroup
   (acad.mean <- tapply(Academic.Comfort,clus.group, mean))
   tapply(Academic.Comfort,clus.group, sd)
   (Intro.mean <- tapply(Introversion.Extraversion,clus.group, mean))
   tapply(Introversion.Extraversion,clus.group, sd)

   ##Total
   mean(Academic.Comfort)
   sd(Academic.Comfort)
   mean(Introversion.Extraversion)
   sd(Introversion.Extraversion)
   
   matplot(rbind(acad.mean,Intro.mean),type="l", lty=1,lwd=2,
   ylim=c(10,70), ylab="Standard Score", xlim=c(.5,2.5),xaxt="n",
   )
   axis(1,at=c(1,2), label=c("Academic", "Intro-Extra"))
   legend(2.2,60,paste("Group",1:4), col=1:4, lty=1)
detach(data1)

graphics.off()