15.3.3 Model-based Methods

# Fig. 15.4.: Mclust results in PCA space library(StatDA) data(chorizon) el=c("Ca","Zn","Mg","Cu","Sr","Na") x=chorizon[,el] sel=((chorizon[,"LITO"]==7) | (chorizon[,"LITO"]==9) | (chorizon[,"LITO"]==82)) grp=chorizon[sel,"LITO"] grp[grp==7] <- 3 grp[grp==9] <- 4 grp[grp==82] <- 20 x=scale(log10(x[sel,])) lab=chorizon[sel,"LITO"] res2=Mclust(x,5:6) res2cl=res2$cl res2cl[res2$cl==5] <- 20 res.pca=princomp(x) expvar=(res.pca$sdev^2)/sum(res.pca$sdev^2) pdf("fig-15-4.pdf",width=9,height=4.5) par(mfrow=c(1,2),mar=c(4,4,2,2)) plot(res.pca$sco[,1:2], xlab=paste("PC1 (",round(expvar[1]*100,1),"% explained)",sep=""), ylab=paste("PC2 (",round(expvar[2]*100,1),"% explained)",sep=""),cex=1.1,cex.lab=1.2, pch=grp) text(1.5,4.5,"Lithologies") plot(res.pca$sco[,1:2], xlab=paste("PC1 (",round(expvar[1]*100,1),"% explained)",sep=""), ylab=paste("PC2 (",round(expvar[2]*100,1),"% explained)",sep=""),cex=1.1,cex.lab=1.2, pch=res2cl) text(1.5,4.5,"Mclust, 5 clusters") dev.off()

# Fig. 15.5.: Mclust results in map library(StatDA) data(chorizon) el=c("Ca","Zn","Mg","Cu","Sr","Na") x=chorizon[,el] sel=((chorizon[,"LITO"]==7) | (chorizon[,"LITO"]==9) | (chorizon[,"LITO"]==82)) grp=chorizon[sel,"LITO"] grp[grp==7] <- 3 grp[grp==9] <- 4 grp[grp==82] <- 20 x=scale(log10(x[sel,])) lab=chorizon[sel,"LITO"] res2=Mclust(x,5:6) res2cl=res2$cl res2cl[res2$cl==5] <- 20 X=chorizon[,"XCOO"] Y=chorizon[,"YCOO"] # true representation of x and y axis of map for plot xwid=diff(range(X))/12e4 ywid=diff(range(Y))/12e4 pdf("fig-15-5.pdf",width=2*xwid,height=1*ywid) par(mfrow=c(1,2),mar=c(1.5,1.5,1.5,1.5)) # NEW MAP plot(X,Y,frame.plot=FALSE,xaxt="n",yaxt="n",xlab="",ylab="",type="n") data(kola.background) plotbg(map.col=c("gray","gray","gray","gray"),add.plot=T) points(chorizon[sel,"XCOO"],chorizon[sel,"YCOO"],col=1,pch=grp) # Legend text(752000,7880000,"Lithologies",cex=1) # scalebar scalebar(761309,7373050,861309,7363050,shifttext=-0.5,shiftkm=37e3,sizetext=0.8) # North arrow Northarrow(362602,7818750,362602,7878750,362602,7838750,Alength=0.15,Aangle=15,Alwd=1.3,Tcex=1.6) # NEW MAP plot(X,Y,frame.plot=FALSE,xaxt="n",yaxt="n",xlab="",ylab="",type="n") plotbg(map.col=c("gray","gray","gray","gray"),add.plot=T) points(chorizon[sel,"XCOO"],chorizon[sel,"YCOO"],col=1,pch=res2cl) # Legend text(752000,7880000,"Mclust, 5 clusters",cex=1) # scalebar scalebar(761309,7373050,861309,7363050,shifttext=-0.5,shiftkm=37e3,sizetext=0.8) # North arrow Northarrow(362602,7818750,362602,7878750,362602,7838750,Alength=0.15,Aangle=15,Alwd=1.3,Tcex=1.6) dev.off()