library(StatDA)
data(moss)
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=log10(moss[,sel])
var=c("Ni","Cu","Mg","Rb","Mn")
x=x[,var]
set.seed(100)
x.mcd=covMcd(x,cor=TRUE)
x.rsc=scale(x,x.mcd$cent,sqrt(diag(x.mcd$cov)))
res0=pfa(x.rsc,factors=2,covmat=x.mcd,scores="regression",rotation="none",
maxit=0,start=rep(0,ncol(x.rsc)))
res1=pfa(x.rsc,factors=2,covmat=x.mcd,scores="regression",rotation="varimax",
maxit=0,start=rep(0,ncol(x.rsc)))
pdf("fig-14-8.pdf",width=10,height=5)
par(mfcol=c(1,2),mar=c(4,4,4,2))
par(pty="s")
biplot(res0$sco[,1:2],res0$loa[,1:2],
xlab="PC1 original",ylab="PC2 original",col=c(gray(0.6),1),xlabs=rep("+",nrow(x)),cex=0.8)
biplot(res1$sco[,1:2],res1$loa[,1:2],
xlab="PC1 rotated",ylab="PC2 rotated",col=c(gray(0.6),1),xlabs=rep("+",nrow(x)),cex=0.8)
loa0=res0$loa
loa1=res1$loa
rotm=solve(t(loa0)%*%loa0)%*%t(loa0)%*%loa1
abline(c(0,rotm[1,2]/rotm[1,1]),lty=2)
abline(c(0,rotm[2,2]/rotm[2,1]),lty=2)
dev.off()
library(StatDA)
data(moss)
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=log10(moss[,sel])
var=c("Ni","Cu","Mg","Rb","Mn","Th")
x=x[,var]
set.seed(100)
x.mcd=covMcd(x,cor=TRUE)
x.rsc=scale(x,x.mcd$cent,sqrt(diag(x.mcd$cov)))
#robust PCA
res0=pfa(x.rsc,factors=2,covmat=x.mcd,scores="regression",rotation="varimax",
maxit=0,start=rep(0,ncol(x.rsc)))
# robust PFA
res1=pfa(x.rsc,factors=2,covmat=x.mcd,scores="regression",rotation="varimax")
pdf("fig-14-9.pdf",width=10,height=5)
par(mfcol=c(1,2),mar=c(4,4,4,2))
par(pty="s")
biplot(res0$sco[,1:2],res0$loa[,1:2],
xlab="Robust PC1",ylab="Robust PC2",col=c(gray(0.6),1),xlabs=rep("+",nrow(x)),cex=0.8)
biplot(res1$sco[,1:2],res1$loa[,1:2],
xlab="Robust F1",ylab="Robust F2",col=c(gray(0.6),1),xlabs=rep("+",nrow(x)),cex=0.8)
dev.off()