电子发烧友App
R语言入门教程四：常用的统计推断
通常一个研究项目能够获得的数据是有限的，以有限的样本特征来推断总体特征就称为统计推断。推断又可细分为区间估计和假设检验，二者虽有区别，但却是一枚硬币的两面，之间有着紧密的关联。
1 对总体均值进行区间估计
假设我们从总体中抽得一个样本，希望根据样本均值判断总体均值的置信区间，如下例所示：
x=rnorm（50，mean=10，sd=5） #随机生成50个均值为10，标准差为5的随机数为作为研究对象
mean（x）-qt（0.975，49）*sd（x）/sqrt（50） #根据统计学区间估计公式，得到95%置信度下的区间下界
mean（x）+qt（0.975，49）*sd（x）/sqrt（50） #95%置信度下的区间上界
也可以直接利用R语言内置函数t.test
t.test（x，conf.level=0.95）
从如下结果可得95%置信区间为（9.56，12.36）
One Sample t-test
data： x
t = 15.7301， df = 49， p-value 《 2.2e-16
alternative hypothesis： true mean is not equal to 0
95 percent confidence interval：
9.563346 12.364729
sample estimates：
mean of x
10.96404
2 对总体均值进行假设检验
还是以上面的X数据作为对象，来检验总体均值是否为10
t.test（x，mu=10，alternative=‘two.sided’） #这里的原假设是总体均值（mu）为10，使用双侧检验，得到P值为0.17，可见P值不够小，不能够拒绝原假设。
T检验是极为常用的检验方法，除了单样本推断之外，t.test命令还可以实现两样本推断和配对样本推断。如果要对总体比率或总体方差进行推断，可以使用prop.test和var.test。
3 正态分布检验
T检验的前提条件是总体服从正态分布，因此我们有必要先检验正态性。而且在评价回归模型时，对残差也需要检验正态性。检验正态性的函数是shapiro.test
shapiro.test（x）
结果所下：
Shapiro-Wilk normality test
data： x
W = 0.9863， p-value = 0.8265
该检验的原假设是服从正态分布，由P值为0.82可判断不能拒绝总体服从正态的假设
4 非参数检验
如果总体不服从正态分布，那么T检验就不再适用，此时我们可以利用非参数方法推断中位数。wilcoxon.test函数可实现符号秩检验。
wilcox.test（x，conf.int=T） #指定conf.int让函数返回中位数的置信区间
wilcox.test（x，mu=1） #指定mu让函数返回中位数为10的检验结果
5 独立性检验（联列表检验）
卡方分布有一个重要应用就是根据样本数据来检验两个分类变量的独立性，我们以CO2数据为例来说明chisq.test函数的使用，help（CO2）可以了解更多信息。
data（CO2） #读入内置的数据包，其中Type和Treatmen是其中两个分类变量。
chisq.test（table（CO2$Type，CO2$Treatment）） #使用卡方检验函数来检验这两个因子之间是否独立
结果显示P值为0.82，因此可以认为两因子之间独立。在样本较小的情况下，还可以使用fisher精确检验，对应的函数是fisher.test。
R语言入门教程五：简单线性回归
线性回归可能是数据分析中最为常用的工具了，如果你认为手上的数据存在着线性定量关系，不妨先画个散点图观察一下，然后用线性回归加以分析。下面简单介绍一下如何在R中进行线性回归。
1 回归建模
我们利用R语言中内置的trees数据，其中包含了Volume（体积）、Girth（树围）、Height（树高）这三个变量，我们希望以体积为因变量，树围为自变量进行线性回归。
plot（Volume~Girth，data=trees，pch=16，col=‘red’）
model=lm（Volume~Girth，data=trees）
abline（model，lty=2）
summary（model）
首先绘制了两变量的散点图，然后用lm函数建立线性回归模型，并将回归直线加在原图上，最后用summary将模型结果进行了展示，从变量P值和F统计量可得回归模型是显著的。但截距项不应该为负数，所以也可以用下面方法将截距强制为0。
model2=lm（Volume~Girth-1，data=trees）
2 模型诊断
在模型建立后会利用各种方式来检验模型的正确性，对残差进行分析是常见的方法，下面我们来生成四种用于模型诊断的图形。
par（mfrow=c（2，2））
plot（model）
par（mfrow=c（1，1））
这里左上图是残差对拟合值作图，整体呈现出一种先下降后下升的模式，显示残差中可能还存在未提炼出来的影响因素。右上图残差QQ图，用以观察残差是否符合正态分布。左下图是标准化残差对拟合值，用于判断模型残差是否等方差。右下图是标准化残差对杠杆值，虚线表示的cooks距离等高线。我们发现31号样本有较大的影响。
3 变量变换
因为31号样本有着高影响力，为了降低其影响，一种方法就是将变量进行开方变换来改善回归结果，从残差标准误到残差图，各项观察都说明变换是有效的。
plot（sqrt（Volume）~Girth，data=trees，pch=16，col=‘red’）
model2=lm（sqrt（Volume）~Girth，data=trees）
abline（model2，lty=2）
summary（model2）
4 模型预测
下面根据上述模型计算预测值以及置信区间，predict函数可以获得模型的预测值，加入参数可以得到预测区间
plot（sqrt（Volume）~Girth，data=trees，pch=16，col=‘red’）
model2=lm（sqrt（Volume）~Girth，data=trees）
data.pre=data.frame（predict（model2，interval=‘prediction’））
lines（data.pre$lwr~trees$Girth，col=‘blue’，lty=2）
lines（data.pre$upr~trees$Girth，col=‘blue’，lty=2）
我们还可以将树围和树高都加入到模型中去，进行多元回归。如果要考虑的变量很多，可以用step函数进行变量筛选，它是以AIC作为评价指标来判断一个变量是否应该加入模型，建议使用这种自动判断函数时要谨慎。对于嵌套模型，还可以使用anova建立方差分析表来比较模型。对于变量变换的形式，则可以使用MASS扩展包中的boxcox函数来进行COX变换。
R语言入门教程六（完）：Logistic回归
让我们用logistic回归来结束本系列的内容吧，本文用例来自于John Maindonald所著的《Data Analysis and Graphics Using R》一书，其中所用的数据集是anesthetic，数据集来自于一组医学数据，其中变量conc表示麻醉剂的用量，move则表示手术病人是否有所移动，而我们用nomove做为因变量，因为研究的重点在于conc的增加是否会使nomove的概率增加。
首先载入数据集并读取部分文件，为了观察两个变量之间关系，我们可以利cdplot函数来绘制条件密度图。
library（DAAG）
head（anesthetic）
cdplot（factor（nomove）~conc，data=anesthetic，main=‘条件密度图’，ylab=‘病人移动’，xlab=‘麻醉剂量’）
从图中可见，随着麻醉剂量加大，手术病人倾向于静止。下面利用logistic回归进行建模，得到intercept和conc的系数为-6.47和5.57，由此可见麻醉剂量超过1.16（6.47/5.57）时，病人静止概率超过50%。
anes1=glm（nomove~conc，family=binomial（link=‘logit’），data=anesthetic）
summary（anes1）
上面的方法是使用原始的0-1数据进行建模，即每一行数据均表示一个个体，另一种是使用汇总数据进行建模，先将原始数据按下面步骤进行汇总
anestot=aggregate（anesthetic［，c（‘move’，‘nomove’）］，by=list（conc=anesthetic$conc），FUN=sum）
anestot$conc=as.numeric（as.character（anestot$conc））
anestot$total=apply（anestot［，c（‘move’，‘nomove’）］，1，sum）
anestot$prop=anestot$nomove/anestot$total
得到汇总数据anestot如下所示
conc move nomove total prop
1 0.8 6 1 7 0.1428571
2 1.0 4 1 5 0.2000000
3 1.2 2 4 6 0.6666667
4 1.4 2 4 6 0.6666667
5 1.6 0 4 4 1.0000000
6 2.5 0 2 2 1.0000000
对于汇总数据，有两种方法可以得到同样的结果，一种是将两种结果的向量合并做为因变量，如anes2模型。另一种是将比率做为因变量，总量做为权重进行建模，如anes3模型。这两种建模结果是一样的。
anes2=glm（cbind（nomove，move）~conc，family=binomial（link=‘logit’），data=anestot）
anes3=glm（prop~conc，family=binomial（link=‘logit’），weights=total，data=anestot）
根据logistic模型，我们可以使用predict函数来预测结果，下面根据上述模型来绘图
x=seq（from=0，to=3，length.out=30）
y=predict（anes1，data.frame（conc=x），type=‘response’）
plot（prop~conc，pch=16，col=‘red’，data=anestot，xlim=c（0.5，3），main=‘Logistic回归曲线图’，ylab=‘病人静止概率’，xlab=‘麻醉剂量’）
lines（y~x，lty=2，col=‘blue’）
工商网监
评论