R语言常用统计方法实现

常用统计方法用R实现描述性统计•位置的度量:均值、顺序统计量、中位数、百分位数。•均值计算：niixnx11•若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(x^2);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)/s^3•计算公式niixxsnnnSk133)()2)(1(峰度系数Kurtosis计算•样本峰度系数ku计算程序n-length(xm-mean(x)s-sd(x)ku-((n*(n+1))/((n-1)*(n-2)*(n-3))*sum((x-m)^4)/s^4-(3*(n-1)^2)/((n-2)*(n-3)))•计算公式)3)(2()1(3)()3)(2)(1()1(2144nnnxxsnnnnnKunii相关分析•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+b2X2•F检验: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(Y~X1+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-predict(lm.sol,new,interval=prediction,level=0.95)•lm.pred回归分析的汇总结果Call:lm(formula=Y~X1+X2,data=blood)Residuals:Min1QMedian3QMax-4.0404-1.01830.46400.69084.3274Coefficients:EstimateStd.ErrortvaluePr(|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.1‘’1Residualstandarderror: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(P{Y=1})=0+X•logit:logit(P{Y=1})=log(P{Y=1}/(1-P{Y=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(Ymat~x,family=binomial,data=norell)#logit回归#glm.sol-glm(Ymat~x,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,2