主成分分析

student.pr <- princomp(student, cor = TRUE)

 summary(student.pr, loadings=TRUE)

predict(student.pr)

screeplot(student.pr,type="lines")

R<-matrix(0, nrow=16, ncol=16, dimnames=list(names, names))
for (i in 1:16){
for (j in 1:i){
R[i,j]<-x[(i-1)*i/2+j]; R[j,i]<-R[i,j]
}

}

pr<-princomp(covmat=R); load<-loadings(pr)

plot(load[,1:2]); text(load[,1], load[,2], adj=c(-0.4, 0.3))

 conomy.pr<-princomp(~x1+x2+x3, data=conomy, cor=T)

 summary(conomy.pr, loadings=TRUE)

pre<-predict(conomy.pr)
 conomy$z1<-pre[,1]; conomy$z2<-pre[,2]
 lm.sol<-lm(y~z1+z2, data=conomy)

 summary(lm.sol)

> beta<-coef(lm.sol); A<-loadings(conomy.pr)
> x.bar<-conomy.pr$center; x.sd<-conomy.pr$scale
> coef<-(beta[2]*A[,1]+ beta[3]*A[,2])/x.sd

> beta0 <- beta[1]- sum(x.bar * coef)

> c(beta0, coef)

  library("car")
      wine <- read.table("http://archive.ics.uci.edu/ml/machine-learning-databases/wine/wine.data",
                         sep=",")
      
      standardisedconcentrations <- as.data.frame(scale(wine[2:14])) # standardise the variables
      wine.pca <- prcomp(standardisedconcentrations)
      summary(wine.pca)
      wine.pca$sdev
      sum((wine.pca$sdev)^2)
      screeplot(wine.pca, type="lines")
      (wine.pca$sdev)^2

 a=eigen(cov(scale(wine)))
      cca=(a$va)/sum(a$va)
      ca=cumsum(cca)
      cca
      ca

      a$va

   library(FactorMineR)

w=read.csv("LA.Neighborhoods.csv")

      u=w[,-c(1,12:15)]
      res.pca=PCA(u,quanti.sup = 7:10)
      plot(res.pca)
      data(glass,package="anacor")
      names(glass)=tolower(names(glass))
      res.ca=CA(glass)

      res.ca=CA(glass,col.sup=5,row.sup=6)

w$density=w$Population/w$Area
      u=w[,-c(12:15)]
      a=eigen(cor(scale(u[-1])))
      (cca=(a$va)/sum(a$va))
     (ca=cumsum(cca))
      par(mfrow=c(1,2))
      plot(1:11,a$va,type="o",pch=17,col=4,main="scree plot",xlab="component number",ylab="eigen value")
      plot(1:11,ca,type="o",pch=17,col=4,main="cumulative contribution",xlab="component number",ylab="cumulative contribution")
      par(mfrow=c(1,1))
      
      (loadings=sweep(a$vec,2,sqrt(a$value),"*"))->b
      par(mfrow=c(1,2))
      plot(b[,1:2],type="n",xlab="component 1(42%",ylab="component 2(17%)",xlim=c(-1.2,1.2),ylim=c(-1.2,1.2),main="loadings")
      text(b[,1],b[,2],names(u[,-1]))
      abline(h=0)
      abline(v=0)
      
      plot(b[,3:4],type="n",xlab="component 3(11%)",ylab="component 4(9%)",xlim=c(-1.2,1.2),ylim=c(-1.2,1.2),main="loading")
      text(b[,3],b[,4],names(u[,-1]))
      abline(h=0)
      abline(v=0)

      par(mfrow=c(1,1))

par(mfrow=c(1,2))
      sc=as.matrix(scale(u[,-1]))%*%a$ve
      
      plot(sc[,1],sc[,2],type="n",ylim=c(-6,6),xlim=c(-7,6),main="sample principle component",xlab="component 1",ylab="component 2")
      text(sc[,1],sc[,2],u[,1],cex=.4)
      abline(v=0,col=2)
      abline(h=0,col=2)
      plot(sc[,3],sc[,4],type="n",ylim=c(-6,6),xlim=c(-7,6),main="sample principal component ",xlab="component 3",ylab="component 4")
      text(sc[,3],sc[,4],u[,1],cex=.4)
      abline(v=0,col=2)
      abline(h=0,col=2)

      par(mfrow=c(1,1))

pgdata <- read.table("c:\\temp\\pgfull.txt",header=T)

names(pgdata)

pgd <- pgdata[,1:54]

model <- prcomp(pgd,scale=TRUE)

summary(model)

plot(model,main="",col="green")

biplot(model)

yv <- predict(model)[,1]
yv2 <- predict(model)[,2]
windows(7,4)
par(mfrow=c(1,2))
plot(pgdata$hay,yv,pch=16,xlab="biomass",ylab="PC 1",col="red")
plot(pgdata$pH,yv2,pch=16,xlab="soil pH",ylab="PC 2",col="blue")