R统计分析

R的统计分析——基于线性模型授课内容授课目的学习如何应用R软件解决统计问题1、统计模型、方法简介2、应用实例3、实验作业统计模型的线性表示•很多统计模型可以用一个线型模型来表示：•写成矩阵形式为2120,1,,;,,,(0,)pijijinjyxiniidN•在R中模型是一种对象，其表达形式叫做一个公式，我们先举几个例子来看一看。•假定下表中y，x，x0，x1，x2，…是数值型变量，X是矩阵，A，B，C，…是因子。YXy~xy~1+x两个式子都表示y对x的简单一元线型回归。第一个式子带有隐含的截距项，而第二个式子把截距项显式地写了出来(并不代表截距是1)。y~-1+xy~x-1都表示y对x的通过原点的回归，即不带截距项的回归(-1只是表示不带截距项)。log(y)~x1+x2表示log(y)对x1和x2的二元回归，带有隐含的截距项。y~poly(x,2)y~1+x+I(x^2)表示y对x的一元二次多项式回归。第一种形式使用正交多项式，第二种形式直接使用x的各幂次。y~X+poly(x,2)因变量为y的多元回归，模型矩阵包括矩阵X，以及x的二次多项式的各项。y~A一种方式分组的方差分析，指标为y，分组因素为A。在R中~运算符用来定义模型公式。一般的线性模型的公式形式为因变量~+-第一项+-第二项+-第三项…其中因变量可以是向量或矩阵，或者结果为向量或矩阵的表达式。加号+或者减号－，表示在模型中加入一项或去掉一项，第一项前面如果是加号可以省略。公式中的各项可以取为：•一个值为向量或矩阵的表达式，或1。•一个因子•一个“公式表达式”，由“公式运算符”把因子、向量、矩阵连接而成。•每一项定义了要加入模型矩阵或从模型矩阵中删除的若干列。一个1表示一个截距项列，除非显式地删除总是隐含地包括在模型公式中。•注意在函数调用的括号内的表达式按普通四则运算解释。函数I()可以把一个计算表达式封装起来作为模型的一项使用。•注意R的模型表示只给出了因变量和自变量及自变量间的关系，这样只确定了线性模型的模型矩阵，而模型参数向量是隐含的，并没有的模型公式中体现出来。线性回归模型•拟合普通的线性模型的函数为lm()，其简单的用法为：fitted.model-lm(formula,data=data.frame)其中formula为模型公式;data.frame为各变量所在的数据框。•mod1=lm(y~x1+x2,data=production)可以拟合一个y对x1和x2的二元回归（带有隐含的截距项），数据来自数据框production。拟合的结果存入了对象mod1中。lm()的基本显示十分简练：•mod1Call:lm(formula=y~x1+x2,data=production)Coefficients:(Intercept)x1x20.01220332.0094758-0.0005314只显示了调用的公式和参数估计结果。summary(mod1)•使用summary()函数可以给出更详尽的结果，包括Call、Residentials(列出残差的最小值，1/4分位数，中位数，3/4分位数，最大值)、Coefficients(列出了模型的系数估计，假设检验统计量和p值)。•类似这样的提取信息的函数还有很多，详见后页提取信息的通用函数•lm()函数的返回值叫做模型拟合结果对象，本质上是一个列表，有model、coefficients、residuals等成员。lm()的结果显示十分简单，为了获得更多的拟合信息，可以使用对lm类对象有特殊操作的通用函数，这些函数包括：add1coefeffectskappapredictresidualsaliasdeviancefamilylabelsprintsummaryanovadrop1formulaplotproj•下表给出了lm类（拟合模型类）常用的通用函数的简单说明。通用函数返回值或效果anova(对象1，对象2)把一个子模型与原模型比较，生成方差分析表。coefficients(对象)返回回归系数（矩阵）。可简写为coef(对象)。deviance(对象)返回残差平方和，如有权重则加权。formula(对象)返回模型公式。plot(对象)生成两张图，一张是因变量对拟合值的图形，一张是残差绝对值对拟合值的图形。predict(对象,newdata=数据框)predict.gam(对象,newdata=数据框)有了模型拟合结果后对新数据进行预报。指定的新数据必须与建模时用的数据具有相同的变量结构。函数结构为对数据框中每一观测的因变量预报结果（为向量或矩阵）。predict.gam()与predict()作用相同但适用性更广，可应用于lm、glm和gam的拟合结果。比如，当多项式基函数用了正交多项式时，加入了新数据导致正交多项式基函数改变，用predict.gam()函数可以避免由此引起的偏差print(对象)简单显示模型拟合结果。一般不用print()而直接键入对象名来显示。residuals(对象)返回模型残差（矩阵），若有权重则适当加权。可简写为resid(对象)。summary(对象)可显示较详细的模型拟合结果。•方差分析(AnalysisofVariance，简称ANOVA)，又称“变异数分析”或“F检验”，是R.A.Fisher发明的，用于两个及两个以上样本均数差别的显著性检验。•方差分析基本步骤1、提出原假设：H0：各组样本均值无差异；H1：各组样本均值有显著差异2、选择检验统计量：方差分析采用的检验统计量是F统计量，即F值检验。3、计算检验统计量的观测值和概率P值：该步骤的目的就是计算检验统计量的观测值和相应的概率P值。4、给定显著性水平，并作出决策方差分析•R进行方差分析的函数是aov()，格式为aov(公式，data=数据框)，用法与lm()类似，提取信息的各通用函数仍有效。•我们以“不同牌子木板磨损比较的数据”为例。假设veneer数据框保存了该数据：veneer-read.table(veneer.txt,header=T)•首先我们把每个牌子的木板的磨损情况画盒形图并且放在同一页面中，作图如下：plot(Wear~Brand,data=veneer)ACMEAJAXCHAMPTUFFYXTRA2.02.22.42.6BrandWear•这种图可以直观地比较一个变量在多个组的分布，或者比较几个类似的变量。从图中可以看出，AJAX牌子较好，TUFFY较差，其它三个牌子差别不明显。•为了检验牌子这个因素对指标磨损量有无显著影响，只要用aov()函数：•aov.veneer=aov(Wear~Brand,data=veneer)•summary(aov.veneer)DfSumSqMeanSqFvaluePr(F)Brand40.617000.154257.4040.001683**Residuals150.312500.02083---Signif.codes:0'***'0.001'**'0.01'*'0.05'.'0.1''1可见因素是显著的。统计分析实例(一)•Forbes数据•19世纪四五十年代，苏格兰物理学家JamesD.Forbes试图通过水的沸点来估计海拔高度。他在阿尔卑斯山及苏格兰收集数据以后得到一些值，发表在论文中，下例收集了论文中的17个地方的数据，进行分析。分析的内容有：•气压和沸点是如何联系的？•关系是强还是弱？•能够根据温度预测气压？如果能，有效性如何？Forbes理论认为：沸点和气压对数值呈线性关系，但是取log以后值较小，因此放大100倍•将数据整理成tab分隔的txt文档，加载每列标题(F,h,log,100log)，读入R中•forbes-read.table(forbes.txt,header=T,sep='\t');•先考察一下散点图，看看F和100log之间的关系•plot(forbes$F,forbes$100log)•lm.sol-lm(log100~F,data=forbes);•summary(lm.sol)•abline(lm.sol,col=red)•lm.res-residuals(lm.sol);•plot(lm.res);•text(locator(1),abnormalpoint,adj=0);•identify(lm.res)[1]12•lm.updated-lm(log100~F,data=forbes,subset=-12);•summary(lm.updated)统计分析实例(二)•下面我们以那个学生班的情况为例进行一些分析。我们希望了解体重、身高、年龄、性别等变量的基本情况及互相之间的关系。数据输入•我们先把数据读入一个R数据框对象中：cl-read.table(cl.txt,header=T,sep='\t');探索性数据分析（EDA）•首先我们先研究各变量的分布情况，看分布是否接近正态，有无明显的异常值，有没有明显的序列相关，等等。•研究连续型变量的分布，可以使用直方图、盒形图、分布密度估计图和正态概率图。研究离散型变量分布只要画其分布频数条形图即可，分布频数用table函数计算。研究序列相关性可以作时间序列图和自相关函数图。因为这些图经常重复使用，我们把它们定义为函数，在同一页面画出：•eda.shape=function(x){oldpar=par(mfrow=c(2,2),mar=c(2,2,0.2,0.2),mgp=c(1.2,0.2,0));hist(x,main=,xlab=,ylab=);boxplot(x);iqd=summary(x)[5]-summary(x)[2];plot(density(x,width=2*iqd),xlab=x,ylab=,type=l,main=);qqnorm(x,main=,xlab=,ylab=);qqline(x);par(oldpar);invisible()}•eda.ts=function(x){oldpar-par(mfrow=c(2,1),mar=c(2,2,1,0.2),mgp=c(1.2,0.2,0));#生成指定格式的图形，并将原来的格式存在oldpar里面plot.ts(x,main=,xlab=);acf(x,main=,xlab=);par(oldpar);#在执行完绘图命令以后，图形格式复原invisible()}•函数中最后的invisible()表示在命令行调用此函数时不要显示任何返回值。•变量iqd计算的是函数的四分位间距。函数density用来作核密度曲线估计，其width参数为核估计的参数。•plot.ts是对时间序列的对象作图。•acf是画出时间序列的自相关系数和自协方差系数；•attach(cl)•summary(cl)•eda.shape(Age)•eda.shape(Height)•eda.shape(Weight)•tab.sex-table(Sex)•barplot(tab.sex)eda.shape(Age)1112131415160123456711121314151681012141618200.000.050.100.15-2-1012111213141516eda.shape(Height)50556065707501234565560657040506070800.000.010.020.030.040.050.06-2-101255606570eda.shape(Weight)406080100120140160024686080100120140501001500.0000.0050.0100.015-2-10126080100120140tab.sex-table(Sex)barplot(tab.sex)FM0246810•因为数据是不同个体的观测所以不可能有序列相关