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

常用统计方法用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允许有缺失值。类似有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),中位数,中位数描述数据中心位置的数字特征.大体上比中位数大或小的数据个数为整个数据的一半.对于对称分布的数据，均值与中位数比较接近;对于偏态分布的数据，均值与中位数不同.中位数的又一显著特点是不受异常值的影响，具有稳健性，因此它是数据分析中相当重要的统计量.在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(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)的计算: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.右侧更分散的数据偏度系数为正,左侧更分散的数据偏度系数为负.当数据的总体分布为正态分布时，峰度系数Skewness近似为0;当峰度系数为正时,两侧极端数据较多;当峰度系数为负时,两侧极端数据较少.,偏度系数Skewness,样本峰度系数sk计算程序n-length(x)m-mean(x)s-sd(x)sk-n/(n-1)*(n-2)*sum(x-m)3)/s3计算公式,峰度系数Kurtosis计算,样本峰度系数ku计算程序n-length(xm-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)计算公式,相关分析,R软件采用用cov()函数计算协方差或协方差阵，用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-计算相关性检验.,相关分析示例,例为了解某种橡胶的性能，今抽取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),回归分析,案例：根据经验，在人的身高相等的情况下，血压的收缩压Y与体重X1(千克)、年龄X2(岁数)有关.现收集了13个男子的数据，见表.试建立Y关于X1，X2的线性回归方程.估计出Y=b0+b1X1+b2X2F检验:H0:b1=b2=0.T检验:H0:bj=0j=0,1,2,求解程序,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,123,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.9633616.99976-3.7040.004083*X12.136560.1753412.1852.53e-07*X20.400220.083214.8100.000713*-Signif.codes:0*0.001*0.01*0.05.0.11Residualstandarderror:2.854on10degreesoffreedomMultipleR-squared:0.9461,AdjustedR-squared:0.9354F-statistic:87.84on2and10DF,p-value:4.531e-07预测结果如下：fitlwrupr1123.9699117.2889130.6509,回归诊断,par(mfrow=c(2,2)#设置画图为2x2的格式plot(lm.sol,which=c(1:4)#模型检验4张图，包括残差图、QQ图和Cook距离图数据太少，上面诊断结果并不理想。library(car)#载入程序包Car,vif()函数在其内round(vif(lm.sol),2)#计算模型的方差膨胀因子，用2位小数点的格式展示各变量的方差膨胀因子情况如下：X1X21.961.96可以看到所有参数估计的VIFj=1/(1-Rj2)值都远远小于10，并且接近1。因此这里我们不用担心多重共线性的问题。,二项选择模型,当我们考虑多个连续解释变量对某个取0-1值的响应变量的影响时，R中常用probit或logit回归来分析。probit:-1(PY=1)=0+Xlogit: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中第一列为对应自变量的响应次数，第2列是不响应的次数；2、y是只由0、1构成的向量，分别表示对应自变量取值是不响应还是相应。,二项选择案例1,研究小电流对农场动物的影响.选择了7头牛,6种电击强度0,1,3,4,5毫安.给出每种电击强度70次试验中牛发生响应的总次数.试分析电击对牛的影响。,案例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=binomial,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);pd-seq(0,5,len=100)pre-predict(glm.sol,data.frame(x=d)p-exp(pre)/(1+exp(pre)norell$y-norell$success/norell$nplot(norell$x,norell$y);lines(d,p),二项选择案例2,50位急性林巴细胞性白血病病人，在入院治疗时取得了外辕血中的细胞数X1、林巴结浸润等级X2(分为0,1,2,3级);出院后有无巩固治疗X3(1”表示有巩固治疗，0”表示无巩固治疗).并取得病人的生存时间,Y=0表示生存时间在1年以内，Y=1表示生存时间在1年或1年以上.试分析病人生存时间长短的概率与X1,X2,X3的关系.,案例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,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