《R语言常用统计方法实现》由会员分享，可在线阅读，更多相关《R语言常用统计方法实现（26页珍藏版）》请在金锄头文库上搜索。
1、常用统计方法用R实现描述性统计位置的度量: 均值、顺序统计量、中位数、百分位数。均值计算：若若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允许有缺失值。允许有缺失值。
2、类似有类似有sum(x)函数可求函数可求x的和。的和。 顺序统计量将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) #若数据中有缺失值百分位数百分位数(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
4、(x,probs=seq(0,1,0.25),na.rm=FALSE)分散程度的度量表示数据分散(或变异)程度的特征量有方差、标准差、极差、四分位极差、变异系数和标准误等.在R软件中，用var()和sd()计算方差、标准差： var(x, na.rm=FALSE,) sd(x,na.rm=FALSE)变异系数、平方和对于变异系数、校正平方和、未校正平方和等指标，需要编写简单的程序.变异系数CV计算: cv-100*sd(x)/mean(x);cv 校正平方和CSS: css-sum(x-mean(x)2);css未校正平方和USS: uss-sum(x2);uss极差与标准误样本极差(记为R)
5、的计算: 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)分布形状的度量偏度系数Kurtosis是刻划数据的对称性指标.关于均值对称的数据其偏度系数为0.右侧更分散的数据偏度系数为正,左侧更分散的数据偏度系数为负.当数据的总体分布为正态分布时
6、，峰度系数Skewness近似为0;当峰度系数为正时,两侧极端数据较多;当峰度系数为负时,两侧极端数据较少.偏度系数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计算公式计算公式峰度系数Kurtosis计算样本样本峰度系数度系数kuku计算程序计算程序n-length(xn-length(xm-mean(x)m-mean(x)s-sd(x)s-sd
7、(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)计算公式计算公式相关分析R R软件采用用软件采用用cov()cov()函数计算协方差或协方差阵，用函数计算协方差或协方差阵，用cor()cor()函数计算相关矩阵函数计算相关矩阵( (相关系数相关系数) )。函数函数cov()cov()和和cor()cor()的使用格式为：的使用格式为：cov(x,y=NULL,use=all
8、.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，缺省，缺省值值) )、向量、矩阵或数据框，但需要与、
9、向量、矩阵或数据框，但需要与x x的维数相一致的维数相一致. . 与与covcov和和corcor有关的函数还有有关的函数还有: : cov.wt-cov.wt-计算加权协计算加权协方差方差( (加权协方差矩阵加权协方差矩阵);cor.test-);cor.test-计算相关性检验计算相关性检验. .相关分析示例例为了解某种橡胶的性能，今抽取10个样品，每个测量三项指标:硬度、变形和弹性(rubber.txt).试计算样本均值、样本协方差阵和样本相关矩阵.并用Pearson相关性检验确认变量X1 , X2,X3是否相关?rubber-read.table(d:/rubber.txt)mean(
10、rubber)cov(rubber)cor(rubber)cor.test(X1+X2,data=rubber)cor.test(X1+X3,data=rubber)cor.test(X2+X3,data=rubber)回归分析案例：根据经验，在人的身高相等的情况下，血压的收缩压Y与体重X1(千克)、年龄X2(岁数)有关.现收集了13个男子的数据，见表 .试建立Y关于X1，X2的线性回归方程.估计出Y=b0+b1X1+b2X2F检验: H0: b1=b2=0. T检验: H0: bj=0 j=0,1,2求解程序blood-data.frame( blood-data.frame( X1=c(7
11、6.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,
12、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-d
13、ata.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
14、 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,
15、 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.6509回归诊断par(mfrow=c(2,
16、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) # #计算模型的方差
17、膨胀因子，用计算模型的方差膨胀因子，用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。因此这里我们不用担心多重共线性的问题。因此这里我们不用担心多重共线性的问题。二项选择模型 当我们考虑多个连续解释变量对某个取当我们考虑多个连续解释变量对某个取0-10-1值的响应变量的影响时，值的响应变量的影响时，R R中常用
18、中常用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-g
19、lm(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列是不响应的列是不
20、响应的次数；次数；2 2、y y是只由是只由0 0、1 1构成的向量，分别表示对应自变量取值是不构成的向量，分别表示对应自变量取值是不响应还是相应。响应还是相应。二项选择案例1研究小电流对农场动物的影响.选择了7头牛,6种电击强度0,1,3,4,5毫安.给出每种电击强度70次试验中牛发生响应的总次数.试分析电击对牛的影响。案例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$Y
21、mat-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,fam
22、ily=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.fram
23、e(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)二项选择案例25050位急性林巴细胞性白血病病人，在入院治疗时取得了外位急性林巴细胞性白血病病人，在入院治疗时取得了外辕血中的细胞数辕血中的细胞数X1X1、林巴结浸润等级、林巴结浸润等级X2(X2(分为分为0,1,2,30,1,2,3级
24、级););出院后有无巩固治疗出院后有无巩固治疗X3(1”X3(1”表示有巩固治疗，表示有巩固治疗，0”0”表示表示无巩固治疗无巩固治疗).).并取得病人的生存时间并取得病人的生存时间,Y=0,Y=0表示生存时间在表示生存时间在1 1年以内，年以内，Y=1Y=1表示生存时间在表示生存时间在1 1年或年或1 1年以上年以上. .试分析病人试分析病人生存时间长短的概率与生存时间长短的概率与X1,X2,X3X1,X2,X3的关系的关系. .案例2的程序life-data.frame(life-data.frame(X1=c(2.5,173,119,10,502,4,14.4,2,40,6.6,X1=c
25、(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,
26、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
27、,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(
28、glm.sol)summary(glm.sol)定序回归模型当我们考察多个连续解释变量对某个取有序离散值的响应当我们考察多个连续解释变量对某个取有序离散值的响应变量的影响时，适用定序回归，变量的影响时，适用定序回归，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)定序回归示例
