对于人口模型可以采用Logistic增长函数形式，它考虑了初期的指数增长以及总资源的限制。其函数形式如下。
首先载入car包以便读取数据，然后使用nls函数进行建模，其中theta1、theta2、theta3表示三个待估计参数，start设置了参数初始值，设定trace为真以显示迭代过程。nls函数默认采用Gauss-Newton方法寻找极值，迭代过程中第一列为RSS值，后面三列是各参数估计值。然后用summary返回回归结果。
library(car)
pop.mod1 <- nls(population ~ theta1/(1+exp(-(theta2+theta3*year))),start=list(theta1 = 400, theta2 = -49, theta3 = 0.025), data=USPop, trace=T)
summary(pop.mod)
在上面的回归过程中我们直接指定参数初始值，另一种方法是采用搜索策略，首先确定参数取值范围，然后利用nls2包的暴力方法来得到最优参数。但这种方法相当费时。
还有一种更为简便的方法就是采用内置自启动模型(self-starting Models)，此时我们只需要指定函数形式，而不需要指定参数初始值。本例的logistic函数所对应的selfstarting函数名为SSlogis
pop.mod2 <- nls(population ~ SSlogis(year,phi1,phi2,phi3),data=USPop)
二、判断拟合效果
非线性回归模型建立后需要判断拟合效果，因为有时候参数最优化过程会捕捉到局部极值点而非全局极值点。最直观的方法是在原始数据点上绘制拟合曲线。
library(ggplot2)
p <- ggplot(USPop,aes(year, population))
主成分分析法是数据挖掘中常用的一种降维算法,是Pearson在1901年提出的,再后来由hotelling在1933年加以发展提出的一种多变量的统计方法，其最主要的用途在于“降维”，通过析取主成分显出的最大的个别差异,也可以用来削减回归分析和聚类分析中变量的数目，与因子分析类似。所谓降维，就是把具有相关性的变量数目减少，用较少的变量来取代原先变量。如果原始变量互相正交，即没有相关性，则主成分分析没有效果。
在生物信息学的实际应用情况中，通常是得到了成百上千个基因的信息，这些基因相互之间会有影响，通过主成分分析后，得到有限的几个主成分就可以代表它们的基因了。也就是所谓的降维。
R语言有非常多的途径做主成分分析，比如自带的princomp()和psych包的principal()函数，还有gmodels包的fast.prcomp函数。
2 拆解主成分分析步骤
实际应用时我们通常会选择主成分分析函数，直接把input数据一步分析到位，只需要看懂输出结果即可。但是为了加深理解，这里一步步拆解主成分分析步骤，讲解原理。
2.1 测试数据
数据集USJudgeRatings包含了律师对美国高等法院法官的评分。数据框包含43个样本，12个变量!
下面简单看一看这12个变量是什么，以及它们的相关性。
library(knitr)
kable(head(USJudgeRatings))
这12个变量的介绍如下：
[,1]CONTNumber of contacts of lawyer with judge.
[,2]INTGJudicial integrity.司法诚实性
[,3]DMNRDemeanor.风度；举止；行为
[,4]DILGDiligence.勤奋，勤勉；注意的程度
[,5]CFMGCase flow managing.
[,6]DECIPrompt decisions.
[,7]PREPPreparation for trial.
[,8]FAMIFamiliarity with law.
[,9]ORALSound oral rulings.
[,10] WRITSound written rulings.
[,11] PHYSPhysical ability.
[,12] RTENWorthy of retention.
这些是专业领域的用词，大家可以先不用纠结具体细节。
2.2 为什么要做主成分分析
不管三七二十一就直接套用统计方法都是耍流氓，做主成分分析并不是拍脑袋决定的。 在这个例子里面，我们拿到了这43个法官的12个信息，就可以通过这12个指标来对法官进行分类，但也许实际情况下收集其他法官的12个信息比较麻烦，所以我们希望只收集三五个信息即可，然后也想达到比较好的分类效果。或者至少也得剔除几个指标吧，这个时候主成分分析就能派上用场啦。这12个变量能得到12个主成分，如果前两个主成分可以揭示85%以上的变异度，也就是说我们可以用两个主成分来代替那12个指标。
在生物信息学领域，比如我们测了1000个病人的2万个基因的表达矩阵，同时也有他们的健康状态信息。那么我们想仔细研究这些数据，想得到基因表达与健康状态的某种关系。这样我就可以对其余几十亿的人检测基因表达来预测其健康状态。如果我们进行了主成分分析，就可以选择解释度比较高的主成分对应的基因，可能就几十上百个而已，大幅度的降低广泛的基因检测成本。
2.3 step1:数据标准化(中心化)
dat_scale=scale(USJudgeRatings,scale=F)
options(digits=4, scipen=4)
kable(head(dat_scale))
scale()是对数据中心化的函数，当参数scale=F时，表示将数据按列减去平均值，scale=T表示按列进行标准化，公式为（x-mean(x)）/sd(x)，本例采用中心化。
2.4 step2:求相关系数矩阵
dat_cor=cor(dat_scale)
options(digits=4, scipen=4)
kable(dat_cor)
从相关系数看，只有 CONT 变量跟其它变量没有相关性。
当然，这样的矩阵不易看清楚规律，很明显，画个热图就一眼看出来了。
2.5 step3:计算特征值和特征向量
利用eigen函数计算相关系数矩阵的特征值和特征向量。这个是主成分分析法的精髓。
dat_eigen=eigen(dat_cor)
dat_var=dat_eigen$values ## 相关系数矩阵的特征值
options(digits=4, scipen=4)
dat_var
## [1] 10.133504 1.104147 0.332902 0.253847 0.084453 0.037286 0.019683
## [8] 0.015415 0.007833 0.005612 0.003258 0.002060
pca_var=dat_var/sum(dat_var)
pca_var
## [1] 0.8444586 0.0920122 0.0277418 0.0211539 0.0070377 0.0031072 0.0016402
## [8] 0.0012846 0.0006528 0.0004676 0.0002715 0.0001717
pca_cvar=cumsum(dat_var)/sum(dat_var)
pca_cvar
## [1] 0.8445 0.9365 0.9642 0.9854 0.9924 0.9955 0.9972 0.9984 0.9991 0.9996
## [11] 0.9998 1.0000
可以看出，PC1(84.4%)和PC2(9.2%)共可以解释这12个变量的93.6的程度，除了CONT外的其他的11个变量与PC1都有较好的相关性，所以PC1与这11个变量基本斜交，而CONT不能被PC1表示，所以基本与PC1正交垂直，而PC2与CONT基本平行，表示其基本可以表示CONT。
2.6 step4:崖低碎石图和累积贡献图
library(ggplot2)
p=ggplot(,aes(x=1:12,y=pca_var))
p1=ggplot(,aes(x=1:12,y=pca_cvar))
p+geom_point(pch=2,lwd=3,col=2)+geom_line(col=2,lwd=1.2)
打开APP查看高清大图
img
p1+geom_point(pch=2,lwd=3,col=2)+geom_line(col=2,lwd=1.2)
打开APP查看高清大图
img
崖低碎石图（scree plot）即贡献率图，是希望图形一开始很陡峭，如悬崖一般，而剩下的数值都很小，如崖底的碎石一样。
崖低碎石图和累积贡献率图是对主成分贡献率和累积贡献率的一种直观表示，用以作为选择主成分个数的参考。本例中第一个主成分解释总变异的84.4%,可以只选择第一个主成分，但第二主成分也不小，因此选择前两个主成分。
主成分的个数选择没有一定之规，需按实际情况具体分析，一般要求累积贡献率大于85%或特征值大于1.
但是在实际的生物信息学应用中，通常达不到这个要求。
2.7 step5:主成分载荷
主成分载荷表示各个主成分与原始变量的相关系数。
pca_vect= dat_eigen$vector ## 相关系数矩阵的特征向量
loadings=sweep(pca_vect,2,sqrt(pca_var),"*")
rownames(loadings)=colnames(USJudgeRatings)
options(digits=4, scipen=4)
kable(loadings[,1:2])
CONT 0.0028 0.2830
INTG -0.2652 -0.0552
DMNR -0.2636 -0.0599
DILG -0.2797 0.0110
CFMG -0.2780 0.0511
DECI -0.2774 0.0388
PREP -0.2843 0.0098
FAMI -0.2818 -0.0004
ORAL -0.2874 -0.0011
WRIT -0.2858 -0.0095
PHYS -0.2580 0.0270
RTEN -0.2847 -0.0119
结果表明，CONT指标跟其它指标表现完全不一样，第一个主成分很明显跟除了CONT之外的所有其它指标负相关，而第二个主成分则主要取决于CONT指标。
2.8 step6:主成分得分计算和图示
将中心化的变量dat_scale乘以特征向量矩阵即得到每个观测值的得分。
pca_score=as.matrix(dat_scale)%*%pca_vect
head(pca_score[,1:2])
## [,1][,2]
## AARONSON,L.H. 0.5098 -1.7080
## ALEXANDER,J.M. -2.2676 -0.8508
## ARMENTANO,A.J. -0.2267 -0.2903
## BERDON,R.I.-3.4058 -0.5997
## BRACKEN,J.J.6.5937 0.2478
## BURNS,E.B. -2.3336 -1.3563
将两个主成分以散点图形式展示,看看这些样本被这两个主成分是如何分开的
p=ggplot(,aes(x=pca_score[,1],y=pca_score[,2]))+geom_point(color=USJudgeRatings[,1],pch=USJudgeRatings[,2])
p
打开APP查看高清大图
img
对于主成分分析，不同数据会有不同的分析方法，应具体情况具体分析。
总结一下PCA的算法步骤：
设有m条n维数据。
1）将原始数据按列组成n行m列矩阵X
2）将X的每一行（代表一个属性字段）进行零均值化，即减去这一行的均值
3）求出协方差矩阵
4）求出协方差矩阵的特征值及对应的特征向量
5）将特征向量按对应特征值大小从上到下按行排列成矩阵，取前k行组成矩阵P
6）Y=PX即为降维到k维后的数据
PCA本质上是将方差最大的方向作为主要特征，并且在各个正交方向上将数据“离相关”，也就是让它们在不同正交方向上没有相关性。
PCA也存在一些限制，例如它可以很好的解除线性相关，但是对于高阶相关性就没有办法了，对于存在高阶相关性的数据，可以考虑Kernel PCA，通过Kernel函数将非线性相关转为线性相关，关于这点就不展开讨论了。另外，PCA假设数据各主特征是分布在正交方向上，如果在非正交方向上存在几个方差较大的方向，PCA的效果就大打折扣了。
最后需要说明的是，PCA是一种无参数技术，也就是说面对同样的数据，如果不考虑清洗，谁来做结果都一样，没有主观参数的介入，所以PCA便于通用实现，但是本身无法个性化的优化。
方差分析与回归分析是有联系又不完全相同的分析方法。方差分析主要研究各变量对结果的影响程度的定性关系，从而剔除对结果影响较小的变量，提高试验的效率和精度。而回归分析是研究变量与结果的定量关系，得出相应的数学模式。在回归分析中，需要对各变量对结果影响进行方差分析，以剔除影响不大的变量，提高回归分析的有效性。方差分析(Analysis of Variance，简称ANOVA)，又称“变异数分析”，是R.A.Fisher发明的，用于两个及两个以上样本均数差别的显著性检验。 由于各种因素的影响，研究所得的数据呈现波动状。造成波动的原因可分成两类，一是不可控的随机因素，另一是研究中施加的对结果形成影响的可控因素。方差分析是从观测变量的方差入手，研究诸多控制变量中哪些变量是对观测变量有显著影响的变量。
回归分析是研究各因素对结果影响的一种模拟经验方程的办法，回归分析（regression analysis)是确定两种或两种以上变量间相互依赖的定量关系的一种统计分析方法。运用十分广泛，回归分析按照涉及的变量的多少，分为一元回归和多元回归分析。
回归分析中，会用到方差分析来判断各变量对结果的影响程度，从而确定哪些因素是应该纳入到回归方程中，哪些由于对结果影响的方差小而不应该纳入到回归方程中。
