R语言常用统计方法实现PPT课件.ppt

常用统计方法用R实现1.描述性统计位置的度量:均值、顺序统计量、中位数、百分位数。均值计算：若x是向量、矩阵，则mean(x)返回其全部元素均值。若要返回数组某一维的均值：apply(x,dim,mean);dim=1计算行均值，dim=2计算列均值。若x是数据框，则mean(x)返回各列的均值Mean的一般用法：mean(x,trim=0,na.rm=FALSE)trim指定去掉x两端数的比例；na.rm=TRUE允许有缺失值。类似有sum(x)函数可求x的和。2.顺序统计量将n个数据(观测值)按从小到大的顺序排列后,称其为顺序统计量.函数sort(x)给出了样本x的顺序统计量order()给出排序后的下标rank()给出了样本x的秩次统计量x-c(75,64,47.4,66.9,62.2,62.2,58.7,63.5)sort(x)order(x)3.中位数中位数描述数据中心位置的数字特征.大体上比中位数大或小的数据个数为整个数据的一半.对于对称分布的数据，均值与中位数比较接近;对于偏态分布的数据，均值与中位数不同.中位数的又一显著特点是不受异常值的影响，具有稳健性，因此它是数据分析中相当重要的统计量.在R软件中，函数median()给观测量的中位数.如x-c(75,64,47.4,66.9,62.2,62.2,58.7,63.5)median(x)median(x,na.rm=TRUE)#若数据中有缺失值4.百分位数百分位数(percentile)是中位数的推广.将数据按从小到大的排列后，0p1，它的p分位点定义为：在R软件中,quantile()函数计算观测量的百分位数.如w-c(75.0,64.0,47.4,66.9,62.2,62.2,58.7,63.5,66.6,64.0,57.0,69.0,56.9,50.0,72.0)quantile(w)一般用法：quantile(x,probs=seq(0,1,0.25),na.rm=FALSE)5.分散程度的度量表示数据分散(或变异)程度的特征量有方差、标准差、极差、四分位极差、变异系数和标准误等.在R软件中，用var()和sd()计算方差、标准差：var(x,na.rm=FALSE,)sd(x,na.rm=FALSE)6.变异系数、平方和对于变异系数、校正平方和、未校正平方和等指标，需要编写简单的程序.变异系数CV计算:cv-100*sd(x)/mean(x);cv 校正平方和CSS:css-sum(x-mean(x)2);css未校正平方和USS:uss-sum(x2);uss7.极差与标准误样本极差(记为R)的计算:R=max(x)-min(x)样本上、下四分位数之差称为四分位差(或半极差)，记为R1.它也是度量样本分散性的重要数字特征，特别对于具有异常值的数据，它作为分散性具有稳健性，因此在稳健性数据分析中具有重要作用.半极差计算：R1=quantile(x,0.75)-quantile(x,0.25)样本标准误(记为sm)定义为s/sqrt(n)样本标准误计算：sm=sd(x)/sqrt(length(x)8.分布形状的度量偏度系数Kurtosis是刻划数据的对称性指标.关于均值对称的数据其偏度系数为0.右侧更分散的数据偏度系数为正,左侧更分散的数据偏度系数为负.当数据的总体分布为正态分布时，峰度系数Skewness近似为0;当峰度系数为正时,两侧极端数据较多;当峰度系数为负时,两侧极端数据较少.9.偏度系数Skewness样本峰度系数sksk计算程序n-length(x)n-length(x)m-mean(x)m-mean(x)s-sd(x)s-sd(x)sksk-n/(n-1)*(n-2)*sum(x-m)3)/s3-n/(n-1)*(n-2)*sum(x-m)3)/s3计算公式10.峰度系数Kurtosis计算样本峰度系数kuku计算程序n-length(xn-length(xm-mean(x)m-mean(x)s-sd(x)s-sd(x)kuku-(n*(n+1)/(n-1)*(n-2)*(n-3)*sum(x-(n*(n+1)/(n-1)*(n-2)*(n-3)*sum(x-m)4)/s4-(3*(n-1)2)/(n-2)*(n-3)m)4)/s4-(3*(n-1)2)/(n-2)*(n-3)计算公式11.相关分析R R软件采用用cov()cov()函数计算协方差或协方差阵，用cor()cor()函数计算相关矩阵(相关系数)。函数cov()cov()和cor()cor()的使用格式为：cov(x,y=NULL,use=all.obs“,method=c(pearson,kecov(x,y=NULL,use=all.obs“,method=c(pearson,kendall,spearman)ndall,spearman)cor(x,y=NULL,use=all.obs“,method=c(pearson,kecor(x,y=NULL,use=all.obs“,method=c(pearson,kendall,spearman)ndall,spearman)其中x x是数值型向量、矩阵或数据框.y.y是空值(NULL(NULL，缺省值)、向量、矩阵或数据框，但需要与x x的维数相一致.与covcov和corcor有关的函数还有:cov.wt-cov.wt-计算加权协方差(加权协方差矩阵);cor.test-);cor.test-计算相关性检验.12.相关分析示例例为了解某种橡胶的性能，今抽取10个样品，每个测量三项指标:硬度、变形和弹性(rubber.txt).试计算样本均值、样本协方差阵和样本相关矩阵.并用Pearson相关性检验确认变量X1,X2,X3是否相关?rubber-read.table(d:/rubber.txt)mean(rubber)cov(rubber)cor(rubber)cor.test(X1+X2,data=rubber)cor.test(X1+X3,data=rubber)cor.test(X2+X3,data=rubber)13.回归分析案例：根据经验，在人的身高相等的情况下，血压的收缩压Y与体重X1(千克)、年龄X2(岁数)有关.现收集了13个男子的数据，见表.试建立Y关于X1，X2的线性回归方程.估计出Y=b0+b1X1+b2X2F检验:H0:b1=b2=0.T检验:H0:bj=0 j=0,1,214.求解程序blood-data.frame(blood-data.frame(X1=c(76.0,91.5,85.5,82.5,79.0,80.5,74.5,79.0,85.0,76.X1=c(76.0,91.5,85.5,82.5,79.0,80.5,74.5,79.0,85.0,76.5,82.0,95.0,92.5),X2=c(50,20,20,30,30,50,60,50,40,55,5,82.0,95.0,92.5),X2=c(50,20,20,30,30,50,60,50,40,55,40,40,20),Y=c(120,141,124,126,117,125,123,125,132,1240,40,20),Y=c(120,141,124,126,117,125,123,125,132,123,132,155,147)#3,132,155,147)#建立数据框lm.sol-lm(YX1+X2,data=blood)lm.sol-lm(YX1+X2,data=blood)#进行回归分析summary(lm.sol)summary(lm.sol)#汇总分析结果Y=-62.96+2.136X1+0.4002X2.Y=-62.96+2.136X1+0.4002X2.预测：X=(80,40)X=(80,40)时，相应Y Y的概率为0.950.95的预测区间.new-data.frame(X1=c(80,75),X2=c(40,38)new-data.frame(X1=c(80,75),X2=c(40,38)lm.pred-lm.pred|t|)Estimate Std.Error t value Pr(|t|)(Intercept)-62.96336 16.99976 -3.704 0.004083*(Intercept)-62.96336 16.99976 -3.704 0.004083*X1 2.13656 0.17534 12.185 2.53e-07*X1 2.13656 0.17534 12.185 2.53e-07*X2 0.40022 0.08321 4.810 0.000713*X2 0.40022 0.08321 4.810 0.000713*-Signif.codes:0*0.001*0.01*0.05.Signif.codes:0*0.001*0.01*0.05.0.1 1 0.1 1 Residual standard error:2.854 on 10 degrees of freedomResidual standard error:2.854 on 10 degrees of freedomMultiple R-squared:0.9461,Adjusted R-squared:0.9354 Multiple R-squared:0.9461,Adjusted R-squared:0.9354 F-statistic:87.84 on 2 and 10 DF,p-value:4.531e-07F-statistic:87.84 on 2 and 10 DF,p-value:4.531e-07预测结果如下：fit fit lwr lwr upr upr1 1 123.9699123.9699 117.2889 117.2889 130.6509130.650916.回归诊断par(mfrow=c(2,2)par(mfrow=c(2,2)#设置画图为2x22x2的格式 plot(lm.sol,which=c(1:4)plot(lm.sol,which=c(1:4)#模型检验4 4张图，包括残差图、QQQQ图和CookCook距离图数据太少，上面诊断结果并不理想。library(car)#library(car)#载入程序包Car,vif()Car,vif()函数在其内round(vif(lm.sol),2)round(vif(lm.sol),2)#计算模型的方差膨胀因子，用2 2位小数点的格式展示各变量的方差膨胀因子情况如下：X1 X2 X1 X2 1.96 1.96 1.96 1.96 可以看到所有参数估计的VIFVIFj j=1/(1-R=1/(1-Rj j2 2)值都远远小于1010，并且接近1 1。因此这里我们不用担心多重共线性的问题。17.二项选择模型 当我们考虑多个连续解释变量对某个取0-10-1值的响应变量的影响时，R R中常用probitprobit或logitlogit回归来分析。probit:probit:-1-1(PY=1)=(PY=1)=0 0+X Xlogit:logit:logit(PY=1)=log(PY=1/(1-PY=1)=logit(PY=1)=log(PY=1/(1-PY=1)=0 0+X X对二项选择的probit/logitprobit/logit回归,R,R软件可用glm()glm()处理.fm-glm(formula,family=binomial(link=probit),data=data.frame)fm-glm(formula,family=binomial(link=probit),data=data.frame)fm-glm(formula,family=binomial(link=logit),data=data.frame)fm-glm(formula,family=binomial(link=logit),data=data.frame)在用glm()glm()函数处理二项选择模型时，公式中响应变量y y的输入形式有两种：1 1、y y中第一列为对应自变量的响应次数，第2 2列是不响应的次数；2 2、y y是只由0 0、1 1构成的向量，分别表示对应自变量取值是不响应还是相应。18.二项选择案例1研究小电流对农场动物的影响.选择了7头牛,6种电击强度0,1,3,4,5毫安.给出每种电击强度70次试验中牛发生响应的总次数.试分析电击对牛的影响。19.案例1的程序norell-norell-data.frame(x=0:5,n=rep(70,6),success=c(0,9,21,47,60,63)data.frame(x=0:5,n=rep(70,6),success=c(0,9,21,47,60,63)norell$Ymat-cbind(norell$success,norell$n-norell$success)norell$Ymat-cbind(norell$success,norell$n-norell$success)glm.sol-glm(Ymatx,family=binomial,data=norell)#logitglm.sol-glm(Ymatx,family=binomial,data=norell)#logit回归#glm.sol-glm(Ymatx,family=binomial(link=probit),data=norell)#glm.sol-glm(Ymatx,family=binomial(link=probit),data=norell)summary(glm.sol)summary(glm.sol)预测:pre-predict(glm.sol,data.frame(x=3.5)pre-predict(glm.sol,data.frame(x=3.5)p-exp(pre)/(1+exp(pre);pp-exp(pre)/(1+exp(pre);pd-seq(0,5,len=100)d-seq(0,5,len=100)pre-predict(glm.sol,data.frame(x=d)pre-predict(glm.sol,data.frame(x=d)p-exp(pre)/(1+exp(pre)p-exp(pre)/(1+exp(pre)norell$y-norell$success/norell$nnorell$y-norell$success/norell$nplot(norell$x,norell$y);lines(d,p)plot(norell$x,norell$y);lines(d,p)20.二项选择案例25050位急性林巴细胞性白血病病人，在入院治疗时取得了外辕血中的细胞数X1X1、林巴结浸润等级X2(X2(分为0,1,2,30,1,2,3级););出院后有无巩固治疗X3(1”X3(1”表示有巩固治疗，0”0”表示无巩固治疗).).并取得病人的生存时间,Y=0,Y=0表示生存时间在1 1年以内，Y=1Y=1表示生存时间在1 1年或1 1年以上.试分析病人生存时间长短的概率与X1,X2,X3X1,X2,X3的关系.21.案例2的程序life-data.frame(life-data.frame(X1=c(2.5,173,119,10,502,4,14.4,2,40,6.6,X1=c(2.5,173,119,10,502,4,14.4,2,40,6.6,21.4,2.8,2.5,6,3.5,62.2,10.8,21.6,2,3.4,21.4,2.8,2.5,6,3.5,62.2,10.8,21.6,2,3.4,5.1,2.4,1.7,1.1,12.8,1.2,3.5,39.7,62.4,2.4,5.1,2.4,1.7,1.1,12.8,1.2,3.5,39.7,62.4,2.4,34.7,28.4,0.9,30.6,5.8,6.1,2.7,4.7,128,35,34.7,28.4,0.9,30.6,5.8,6.1,2.7,4.7,128,35,2,8.5,2,2,4.3,244.8,4,5.1,32,1.4),2,8.5,2,2,4.3,244.8,4,5.1,32,1.4),X2=rep(c(0,2,0,2,0,2,0,2,0,2,0,2,0,2,0,2,0,2,0,2,0,2,0),X2=rep(c(0,2,0,2,0,2,0,2,0,2,0,2,0,2,0,2,0,2,0,2,0,2,0),c(1,4,2,2,1,1,8,1,5,1,5,1,1,1,2,1,1,1,3,1,2,1,4),c(1,4,2,2,1,1,8,1,5,1,5,1,1,1,2,1,1,1,3,1,2,1,4),X3=rep(c(0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1),X3=rep(c(0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1),c(6,1,3,1,3,1,1,5,1,3,7,1,1,3,1,1,2,9),c(6,1,3,1,3,1,1,5,1,3,7,1,1,3,1,1,2,9),Y=rep(c(0,1,0,1),c(15,10,15,10)Y=rep(c(0,1,0,1),c(15,10,15,10)glm.sol-glm(YX1+X2+X3,family=binomial,data=life)glm.sol-glm(YX1+X2+X3,family=binomial,data=life)summary(glm.sol)summary(glm.sol)22.定序回归模型当我们考察多个连续解释变量对某个取有序离散值的响应变量的影响时，适用定序回归，R R中常用MASSMASS包中polrpolr函数来实现。ProbitProbit定序回归:-1-1(PY(PY k)=k)=k k-k kX XLogitLogit定序回归:logit(PY=1)=logit(PY=1)=k k-k kX Xpolr(formula,data,weights,start,.,subset,na.action,contrasts=NULL,Hess=FALSE,model=TRUE,method=c(logistic,probit,cloglog,cauchit)23.定序回归示例24.25.26.