14.2.3 Mapping the Principal Components

# Fig 14.5: Mapping the scores of log-transf. and log-centred PCA library(StatDA) data(moss) data(kola.background) X=moss[,"XCOO"] Y=moss[,"YCOO"] sel=c("Ag","Al","As","B","Ba","Bi","Ca","Cd","Co","Cr","Cu","Fe","Hg","K","Mg", "Mn","Mo","Na","Ni","P","Pb","Rb","S","Sb","Si","Sr","Th","Tl","U","V","Zn") x=moss[,sel] # Closure problem with log-centring transformation xgeom=10^apply(log10(x),1,mean) x2=x/xgeom x2.obj=log10(x2) res=princomp(x2.obj,cor=TRUE) # true representation of x and y axis of map for plot xwid=diff(range(X))/12e4 ywid=diff(range(Y))/12e4 pdf("fig-14-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: el=res$sco[,1] # generate plot with background 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) SymbLegend(X,Y,el,type="percentile",qutiles<-c(0,0.05,0.25,0.75,0.95,1),symbtype="EDA",symbmagn=0.8, leg.position="topright",leg.title="PC1 scores",leg.title.cex=0.72,leg.round=1,leg.wid=3,leg.just="right") # Percentile legend text(min(X)+diff(range(X))*4/7.1,max(Y),paste(qutiles*100,collapse=","),cex=0.7) text(min(X)+diff(range(X))*4/7,max(Y)-diff(range(Y))/25,"Percentiles",cex=0.75) # 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: el=res$sco[,2] # generate plot with background 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) SymbLegend(X,Y,el,type="percentile",qutiles<-c(0,0.05,0.25,0.75,0.95,1),symbtype="EDA",symbmagn=0.8, leg.position="topright",leg.title="PC2 scores",leg.title.cex=0.72,leg.round=1,leg.wid=3,leg.just="right") # Percentile legend text(min(X)+diff(range(X))*4/7.1,max(Y),paste(qutiles*100,collapse=","),cex=0.7) text(min(X)+diff(range(X))*4/7,max(Y)-diff(range(Y))/25,"Percentiles",cex=0.75) # 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()

# Fig 14.6: Biplot for log-centred PCA 1-3, map PC3 library(StatDA) data(moss) data(kola.background) X=moss[,"XCOO"] Y=moss[,"YCOO"] sel=c("Ag","Al","As","B","Ba","Bi","Ca","Cd","Co","Cr","Cu","Fe","Hg","K","Mg", "Mn","Mo","Na","Ni","P","Pb","Rb","S","Sb","Si","Sr","Th","Tl","U","V","Zn") x=moss[,sel] # Closure problem with log-centring transformation xgeom=10^apply(log10(x),1,mean) x2=x/xgeom x2.obj=log10(x2) res2=princomp(x2.obj,cor=TRUE) # true representation of x and y axis of map for plot xwid=diff(range(X))/12e4 ywid=diff(range(Y))/12e4 pdf("fig-14-6.pdf",width=2*xwid,height=1*ywid) par(mfcol=c(1,2),mar=c(4,4,4,2)) biplot(res2,xlab="PC1 (clr-transformed)",ylab="PC3 (clr-transformed)",col=c(gray(0.6),1),xlabs=rep("+",nrow(x)),cex=0.8, choices=c(1,3)) par(mar=c(1.5,1.5,1.5,1.5)) ### NEW MAP: el=res2$sco[,3] # generate plot with background 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) SymbLegend(X,Y,el,type="percentile",qutiles<-c(0,0.05,0.25,0.75,0.95,1),symbtype="EDA",symbmagn=0.8, leg.position="topright",leg.title="PC3 scores",leg.title.cex=0.72,leg.round=1,leg.wid=3,leg.just="right") # Percentile legend text(min(X)+diff(range(X))*4/7.1,max(Y),paste(qutiles*100,collapse=","),cex=0.7) text(min(X)+diff(range(X))*4/7,max(Y)-diff(range(Y))/25,"Percentiles",cex=0.75) # 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()