########################################################################################
########################################################################################
#I. Cluster analysis:  
#    Perform a product segmentation (group products by similar characteristics),
#    youll need the "cluster" package which has routines of source code that will perform
#    the analysis 
########################################################################################
########################################################################################

#A. Packages can be installed from the menu at the top of your screens R session:
#	Packages / set CRAN mirror (choose first entry, 0-Cloud)
#   Packages / install packages (choose "cluster" from the alphabetical list)
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

#	1. OTHER PACKAGE(S) WE will use... for the second package, use the 'installLoadPkg' function above 
#		a. read in source code
source("http://www.ceresanalytics.com/R_for_IBF_BootCamp/installLoadPkg.r")
source("http://www.ceresanalytics.com/R_for_IBF_BootCamp/clusEval.R")

#		b. run the function 
installLoadPkg("cluster")
installLoadPkg("fpc")


#B. read the data 
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

#A. From the web, load this file: ________________
#	1. for ease of access, assign the filename to a character variable 
#	   this can be a web reference (as shown) or a filename using unix-like syntax, 
#      e.g. "C:/users/sara/myFile.csv"
clusterFilename<-"http://ceresanalytics.com/R_for_IBF_BootCamp/clusterData.txt"

#	2. then read it in 
#      the read.delim functions first argument is file... 
#      you can either specify parameters by name, or place them in the order read

tryCatch({
	clusterData<-read.table(clusterFilename,sep="\t",header=TRUE)
}, error=function(e) {
	print(paste("YOU HAVE NO FILE NAMED",clusterFilename))
})

#	1. see whats in your workspace.  This command is like the unix/linux "ls" command
ls() 

#C. explore the data 
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

#	1. frequency counts--as above, where noted like SAS PROC FREQ 
#		a. for exploratory data analysis, tabulate the dataset by country
#          the "table" function
table(clusterData$OldGroup)

#		b. now a cross-tabulation of frequency counts
table(clusterData$OldGroup,clusterData$OldGroupNum,deparse.level=2)


#	3. frequency histograms of numeric columns 
#		a. how many old groups are there?
myNumOldGroups<-length(table(clusterData$OldGroup))
print(myNumOldGroups)

#		b. make an array of colors ... this requires package "grDevices", which *should* be part of default installation
myNumericColors<-rainbow(myNumOldGroups)
print(myNumericColors)

#		c. just for kicks, show your colors in a pie chart, with all slices of equal width
#			ii. run it, with 'col' parameter and 'main' parameter for titles
pie(x=rep(1,3),col=myNumericColors,main="A pie chart with my rainbow of colors!")


#D. now the cluster analysis... use a hierarchical method
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
#	1. check out what arguments youll need for your method
?agnes

#	2. apply the method 
#		a. vars to use
varsToUse<-c("Mon_AvgSales","Freq_NumPurchPastYear","Recency_TimeSinceLastPchs")
clusResult<-agnes(
	x=clusterData[,varsToUse],
	diss=FALSE,
	keep.diss=TRUE, #the dissimilarity matrix contains the distances between observations
	metric="euclidean",  #sum-of-squares
	stand=TRUE,    		 #make all variables carry equal weight
	method="ward"  		 #based in analysis of variance (ANOVA), good with binary variables
)
#	3. see the tree ... a conversion to "hclust" type simplifies the plot 
#		the height of the legs reflects the strength of the division 
plot(as.hclust(clusResult))

#	4. now for labels ... this might resolve differently for different versions of R
plot(as.hclust(clusResult),labels=clusterData$CustomerID)

#	5. "cut" the tree to assign clusters at your chosen level of aggregation
#		note that the order here is that of the original data frame
myClusAssign<-cutree(tree=clusResult,k=3)  #"k" parameter is how many clusters.  You can see alternatives with ?cutree
print(myClusAssign)
#		a. to see this vertically
	
#	5. 	double-check statistics
cluster.stats(
	d = clusResult$diss, 
	clustering=myClusAssign
)

#E. function to test a different number of clusters, and create 
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
myClusEval<-clusEval(clusResult,2,8)
print(myClusEval)
	
#F. Graphics to evaluate clusters 	
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
#	0. solution names
all_solutionNames<-names(myClusEval)

#	1. knee of curve--within:between ratio 
all_WBRatio<-sapply(myClusEval,function(anEval) anEval$WBRatio)
plot(all_WBRatio,
	xaxt="n",
	pch=19,
	type="l",
	xlab="",
	main="Within:Between Ratio--Choose Knee of Curve")
axis(1,at=1:length(all_WBRatio),labels=all_solutionNames,las=2)

#	2. avg silhouette
all_avgSil<-sapply(myClusEval,function(anEval) anEval$AvgSilWidth)
plot(all_avgSil,
	xaxt="n",
	pch=19,
	type="l",
	xlab="",
	main="Average Silhouette Score--Prefer Maximum")
axis(1,at=1:length(all_WBRatio),labels=all_solutionNames,las=2)
	
#	3. individual silhouette
mySil<- mySil<-silhouette(x=myClusAssign,dist=clusResult$diss,plot=T)
plot(mySil)

#G. Learn about your clusters 
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 

#	1. append the cluster to the record
afterClusterData<-cbind(clusterData,myClusAssign)
#	2. see your data 
head(afterClusterData)

#	3. set up for plots
#		a. rows will be clusters 
#			i. how many clusters  
howManyClusters<-max(myClusAssign)
myClusNums<-c(1:howManyClusters)
#			ii. give each cluster a different color 
myClusColors<-rainbow(howManyClusters)
#		b. columns will be vars 
myClusVars<-colnames(clusResult$data)
#		c. pie plot
clusterSizes<-table(afterClusterData$myClusAssign)
clusterSizePcts<-prop.table(table(afterClusterData$myClusAssign))
pie(clusterSizes,
	col=myClusColors,
	main="Cluster Sizes",
	labels=paste0("Cluster ",names(clusterSizes),":  ",clusterSizes," (",100*round(clusterSizePcts,3),"% )"),
	clockwise=TRUE
)

#	4. plot altogether

#		a. how many bars in histograms 
numBars<-20
#		b. the plot 
par(mfrow=c(length(myClusVars),max(myClusAssign)))
#		c. iterate thru clusters 
for (i in 1:howManyClusters) {
#			i. separate data only for one cluster at a time 
	tempData<-afterClusterData[afterClusterData$myClusAssign==i,]
#			ii. make a histogram of that data 
	j=0
	for (varName in myClusVars) {
#				a. to compare clusters, horizontal axis must be same for all 
		allMin<-min(afterClusterData[,varName])
		allMax<-max(afterClusterData[,varName])
		j=j+1
#				b. now the histogram 
		hist(
			tempData[,varName],
			col=myClusColors[j],
			xlim=c(allMin,allMax),
			main=paste(varName,"Cluster",i),
			cex.main=1.4,
			breaks=seq(allMin,allMax,(allMax-allMin)/numBars)
		)
	}
}