R语言常用统计方法实现

回归分析的汇总结果
Call: lm(formula = Y ~ X1 + X2, data = blood) Residuals: Min 1Q Median 3Q Max -4.0404 -1.0183 0.4640 0.6908 4.3274 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -62.96336 16.99976 -3.704 0.004083 ** X1 2.13656 0.17534 12.185 2.53e-07 *** X2 0.40022 0.08321 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
相关分析示例
• 例为了解某种橡胶的性能，今抽取10个样品，每个 测量三项指标:硬度、变形和弹性(rubber.txt).试计算 样本均值、样本协方差阵和样本相关矩阵.并用 Pearson相关性检验确认变量X1 , X2,X3是否相关?
• rubber<-read.table("d:/rubber.txt")
极差与标准误
• 样本极差(记为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;当峰度系数为正时,两侧 极端数据较多;当峰度系数为负时,两侧极端 数据较少.
二项选择案例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));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,norell$y);lines(d,p)
常用统计方法用R实现
描述性统计
• 位置的度量: 均值、顺序统计量、中位数、百分位数。
• 均值计算：
x
1
n
n
xi
i 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的和。
• • • • • • 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=0 j=0,1,2
二项选择模型
• 当我们考虑多个连续解释变量对某个取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构成的向量，分别表示对应自变量取值是 不响应还是相应。
百分位数
• 百分位数(percentile)是中位数的推广.将数据按从 小到大的排列后，0<p<1，它的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)
中位数
• 中位数描述数据中心位置的数字特征.大体上比中 位数大或小的数据个数为整个数据的一半.对于对 称分布的数据，均值与中位数比较接近;对于偏态 分布的数据，均值与中位数不同.中位数的又一显 著特点是不受异常值的影响，具有稳健性，因此 它是数据分析中相当重要的统计量. • 在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) #若数据中有缺失值
回归诊断
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位小数点的格式展示 • 各变量的方差膨胀因子情况如下： X1 X2 1.96 1.96 • 可以看到所有参数估计的VIFj=1/(1-Rj2)值都远远小于10， 并且接近1。因此这里我们不用担心多重共线性的问题。
求解程序
• blood<data.frame( X1=c0,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
4

n
( xi x )
4
3 ( n 1)
2
i 1
( 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---计算相关性检验.
顺序统计量
• 将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)
• 样本峰度系数ku计算程序 n<-length(x m<-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))) • 计算公式
Ku n ( n 1) ( n 1)( n 2 )( n 3 ) s
变异系数、平方和
• 对于变异系数、校正平方和、未校正平方和等指 标，需要编写简单的程序. • 变异系数CV计算: cv<-100*sd(x)/mean(x);cv • 校正平方和CSS: css<-sum((x-mean(x))^2);css • 未校正平方和USS: uss<-sum(x^2);uss
偏度系数Skewness
• 样本峰度系数sk计算程序 n<-length(x ) m<-mean(x) s<-sd(x) • sk<-n/((n-1)*(n-2))*sum((x-m)^3)/s^3 • 计算公式
Sk n ( n 1)( n 2 ) s
3
(x
i 1
n
i
x)
3
峰度系数Kurtosis计算