R语言下的线性回归模型

一、一元线性回归

1、定义

只考虑一个因变量和一个自变量之间的关系，其数学表达式为：

其中beta0为模型的截距项，beta1为模型的回归系数，eps为模型的随机误差，一般假定其服从正态分布。

2、回归系数的估计

线性模型中最重要的就是在给定的数据下，如何求得线性模型的回归系数，即上式中的beta0和beta1值。统计学中通过最小二乘法求得线性回归系数的估计值，其基本思想就是在模型的随机误差平方和最小的情况下，求得beta0和beta1值。

随机误差平方和:

最小化随机误差平方和（通过偏微分方程的计算）：

求得回归系数：

其中：

在求得回归系数后，进一步可以计算回归系数的标准差：

其中sigma可以通过下面的式子求得：

根据回归系数的标准差还可以计算其置信区间：

为了更好地理解线性回归模型的最小二乘法思想，这里自定义函数求解一元线性回归模型的参数估计、标准差和置信区间，脚本如下：

estmate <- function(y,x){

mean.x <- mean(x)

mean.y <- mean(y)

sxx <- sum((x-mean.x)^2)

sxy <- sum((y-mean.y)*(x-mean.x))

#计算回归系数

alpha2 <- sxy/sxx

alpha1 <- mean.y-beta*mean.x

#计算回归系数的标准差

n <- length(x)

sigma <- sqrt(sum((y-alpha-beta*x)^2)/(n-2))

sd.alpha1 <- sigma*sqrt(1/n+mean.x^2/sxx)

sd.alpha2 <- sigma/sqrt(sxx)

#计算回归系数的置信区间

p <- 0.05

alpha1l <- alpha1-sd.alpha1*qt(1-p/2,df = n-2)

alpha1u <- alpha1+sd.alpha1*qt(1-p/2,df = n-2)

alpha2l <- alpha2-sd.alpha2*qt(1-p/2,df = n-2)

alpha2u <- alpha2+sd.alpha2*qt(1-p/2,df = n-2)

#返回参数估计、标准差和置信区间

matrix(c(alpha1,alpha2,sd.alpha1,sd.alpha2,alpha1l,alpha2l,alpha1u,alpha2u),nrow = 2, dimnames = list(c('alpha1','alpha2'),c('Estmate','SD','LCI','UCI')))

}

以R中自带的women数据集做一个测试：

head(women)

estmate(women$weight,women$height)

其实，R中可以直接通过lm函数和confint函数求得模型的参数估计、标准差和置信区间。

看一下lm函数语法和参数：

lm(formula, data, subset, weights, na.action,

method = "qr", model = TRUE, x = FALSE,

y = FALSE, qr = TRUE,singular.ok = TRUE,

contrasts = NULL, offset, ...)

formula指定回归模型的公式表达式，类似于y=a+bx；

data指定要分析的数据对象；

subset为可选向量，表示观测值的子集；

weights为可选向量，用于数据拟合的权重；

method指定模型拟合的方法。

该结果与自定义函数返回的结果一致。

3、回归系数和方程的显著性检验

一旦模型的参数估计计算出来之后，需要进一步考察模型的回归系数是否显著、模型是否显著，自变量对因变量的解释程度如何？下面对这三个问题做一个解释：

回归系数是否显著

对回归系数进行假设检验时，通常原假设为“某一个回归系数为0”。可以通过t检验判决该回归系数是否显著，从lm函数的结果中可以查看到每个回归系数对应的P值，一般P值越小，则原假设成立的概率也就越小，上图中显示截距项和height是非常显著的，P值均远远小于0.05。

回归方程是否显著

对回归方程进行假设检验时，通常原假设为“所有回归系数均为为0”。可以通过F检验判断模型是否显著，从lm函数结果中的最后一行可知该模型是非常显著的。

自变量对因变量的解释程度

可以通过R-Square(相关系数平方值)得到模型中的自变量多大程度上解释了因变量，但R-Square存在一个缺陷，就是自变量越多时，R-Square一定越大，为克服该缺点，可以引入修正的R-Square值，从lm结果中得知R方为99.1%。

4、模型的预测

回归模型的最大用处就是预测，只有在模型通过各种检验和假设条件下，方程才是有意义的，其预测结果才会准确。计算预测值及其置信区间的数学表达式为：

R中可以通过predict函数实现模型预测值和置信区间，假设一个女士的身高为66，那么其体重可能为多少呢？置信区间又是多少呢？

newdata <- data.frame(height = 66)

weight.p <- predict(fit,newdata=newdata,interval='prediction')

weight.p

二、多元线性回归

在实际问题中影响因变量的因素不止一个，多元线性回归模型的数学表达式可表示为：

其余关于模型的参数估计、标准差、置信区间等都类似于一元线性回归，只是将单变量变换为了多元变量。具体的数学表达式由下方给出：

参数估计

参数标准差

参数置信区间

应用：

下面采用lm函数拟合一个多元线性回归模型，数据采用R包自带的state.x77，对美国各个洲谋杀率建立多元线性回归模型。

str(as.data.frame(state.x77))

my.data <- as.data.frame(state.x77[,c(1,2,3,4,5,6)])

fit <- lm(Murder~., data = my.data)

summary(fit)

从上面的拟合结果可知各个洲的人口数、文盲率、寿命三个变量是显著的，P值均小于0.05；模型也是显著的；同时这5个自变量对谋杀率的解释率达到76.4%。但是模型中还存在两个变量是不显著的，即收入和高中毕业率，这是否跟模型显著存在悖论呢？模型的好坏是否还有其他指标可以衡量呢？在下一个系列的文章中将重点介绍线性模型的诊断。

参考书目

R语言实战

统计建模与R语言

R语言与网站分析

总结：文中涉及到的R包和函数

stats包

lm()

summary()

confint()

predict()