R语言学习之数据分析实战(一)

一、线性回归

回归（regression）：通常指那些用一个或多个预测变量，也称自变量或解释变量，来预测响应变量，也称为因变量、效标变量或结果变量的方法。

普通最小二乘回归法（OLS）

以women数据集为例：

lm()函数

formula：需要拟合的形式，y~x
data：需要使用的数据集，数据框的形式

R表达式中常用的符号

使用lm()函数对women数据集进行回归分析

> fit <- lm(weight ~ height, data=women)
> summary(fit)

Call这一列列出了使用的回归的公式

Residuals为残差的分布，即真实值和预测值之间的差

Coefficients为系数项，其中Intercept为截距项（即当x=0时与y轴的交点），Estimate为系数，因此关系方程式应为

weight = 3.45*Height - 87.52

Pr为估计系数为零假设的概率，小于0.05较好

估计模型时首先查看F统计量，若F的p-value大于0.05，则该模型无价值，需重新进行拟合。如果小于0.05，再观察r方，模型可以解释多少变量。

线性拟合常用的函数

拟合weight和height的平方的线性函数，不能直接用^,需要将height添加到I()函数中

fit2 <- lm(weight ~ height + I(height^2), data=women)
summary(fit2)

上述结果与之前结果比较，发现残差值明显小了许多，p值也更小，说明这个结果更好，R方值为99.95%，说明这个模型可以解释99.95%的情况。

下面将两次拟合的效果反应在图中

plot(women$height,women$weight,main="Women Age 30-39",xlab="Height (in inches)",ylab="Weight (in lbs)")
lines(women$height,fitted(fit2))
abline(fit)
lines(women$height,fitted(fit2),col="red")

二、多元线性回归

本次利用state.x77数据集来研究犯罪率和其他因素的关系

首先将数据集转为数据框（因为lm函数需要输入数据框格式）

states <- as.data.frame(state.x77[,c("Murder", "Population","Illiteracy", "Income", "Frost")])

利用lm()进行回归分析

fit <- lm(Murder ~ Population + Illiteracy + Income + Frost, data=states)
summary(fit)

可以使用coef()函数查看变量的系数

coef(fit)

在很多研究中，变量之间有交互项，也就是变量之间并不是独立的，比如mtcars数据集，汽车的重量和马力之间存在交互，下面拟合一下mtcars数据集中每加仑行驶里程数mpg变量与马力hp变量和车重vt变量之间的关系。

#由于hp和wt存在交互关系但并不清楚是何种关系，可以用hp:wt表示
fit <- lm(mpg ~ hp + wt + hp:wt, data=mtcars)
summary(fit)

在hp:wt的交互项中为3颗星，说明马力和车重的交互项是非常显著的，这就意味着响应变量mpg和其中一个预测变量的关系依赖于另一个预测变量的水平。

如果存在多个变量，除去响应变量之外的自变量有多种组合形式，那么如何从众多可能的模型中选择最佳的模型呢？

对于上述问题，可以使用AIC()比较模型

fit1 <- lm (Murder ~ Population+Illiteracy+Income+Frost,data=states)
fit2 <- lm (Murder ~ Population+Illiteracy,data=states)
AIC(fit1,fit2)

上述结果显示，fit2的拟合效果更好。

当模型较多时，用AIC()两两比较就显然不行，这时对于变量的选择可以用逐步回归法和全子集回归法。逐步回归法中，模型依次增加或删除变量，直到达到某个节点为止，即继续添加或删除变量模型不再变化，若是每次增加变量，则是向前回归，若是每次删除变量，则是向后回归。全子集回归法则是取出所有可能的模型，从中计算出最佳模型。这两种方法都需要通过R扩展包来实现，逐步回归法可以通过MASS包中的stepAIC()函数进行计算，全子集回归法可以通过leaps包中的regsubsets()函数进行计算。

# Backward stepwise selection
library(MASS)
states <- as.data.frame(state.x77[,c("Murder", "Population","Illiteracy", "Income", "Frost")])
fit <- lm(Murder ~ Population + Illiteracy + Income + Frost,data=states)
stepAIC(fit, direction="backward")#  All subsets regression
library(leaps)
states <- as.data.frame(state.x77[,c("Murder", "Population","Illiteracy", "Income", "Frost")])
leaps <-regsubsets(Murder ~ Population + Illiteracy + Income +Frost, data=states, nbest=4)
plot(leaps, scale="adjr2")

三、回归诊断

满足OLS模型统计假设

利用women数据集通过lm()函数进行分析

fit2 <- lm(weight ~ height + I(height^2), data=women)
opar <- par(no.readonly=TRUE)
par(mfrow=c(2,2))
plot(fit2)

为理解这些图形，我们来回顾一下OLS回归的统计假设。
(1)正态性（主要使用QQ图）当预测变量值固定时，因变量成正态分布，则残差值也应该是一个均值为0的正态分布。正态Q-Q图（Normal Q-Q，右上）是在正态分布对应的值下，标准化残差的概率图。若满足正态假设，那么图上的点应该落在呈45度角的直线上；若不是如此，那么就违反了正态性的假设。
(2)独立性你无法从这些图中分辨出因变量值是否相互独立，只能从收集的数据中来验证。上面的例子中，没有任何先验的理由去相信一位女性的体重会影响另外一位女性的体重。假若你发现数据是从一个家庭抽样得来的，那么可能必须要调整模型独立性的假设。
(3)线性（使用左上角的图，该曲线尽量拟合所有点）若因变量与自变量线性相关，那么残差值与预测（拟合）值就没有任何系统关联。换句话说，除了白噪声，模型应该包含数据中所有的系统方差。在“残差图与拟合图”Residuals vs Fitted，左上）中可以清楚的看到一个曲线关系，这暗示着你可能需要对回归模型加上一个二次项。
(4)同方差性（左下角，点随机分布在曲线的周围）若满足不变方差假设，那么在位置尺度图（Scale-Location Graph，左下）中，水平线周围的点应该随机分布。该图似乎满足此假设。最后一幅“残差与杠图”（Residuals vs Leverage，右下）提供了你可能关注的单个观测点的信息。从图形可以鉴别出离群点、高杠杆值点和强影响点

抽样法验证

1、数据集中有1000个样本，随机抽取500个数据进行回归分析；

2、模型建好后，利用predict函数，对剩余500个样本进行预测，比较残差值；

3、如果预测准确，说明模型可以，否则就需要调整模型；

四、方差分析

方差分析，称为Analysis of Variance,简称ANOVA,也称为“变异数分析”，用于两个及两个以上样本均数差别的显著性检验。从广义上来讲，方差分析也属于回归分析的一种。只不过线性回归的因变量一般是连续型变量。而当自变量是因子时，研究关注的重点通常会从预测转向不同组之间差异的比较。这就是方差分析。

R中因子的应用

计算频数

独立性检验

相关性检验

方差分析

主成分分析

因子分析

……

方差分析分为多种：

单因素方差分析ANOVA(组内，组间)
双因素方差分析ANOVA
协方差分析ANCOVA
多元方差分析MANOVA
多元方差分析MANCOVA

方差公式中特殊符号

各种方差分析公式写法

表达式中变量的顺序很重要

单因素方差分析

这里使用multcomp包中的cholesterol数据集进行演示，该数据集是50个患者接受降低胆固醇药物治疗的五种疗法的数据。

首先利用aggregate()函数分组统计

aggregate(weight, by=list(dose), FUN=mean)

然后使用aov函数进行方差分析

fit <- aov(weight ~ gesttime + dose)
summary(fit)

上述结果p-value为9.8e-13,说明五种治疗效果不同。

五、功效分析

功效分析，power analysis，可以帮助在给定置信度的情况下，判断检测到给定效应值时所需的样本量。反过来，它也可以在给定置信度水平的情况下，计算在某样本量内能检测到给定效应值的概率。

功效分析函数

功效分析的理论基础

1、样本大小指的是实验设计中每种条件/组中观测的数目。

2、显著性水平（也称为alpha）由I型错误的概率来定义。也可以把它看做是发现效应不发生的概率。

3、功效通过I减去II型错误的概率来定义。我们可以把它看做是真实效应发生的概率。

4、效应值指的是在备择或研究假设下效应的量。效应值的表达式依赖于假设检验中使用的统计方法。

R中使用pwr包进行功效分析

线性回归的功效分析可以使用pwr包中的pwr.f2.text函数进行分析

pwr.f2.test(u = NULL, v = NULL, f2 = NULL, sig.level = 0.05, power = NULL)
Arguments

其中sig.level 表示显著性水平，默认为0.05

power为功效水平

u和v分别是分子自由度和分母自由度

f2是效应值,是根据一些公式计算出来的

pwr.f2.test(u=3, f2=0.0769, sig.level=0.05, power=0.90)

得到结果v = 184.2426，也就是假定显著性水平为0.05，在90%置信度的情况下，至少需要185个受治者才可以

平衡单因素方差的功效分析可以使用pwr.anova.test()函数进行分析

pwr.anova.test(k = NULL, n = NULL, f = NULL, sig.level = 0.05, power = NULL)

其中k是组的个数，n是各组的样本大小，也就是我们需要求的量

sig.level 表示显著性水平，默认为0.05power为功效水平

pwr.anova.test(k=2,f=.25,sig.level=.05,power=.9)

功效分析可以使用pwr.anova.test()函数进行分析**

pwr.anova.test(k = NULL, n = NULL, f = NULL, sig.level = 0.05, power = NULL)

其中k是组的个数，n是各组的样本大小，也就是我们需要求的量

sig.level 表示显著性水平，默认为0.05power为功效水平

pwr.anova.test(k=2,f=.25,sig.level=.05,power=.9)

求得n=85.03128，也就是说每组至少需要86个样本，总体样本也就是它的2（k）倍，即172个。