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

1、a,1,常用统计方法用R实现,a,2,描述性统计,位置的度量: 均值、顺序统计量、中位数、百分位数。 均值计算：,若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的和。,a,3,顺序统计量,将n个数据(观测值)按从小到大的顺序排列后,称其为顺序统计量. 函数sor

2、t(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),a,4,中位数,中位数描述数据中心位置的数字特征.大体上比中位数大或小的数据个数为整个数据的一半.对于对称分布的数据，均值与中位数比较接近;对于偏态分布的数据，均值与中位数不同.中位数的又一显著特点是不受异常值的影响，具有稳健性，因此它是数据分析中相当重要的统计量. 在R软件中，函数median()给观测量的中位数.如 x-c(75,64,47.4,66.9,62.2,6

3、2.2,58.7,63.5) median(x) median(x,na.rm=TRUE) #若数据中有缺失值,a,5,百分位数,百分位数(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),a,6,分散程度的

4、度量,表示数据分散(或变异)程度的特征量有方差、标准差、极差、四分位极差、变异系数和标准误等. 在R软件中，用var()和sd()计算方差、标准差： var(x, na.rm=FALSE,) sd(x,na.rm=FALSE),a,7,变异系数、平方和,对于变异系数、校正平方和、未校正平方和等指标，需要编写简单的程序. 变异系数CV计算: cv-100*sd(x)/mean(x);cv 校正平方和CSS: css-sum(x-mean(x)2);css 未校正平方和USS: uss-sum(x2);uss,a,8,极差与标准误,样本极差(记为R)的计算: R=max(x)-min(x) 样本上

5、、下四分位数之差称为四分位差(或半极差)，记为R1.它也是度量样本分散性的重要数字特征，特别对于具有异常值的数据，它作为分散性具有稳健性，因此在稳健性数据分析中具有重要作用. 半极差计算：R1= quantile(x,0.75)- quantile(x,0.25) 样本标准误(记为sm)定义为s/sqrt(n) 样本标准误计算：sm=sd(x)/sqrt(length(x),a,9,分布形状的度量,偏度系数Kurtosis是刻划数据的对称性指标.关于均值对称的数据其偏度系数为0.右侧更分散的数据偏度系数为正,左侧更分散的数据偏度系数为负. 当数据的总体分布为正态分布时，峰度系数Skewness

6、近似为0;当峰度系数为正时,两侧极端数据较多;当峰度系数为负时,两侧极端数据较少.,a,10,偏度系数Skewness,样本峰度系数sk计算程序 n-length(x ) m-mean(x) s-sd(x) sk-n/(n-1)*(n-2)*sum(x-m)3)/s3 计算公式,a,11,峰度系数Kurtosis计算,样本峰度系数ku计算程序 n-length(x m-mean(x) s-sd(x) ku-(n*(n+1)/(n-1)*(n-2)*(n-3)*sum(x-m)4)/s4 -(3*(n-1)2)/(n-2)*(n-3) 计算公式,a,12,相关分析,R软件采用用cov()函数计算

7、协方差或协方差阵，用cor()函数计算相关矩阵(相关系数)。 函数cov()和cor()的使用格式为： cov(x,y=NULL,use=all.obs“,method=c(pearson,kendall,spearman) cor(x,y=NULL,use=all.obs“,method=c(pearson,kendall,spearman) 其中x是数值型向量、矩阵或数据框.y是空值(NULL，缺省值)、向量、矩阵或数据框，但需要与x的维数相一致. 与cov和cor有关的函数还有: cov.wt-计算加权协方差(加权协方差矩阵);cor.test-计算相关性检验.,a,13,相关分析示例,

8、例为了解某种橡胶的性能，今抽取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),a,14,回归分析,案例：根据经验，在人的身高相等的情况下，血压的收缩压Y与体

9、重X1(千克)、年龄X2(岁数)有关.现收集了13个男子的数据，见表 .试建立Y关于X1，X2的线性回归方程. 估计出Y=b0+b1X1+b2X2 F检验: H0: b1=b2=0. T检验: H0: bj=0 j=0,1,2,a,15,求解程序,blood-data.frame( 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,40,40,20),Y=c(120,141,124,126,117,125,123,125,132,12

10、3,132,155,147) ) #建立数据框 lm.sol-lm(YX1+X2,data=blood) #进行回归分析 summary(lm.sol) #汇总分析结果 Y=-62.96+2.136X1+0.4002X2. 预测：X=(80, 40)时，相应Y的概率为0. 95的预测区间. new-data.frame(X1=c(80,75),X2=c(40,38) lm.pred|t|) (Intercept) -62.96336 16.99976 -3.704 0.004083 * X1 2.13656 0.17534 12.185 2.53e-07 * X2 0.40022 0.0832

11、1 4.810 0.000713 * - Signif. codes: 0 * 0.001 * 0.01 * 0.05 . 0.1 1 Residual standard error: 2.854 on 10 degrees of freedom Multiple R-squared: 0.9461, Adjusted R-squared: 0.9354 F-statistic: 87.84 on 2 and 10 DF, p-value: 4.531e-07 预测结果如下： fit lwr upr 1 123.9699 117.2889 130.6509,a,17,回归诊断,par(mfro

12、w=c(2,2) #设置画图为2x2的格式 plot(lm.sol,which=c(1:4) #模型检验4张图，包括残差图、QQ图和Cook距离图 数据太少，上面诊断结果并不理想。 library(car) #载入程序包Car,vif()函数在其内 round(vif(lm.sol),2) #计算模型的方差膨胀因子，用2位小数点的格式展示 各变量的方差膨胀因子情况如下： X1 X2 1.96 1.96 可以看到所有参数估计的VIFj=1/(1-Rj2)值都远远小于10，并且接近1。因此这里我们不用担心多重共线性的问题。,a,18,二项选择模型,当我们考虑多个连续解释变量对某个取0-1值的响应变

13、量的影响时，R中常用probit或logit回归来分析。 probit: -1(PY=1)=0+X logit: logit(PY=1)=log(PY=1/(1-PY=1)=0+X 对二项选择的probit/logit回归,R软件可用glm()处理. fm-glm(formula,family=binomial(link=probit),data=data.frame) fm-glm(formula,family=binomial(link=logit),data=data.frame) 在用glm()函数处理二项选择模型时，公式中响应变量y的输入形式有两种：1、y中第一列为对应自变量的响应次

14、数，第2列是不响应的次数；2、y是只由0、1构成的向量，分别表示对应自变量取值是不响应还是相应。,a,19,二项选择案例1,研究小电流对农场动物的影响.选择了7头牛,6种电击强度0,1,3,4,5毫安.给出每种电击强度70次试验中牛发生响应的总次数.试分析电击对牛的影响。,a,20,案例1的程序,norell-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) glm.sol-glm(Ymatx,family=binom

15、ial,data=norell) # logit回归 #glm.sol-glm(Ymatx,family=binomial(link=probit),data=norell) summary(glm.sol) 预测: pre-predict(glm.sol,data.frame(x=3.5) p-exp(pre)/(1+exp(pre);p d-seq(0,5,len=100) pre-predict(glm.sol,data.frame(x=d) p-exp(pre)/(1+exp(pre) norell$y-norell$success/norell$n plot(norell$x,nor

16、ell$y);lines(d,p),a,21,二项选择案例2,50位急性林巴细胞性白血病病人，在入院治疗时取得了外辕血中的细胞数X1、林巴结浸润等级X2(分为0,1,2,3级);出院后有无巩固治疗X3(1”表示有巩固治疗，0”表示无巩固治疗).并取得病人的生存时间,Y=0表示生存时间在1年以内，Y=1表示生存时间在1年或1年以上.试分析病人生存时间长短的概率与X1,X2,X3的关系.,a,22,案例2的程序,life-data.frame( 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,

17、3.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, 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), 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), 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) glm.sol-glm(YX1+X2+X3,family=binomial,data=life) summary(glm.sol),a,23,定序回归模型,当我们考察多个连续解释变