给女朋友写的生统资料_Part16

多元线性回归

多元线性回归的方程写为:
Y = \beta_0 + \beta_1 X_1 + \beta_2 X_2 + \cdot\cdot\cdot + \beta_p X_p + \epsilon
其中X_j代表第j个预测变量,\beta_j是对应的模型参数。\beta_j可以解释为在所有其他预测变量保持不变的情况下,X_j增加一个单位对Y产生的平均影响。(这段话其实很重要,我觉得有可能要考的)

和前面的一元线性回归类似,我们需要用\hat{\beta}_0,\hat{\beta}_1,……,\hat{\beta}_p来估计\beta_0,\beta_1,……,\beta_p。对于给定的\hat{\beta}_0,\hat{\beta}_1,……,\hat{\beta}_p,我们的回归模型就变成了
\hat{y} = \hat{\beta}_0 + \hat{\beta}_1 x_1 + \hat{\beta}_2 x_2 + \cdot\cdot\cdot + \hat{\beta}_p x_p
与一元线性回归类似,多元线性回归的参数也是用最小二乘法进行估计的,也是使得残差平方和(RSS)最小。计算方法的话,似乎是矩阵代数最为方便。一来我也不会╮(╯_╰)╭,二来老师PPT也没有放,所以这里就不算了。

然后我们还是按照一元线性回归的步骤,首先去检验整个方程的显著性。

模型的显著性检验(模型的意义)

假设检验第一步首先是建立零假设和备则假设(这里跟一元就不一样了)。
H_{0}: \beta_1 = \beta_2 =... = \beta_p = 0

H_{1}: 至少一个\beta_j不等于0

统计检验的F统计量为
F = \frac{(TSS-RSS)/p}{RSS/(n-p-1)} \sim F(p,n-p-1)
然后计算p-value就可以了。

模型参数的显著性

还记得我们在简单线性回归中提到的,模型参数的准确性么。我们利用了 t 统计量,即可以用来做模型参数的置信区间,也可以用来做模型参数的统计检验。但在一元线性回归中,t 检验跟 F 检验几乎可以说是等价的,因为只有一个参数,你检验了那一个参数的显著性,就是等于检验了整个方程的显著性。但在这里,t 检验就非常有用了。其可以检验每个模型参数的统计显著性,而F检验做的,是整个模型的统计检验性

假设检验:
H0:\beta_j=0,其他模型参数不为0

H1:\beta_j \ne0,其他模型参数不为0

然后统计量为:
t = \frac{\hat{\beta_1} - 0}{SE(\hat{\beta_1})} \sim t(n-p-1)
然后计算 p-value 就行了。

模型的比较

我们先前讲过了 F 统计量是检验整体的统计量,t 统计量可以检验单参数的统计量,其实我们还可以检验部分参数的统计量(这个考试的时候,应该是用anova函数来比较含有不同参数的模型)。

在前面 F 统计量的假设检验是假设所有系数都是零。但这里我们想要检验的是规模为q的特定子集为0。
H0:\beta_{p-q+1}=\beta_{p-q+2}……=\beta_{p}=0,其他剩余的p-q个不等于0

H1:规模为q的特定子集不都为0,其他剩余的p-q个不等于0

为方便起见,我们用除了q个变量之外的所有变量建立第二个模型。假设该模型的残差平方和为RSS_0,相应的F统计量为
F=\frac{(RSS_0-RSS)/q}{RSS/(n-p-1)} \sim F(q,n-p-1)
这部分我觉得不太需要掌握。唯一需要掌握的是利用anova函数比较去掉变量前后的模型,这个后面会给出例子。

评估模型的准确性

这部分还是用残差标准误和R^2​统计量暂时应该就可以了。残差标准误稍微注意下
RSE = \sqrt{\frac{RSS}{n-p-1}}
参考文章

  • 《统计学习导论》

R语言实现

R 语言实现的函数就是一个 lm 函数,用法见《R语言实战》第二版162页。我这里放个截图


15_1.png

15_2.png

我这里用的是我们生统第六次作业的第二题,我将这个数据集命名为test3(数据集已放入part0部分)。

某农场通过试验取得早稻收获量与春季降雨量和春季温度的数据如下:

# 准备数据
> test3 <- read.table("rawdata/homework-6.2.txt",header = T)
> colnames(test3) <- c("weight","rainfall","tempature")
> test3
  weight rainfall tempature
1   2250       25         6
2   3450       33         8
3   4500       45        10
4   6750      105        13
5   7200      110        14
6   7500      115        16
7   8250      120        17
  • 建立早稻收获量对春季降雨量和春季温度的二元线性回归方程。

    用的是lm函数,用summary函数展现结果

    > test3_lm <- lm(weight ~ ., data = test3)
    > summary(test3_lm)
    
    Call:
    lm(formula = weight ~ ., data = test3)
    
    Residuals:
           1        2        3        4        5        6        7 
    -275.101   90.464  216.483  140.280  150.676 -316.599   -6.203 
    
    Coefficients:
                Estimate Std. Error t value Pr(>|t|)  
    (Intercept)   -0.591    505.004  -0.001   0.9991  
    rainfall      22.387      9.601   2.332   0.0801 .
    tempature    327.672     98.798   3.317   0.0295 *
    ---
    Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
    
    Residual standard error: 261.4 on 4 degrees of freedom
    Multiple R-squared:  0.9913,  Adjusted R-squared:  0.987 
    F-statistic: 228.4 on 2 and 4 DF,  p-value: 7.532e-05
    

    写出方程
    y= -0.591+22.387x1+327.672x2

  • 检验模型的显著性(模型的意义)**

    第一步 提出假设,设置显著性水平:
    H_{0}: \beta_1 = \beta_2 = 0

    H1:不都为0

    显著性水平\alpha设置为0.05

    第二步 计算检验统计量 F

    F统计量为
    F = \frac{(TSS-RSS)/p}{RSS/(n-p-1)} \sim F(p,n-p-1)
    第三步 做出决策。按 R 输出结果为 228.4,显著性水平 P=7.532e-05 < 0.05。所以拒绝 H_0,即收获量y与降雨量和温度之间的线性关系显著。

  • 判断每个自变量对因变量的影响是否都是显著的

    需要对各回归系数\beta_j分别进行 t 检验,假设检验和t统计量见前面的模型参数的显著性。根据前面线性回归输出的结果,可以看到降雨量和温度的回归系数对应的显著性水平分别为0.0801和0.0295。只有温度对应的显著性水平小于0.05通过检验,这表明影响收获量的自变量中,只有温度对收获量的影响显著,而降雨量对收获量的影响不显著。

  • 各回归系数的置信区间(根据前面截图的函数)

    > confint(test3_lm)
                       2.5 %     97.5 %
    (Intercept) -1402.707516 1401.52552
    rainfall       -4.268921   49.04184
    tempature      53.364699  601.97873
    
  • 检验两个模型的区别

    因为先前我们发现降雨量对收获量的影响不显著,所以我们剔除到降雨量,再重新构建一个模型。然后拿新模型和旧模型比较,看去掉的那个变量是否重要

    # 重新构建模型
    > test3_lm_new <- lm(weight ~ tempature, data = test3)
    > summary(test3_lm_new)
    
    Call:
    lm(formula = weight ~ tempature, data = test3)
    
    Residuals:
       1    2    3    4    5    6    7 
    -150  -50 -100  500  400 -400 -200 
    
    Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
    (Intercept)  -900.00     447.82   -2.01    0.101    
    tempature     550.00      35.56   15.47 2.05e-05 ***
    ---
    Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
    
    Residual standard error: 359.2 on 5 degrees of freedom
    Multiple R-squared:  0.9795,  Adjusted R-squared:  0.9754 
    F-statistic: 239.2 on 1 and 5 DF,  p-value: 2.052e-05
    
    # 利用anova函数比较
    > anova(test3_lm,test3_lm_new)
    Analysis of Variance Table
    
    Model 1: weight ~ rainfall + tempature
    Model 2: weight ~ tempature
      Res.Df    RSS Df Sum of Sq      F  Pr(>F)  
    1      4 273385                              
    2      5 645000 -1   -371615 5.4372 0.08009 .
    ---
    Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
    

    稍微要注意的一点就是,在后面逻辑斯蒂回归的(广义线性模型)比较中,anova要设置test="Chisq"。

  • 其余一些结果

    • 残差标准误(似乎也有叫剩余标准差的):261.4

    • R^2统计量:0.9913。

      这里你会发现还有一个 adjust\ R^2。因为实际上只要你的变量加的足够多,你的R^2就肯定会增加,哪怕你的变量跟你的Y值是无关的。所以就会引入adjust,来适当调整R^2。我觉得考试时候两个值可以都写。

置信区间与预测区间

这部分我觉得大家不用看,我只是写一下,考试肯定也不会考。

先提下可约和不可约误差。

我们利用统计模型,对Y做了预测,得到预测值\hat Y。预测值的精确度取决于两个方面,一方面是可约误差,另一方面是不可约误差。其中可约误差是来自于我们所构建的模型。如果我们构建的模型不是那么好,可约误差就会大一点。但我们可以通过优化我们的模型,比如增加一些变量或者把线性变成非线性等等,来降低这种误差。然而,尽管我们可以得到一个精确的模型,我们还是不能减少不可约误差。因为不可约误差是跟误差项,即我们在最一开始的时候写过
Y = \beta_0 + \beta_1 X_1 + \beta_2 X_2 + \cdot\cdot\cdot + \beta_p X_p + \epsilon
里面的\epsilon有关的。按照定义,\epsilon 是不能用X去预测的。所以我们是无法去减少这部分误差的。那么其实最后的误差来源就表示为
E(Y-\hat{Y})^2=E[f(X)+\epsilon-\hat f(X)]^2

=[f(X)-\hat f(X)]^2+Var(\epsilon)

我为啥要提上面这个呢。因为其实在针对线性回归有一个函数叫predict,里面有个参数叫interval。大家可以看下。

## S3 method for class 'lm'
predict(object, newdata, se.fit = FALSE, scale = NULL, df = Inf,
        interval = c("none", "confidence", "prediction"),
        level = 0.95, type = c("response", "terms"),
        terms = NULL, na.action = na.pass,
        pred.var = res.var/weights, weights = 1, ...)

书里面一般写的是prediction。prediction其实就是在构建区间的时候考虑了误差项的影响,而confidence则没有考虑误差项的影响。这就会导致预测区间比置信区间要宽很多。

顺便提一下,y_0=E(y_0)+\epsilonE(y_0)就是我们通常用线性回归预测出来的值。而confidence构建的就是E(y_0)的区间,而prediction构建的就是y_0的区间。

# 利用先前的结果,看下区别。可以看到预测区间更宽一点。
> predict(test3_lm,test3,interval = "confidence")
       fit      lwr      upr
1 2525.101 1992.745 3057.457
2 3359.536 2928.396 3790.676
3 4283.517 3795.736 4771.298
4 6609.720 6096.023 7123.417
5 7049.324 6620.305 7478.342
6 7816.599 7407.001 8226.198
7 8256.203 7748.618 8763.789
> predict(test3_lm,test3,interval = "prediction")
       fit      lwr      upr
1 2525.101 1624.957 3425.245
2 3359.536 2515.298 4203.774
3 4283.517 3408.996 5158.038
4 6609.720 5720.483 7498.956
5 7049.324 6206.167 7892.481
6 7816.599 6983.156 8650.043
7 8256.203 7370.483 9141.923

参考文献

©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 205,132评论 6 478
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 87,802评论 2 381
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 151,566评论 0 338
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 54,858评论 1 277
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 63,867评论 5 368
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 48,695评论 1 282
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 38,064评论 3 399
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 36,705评论 0 258
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 42,915评论 1 300
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 35,677评论 2 323
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 37,796评论 1 333
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 33,432评论 4 322
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 39,041评论 3 307
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 29,992评论 0 19
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 31,223评论 1 260
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 45,185评论 2 352
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 42,535评论 2 343

推荐阅读更多精彩内容