library(xlsx)
data<-read.xlsx("E://code//food.xlsx",sheetIndex = 1,header = T,
                encoding = "UTF-8",)
data<-data[,-1]#选择除去第一行的数据
library(labdsv)

Loading required package: mgcv
Loading required package: nlme
This is mgcv 1.8-15. For overview type 'help("mgcv-package")'.
Loading required package: MASS
Loading required package: cluster

Attaching package: ‘labdsv’

The following object is masked from ‘package:stats’:

    density

Warning message:
package ‘labdsv’ was built under R version 3.3.3 

data.pca=pca(data,dim=4,cor = T)#利用相关系数矩阵计算
    第一步，将pca（mat，cor，dim）中的参数cor设为T，表示利用相关系数矩阵作主成分分析，相当于使用了标准化数据，因此无需做事前的标准化；而dim表示变量综合后的维数，生成4个主成分。第二步，提取分析结果。4个主成分的标准差、方差贡献以及累计方差的贡献率。loadings.pca()提取载荷系数。

> summary(data.pca)
Importance of components:
                            [,1]      [,2]      [,3]       [,4]
Standard deviation     1.7032337 1.4413673 1.1591787 0.91233021
Proportion of Variance 0.3223339 0.2308377 0.1492995 0.09248294
Cumulative Proportion  0.3223339 0.5531716 0.7024711 0.79495405

> loadings.pca(data.pca)

Loadings:
             PC1    PC2    PC3    PC4   
粮食                -0.557              
蔬菜                -0.313  0.630 -0.434
食油                 0.209  0.686  0.475
猪牛羊肉      0.462         0.168 -0.309
家禽          0.522                0.189
蛋类及其制品 -0.364  0.383  0.261       
水产品        0.491  0.179  0.106  0.218
食糠          0.362  0.393 -0.145       
酒                   0.453        -0.627

    第三步，使用函数varplot.pca()绘制4个因子的碎石图以及累计方差图。

op=par(mfrow=c(1,2))#分割图形区域

varplot.pca(data.pca)

http://s1/mw690/006ySAqhzy7adcsIIIo90&690

#princomp函数的主成分分析

data.pr=princomp(data,cor=T)#用相关阵计算

options(digits=4)

summary(data.pr,loadings = T)#loadings=T选项列出主成分对应原始变量的系数

Importance of components:
                       Comp.1 Comp.2 Comp.3  Comp.4  Comp.5  Comp.6  Comp.7
Standard deviation     1.7032 1.4414 1.1592 0.91233 0.87345 0.62696 0.60360
Proportion of Variance 0.3223 0.2308 0.1493 0.09248 0.08477 0.04368 0.04048
Cumulative Proportion  0.3223 0.5532 0.7025 0.79495 0.87972 0.92340 0.96388
                        Comp.8 Comp.9
Standard deviation     0.47261 0.3190
Proportion of Variance 0.02482 0.0113
Cumulative Proportion  0.98870 1.0000

Loadings:
             Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6 Comp.7 Comp.8 Comp.9
粮食                -0.557               -0.568  0.205 -0.252  0.486  0.112
蔬菜                -0.313 -0.630  0.434  0.151  0.296        -0.349 -0.277
食油                 0.209 -0.686 -0.475 -0.124 -0.222 -0.380         0.226
猪牛羊肉     -0.462        -0.168  0.309  0.443 -0.321  0.130  0.509  0.297
家禽         -0.522               -0.189 -0.193  0.259  0.306 -0.416  0.556
蛋类及其制品  0.364  0.383 -0.261                0.477  0.494  0.363  0.198
水产品       -0.491  0.179 -0.106 -0.218 -0.335         0.348  0.173 -0.630
食糠         -0.362  0.393  0.145         0.182  0.581 -0.547  0.121       
酒                   0.453         0.627 -0.504 -0.278 -0.130 -0.182  0.134

    分析结果standard deviation表示主成分的标准差，即相应特征值（方法）的平方根；Proportion of variance表示方差贡献率，cumulative proportion表示累计方差贡献率。第一主成分的F1的贡献率为32.23%，第二主成分F2的贡献率为23.08%，第三主成分F3的贡献率为14.93%，前5个主成分的累计方差贡献率达到97.97%，选取5个主成分，其余舍去。

    写出主成分与原变量的线性关系。。。【并给出合理的解释】

    碎石图中由陡峭变为平坦的转折点即为主成分选择的最佳个数。

screeplot(data.pr,type="line",main = "碎石图")

http://s1/mw690/006ySAqhzy7adh1jeOQ20&690

    函数biplot（）可绘制数据关于主成分的散点图，并自动表名原坐标在主成分下的方向，调用的格式为biplot（x,choices=1:2），x是princomp（）返回的对象，参数choices选择主成分，默认前两个。



biplot(data.pr)

http://s15/mw690/006ySAqhzy7adgYZ0ZM0e&690





    R中自带的分析函数factanal（）采用极大似然估计方法估计因子载荷，适用于大样本量的数据分析。

这里自己写出主成分法的因子分析函数。





#主成分法的因子分析factor.analysis()

factor.analysis=function(x,m){

  p=nrow(x);x.diag=diag(x);sum.rank=sum(x.diag)

  rowname=paste("X",1:p,sep="")#设置行名、列名

  colname=paste("Factor",1:m,sep="")

  #构造因子载荷矩阵A，初值设为0

  A=matrix(0,nrow = p,ncol = m,dimnames = list(rowname,colname))

  eig=eigen(x)#eig包含两个元素，values为特征根，vectors为特征向量

  for(i in 1:m)

    A[,i]=sqrt(eig$values[i])*eig$vectors[,i]#填充矩阵A的值

  var.A=diag(A%*%t(A))#公共因子的方差

  rownamel=c("SS loadings","ProportionVar","Cumulative Var")

  #构造输出结果的矩阵，初值设为0

  result=matrix(0,nrow=3,ncol = m,dimnames = list(rownamel,colname))

  for(i in 1:m){

    result[1,i]=sum(A[,i]^2)#计算各因子的方差（对变量的贡献）

    result[2,i]=result[1,i]/sum.rank#计算方差贡献率

    result[3,i]=sum(result[1,1:i])/sum.rank#计算方差贡献率

  }

  method=c("Principal Component Method")

  #计算输出结果

  list(method=method,loadings=A,var=cbind(commom=var.A,specific=x.diag-var.A),

       result=result)

}

    函数的输入参数x用于因子分析的样本方法矩阵或相关矩阵；m表示因子个数。函数输出一个列表，包括主成分法计算得到的因子载荷、公共因子方差和特殊因子方差，以及因子的方差、贡献率和累计方差贡献率等，，



R=cor(data)#计算相关系数矩阵

options(digits = 3)#显示3位有效数字

factor.analysis(R,5)

$method
[1] "Principal Component Method"

$loadings
       Factor1     Factor2      Factor3     Factor4    Factor5
X1 -0.07719473 -0.80349179 -0.020348673  0.09013283 -0.4961693
X2 -0.13664156 -0.45092557 -0.730792229  0.39597444  0.1316376
X3  0.01075994  0.30139712 -0.795248274 -0.43334715 -0.1083670
X4 -0.78647898 -0.03849984 -0.194411310  0.28158148  0.3866967
X5 -0.88863483 -0.13149475  0.064827566 -0.17236815 -0.1684075
X6  0.61921226  0.55222537 -0.302596192  0.05267647 -0.0742369
X7 -0.83600069  0.25778990 -0.122585660 -0.19907381 -0.2923060
X8 -0.61598943  0.56603896  0.167954242  0.03201468  0.1589624
X9 -0.07921100  0.65359372 -0.002331231  0.57202093 -0.4399285

$var
      commom   specific
X1 0.9062800 0.09371998
X2 0.9301863 0.06981372
X3 0.9229090 0.07709101
X4 0.8866496 0.11335038
X5 0.8692372 0.13076280
X6 0.7882271 0.21177293
X7 0.9054532 0.09454681
X8 0.7543457 0.24565430
X9 0.9542096 0.04579043

$result
                 Factor1   Factor2   Factor3    Factor4    Factor5
SS loadings    2.9010051 2.0775396 1.3436953 0.83234642 0.76291121
ProportionVar  0.3223339 0.2308377 0.1492995 0.09248294 0.08476791
Cumulative Var 0.3223339 0.5531716 0.7024711 0.79495405 0.87972196