R lm 模型中的变量顺序

机器算法验证 r 回归 分类数据 固定效应模型 流明
2022-04-12 13:52:41

R模型中变量的顺序是否应该很重要?出于某种原因,以下两个模型导致与 fm 和 yr 相关的不同系数(它们应该分别模拟与 fm 和 yr 相关的固定效应):

set.seed(0)

fm = c(rep("A", 5), rep("B", 5), rep("C", 5), rep("D", 5))
yr = rep(c(1985,1986,1987,1988,1989), 4)
d = data.frame(yr, fm, y=rnorm(length(yr)), x=rnorm(length(yr)))

out_1 = lm(y~x+factor(yr)+factor(fm)-1, data= d)
out_2 = lm(y~x+factor(fm)+factor(yr)-1, data=d)
2个回答

更新:

方程中只有一个截距。截距包括与因子 A 和 1985 年相关的观测值(模型 3 就是这种情况)。但是,在您的第一种情况下,您忽略了因子 A(因此它作为基数),而在模型 2 中,您使用 1985 年作为基数。因此,系数必须不同,因为您要与不同的基数进行比较。如果要运行固定效果,那么应该只有 7 个假人(4 个代表年,3 个代表 fm,如 out_3 中的)。

摘要(out_1)

Call:
lm(formula = y ~ x + factor(yr) + factor(fm) - 1, data = d)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.5538 -0.5543 -0.1142  0.3762  2.1349 

Coefficients:
               Estimate Std. Error t value Pr(>|t|)
x                0.1723     0.4289   0.402    0.696
factor(yr)1985   0.7434     0.7129   1.043    0.319
factor(yr)1986   0.2433     0.7349   0.331    0.747
factor(yr)1987   0.5862     0.7015   0.836    0.421
factor(yr)1988   1.1247     0.6945   1.619    0.134
factor(yr)1989   1.0779     0.6981   1.544    0.151
factor(fm)B     -0.8163     0.7025  -1.162    0.270
factor(fm)C     -1.0703     0.7171  -1.493    0.164
factor(fm)D     -1.2179     0.7067  -1.723    0.113

Residual standard error: 1.095 on 11 degrees of freedom
Multiple R-squared: 0.3347, Adjusted R-squared: -0.2096 
F-statistic: 0.615 on 9 and 11 DF,  p-value: 0.7628 

> summary(out_2)

Call:
lm(formula = y ~ x + factor(fm) + factor(yr) - 1, data = d)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.5538 -0.5543 -0.1142  0.3762  2.1349 

Coefficients:
               Estimate Std. Error t value Pr(>|t|)
x               0.17232    0.42894   0.402    0.696
factor(fm)A     0.74338    0.71286   1.043    0.319
factor(fm)B    -0.07289    0.69443  -0.105    0.918
factor(fm)C    -0.32696    0.69275  -0.472    0.646
factor(fm)D    -0.47447    0.75862  -0.625    0.544
factor(yr)1986 -0.50005    0.77809  -0.643    0.534
factor(yr)1987 -0.15720    0.82354  -0.191    0.852
factor(yr)1988  0.38128    0.78301   0.487    0.636
factor(yr)1989  0.33454    0.77855   0.430    0.676

Residual standard error: 1.095 on 11 degrees of freedom
Multiple R-squared: 0.3347, Adjusted R-squared: -0.2096 
F-statistic: 0.615 on 9 and 11 DF,  p-value: 0.7628 



 out_3 = lm(y~x+factor(fm)+factor(yr), data=d)
    > summary(out_3)

Call:
lm(formula = y ~ x + factor(fm) + factor(yr), data = d)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.2632 -0.2938 -0.1132  0.2488  1.2838 

Coefficients:
               Estimate Std. Error t value Pr(>|t|)   
(Intercept)    -1.04018    0.49474  -2.102  0.05935 . 
x               0.03699    0.22976   0.161  0.87501   
factor(fm)B     0.06037    0.48535   0.124  0.90325   
factor(fm)C    -0.27198    0.49027  -0.555  0.59017   
factor(fm)D     0.34259    0.50111   0.684  0.50833   
factor(yr)1986  1.14448    0.62217   1.839  0.09296 . 
factor(yr)1987  1.39348    0.54603   2.552  0.02690 * 
factor(yr)1988  1.95007    0.54562   3.574  0.00436 **
factor(yr)1989  1.34118    0.57869   2.318  0.04075 * 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 0.7643 on 11 degrees of freedom
Multiple R-squared: 0.5898, Adjusted R-squared: 0.2915 
F-statistic: 1.977 on 8 and 11 DF,  p-value: 0.146 
# install.packages("tidyverse")
library(dplyr)
library(tibble)

set.seed(0)

d <- tibble(
  fm = c(rep("A", 5), rep("B", 5), rep("C", 5), rep("D", 5)),
  yr = rep(c(1985, 1986, 1987, 1988, 1989), 4),
  y = rnorm(length(yr)), 
  x = rnorm(length(yr))
  ) %>%
  mutate(shock = ifelse(yr == 1988, 1, 0))

鉴于最近关于这个主题的帖子,我想我会强调 R 使用变量排序来打破共线性。例如,factor(fm) - 1在第一次运行和factor(yr) - 1第二次运行中可能给人的印象是我们正在明确决定要降低哪个因素水平;不是这种情况。软件根据可变排序做出妥协。这是您的第一个模型:

summary(lm(y ~ -1 + x + factor(fm) + factor(yr), data = d))

为了避免混淆,我将 移到-1了波浪号的右侧。在这里,-1删除截距(即,1 的列)。下一个变量,从左到右,是用 表示的公司固定效应factor(fm)现在将估计所有级别。第三个变量是用 表示的年份固定效应factor(yr)必须删除一个级别(即 1985 年) 。请记住,代表您的公司虚拟对象的所有列总和为单位,因此 R 将下降一年效应以打破共线性。请参阅下面的删节输出:

Call:
lm(formula = y ~ -1 + x + factor(fm) + factor(yr), data = d)

Coefficients:
               Estimate Std. Error t value Pr(>|t|)
x               0.17232    0.42894   0.402    0.696
factor(fm)A     0.74338    0.71286   1.043    0.319
factor(fm)B    -0.07289    0.69443  -0.105    0.918
factor(fm)C    -0.32696    0.69275  -0.472    0.646
factor(fm)D    -0.47447    0.75862  -0.625    0.544
factor(yr)1986 -0.50005    0.77809  -0.643    0.534
factor(yr)1987 -0.15720    0.82354  -0.191    0.852
factor(yr)1988  0.38128    0.78301   0.487    0.636
factor(yr)1989  0.33454    0.77855   0.430    0.676

现在假设我们交换公司效应和年份效应。该模型具有与之前相同的结构:

summary(lm(y ~ -1 + x + factor(yr) + factor(fm), data = d))

我们仍然-1在函数内部引入,lm()但现在固定效应排在第三位。结果,软件将排除公司。和以前一样,您的年份虚拟变量系列总和为单位,所以现在从估计中删除了一个确定的效果。请参阅下面的删节输出:

Call:
lm(formula = y ~ -1 + x + factor(yr) + factor(fm), data = d)

Coefficients:
               Estimate Std. Error t value Pr(>|t|)
x                0.1723     0.4289   0.402    0.696
factor(yr)1985   0.7434     0.7129   1.043    0.319
factor(yr)1986   0.2433     0.7349   0.331    0.747
factor(yr)1987   0.5862     0.7015   0.836    0.421
factor(yr)1988   1.1247     0.6945   1.619    0.134
factor(yr)1989   1.0779     0.6981   1.544    0.151
factor(fm)B     -0.8163     0.7025  -1.162    0.270
factor(fm)C     -1.0703     0.7171  -1.493    0.164
factor(fm)D     -1.2179     0.7067  -1.723    0.113

如果包含拦截,则在这种情况下排序无关紧要。应该清楚的是,其中一个公司假人和一个年份假人将与附加到模型矩阵的 1 列共线。运行第三个模型,但在此试验中不要放弃截距。它应该产生 3 个独立的公司效应和 4 个独立的年份效应。我不确定为什么x上一个答案中的所有估计都不相等(@SextusEmpiricus 很好地指出);它们应该是相似的。

当您的模型在i或者t. 假设您怀疑 1988 年对公司造成了宏观层面的冲击。您通过在面板中包含特定于所有实体的时间虚拟变量来模拟冲击。我将在时间效应之前运行一个模型,然后再运行第二个模型,包括时间效应之后的冲击。这是第一次运行:

Call:
lm(formula = y ~ x + factor(fm) + shock + factor(yr), data = d)

Coefficients: (1 not defined because of singularities)
               Estimate Std. Error t value Pr(>|t|)
(Intercept)      0.7434     0.7129   1.043    0.319
x                0.1723     0.4289   0.402    0.696
factor(fm)B     -0.8163     0.7025  -1.162    0.270
factor(fm)C     -1.0703     0.7171  -1.493    0.164
factor(fm)D     -1.2179     0.7067  -1.723    0.113
shock            0.3813     0.7830   0.487    0.636
factor(yr)1986  -0.5001     0.7781  -0.643    0.534
factor(yr)1987  -0.1572     0.8235  -0.191    0.852
factor(yr)1988       NA         NA      NA       NA
factor(yr)1989   0.3345     0.7785   0.430    0.676

请注意,1988 年时间冲击于变量排序中的年份固定效应。必须去掉一年效应才能打破共线性;R 进入时间效应。现在我将交换排序并重新估计模型:

Call:
lm(formula = y ~ x + factor(fm) + factor(yr) + shock, data = d)

Coefficients: (1 not defined because of singularities)
               Estimate Std. Error t value Pr(>|t|)
(Intercept)      0.7434     0.7129   1.043    0.319
x                0.1723     0.4289   0.402    0.696
factor(fm)B     -0.8163     0.7025  -1.162    0.270
factor(fm)C     -1.0703     0.7171  -1.493    0.164
factor(fm)D     -1.2179     0.7067  -1.723    0.113
factor(yr)1986  -0.5001     0.7781  -0.643    0.534
factor(yr)1987  -0.1572     0.8235  -0.191    0.852
factor(yr)1988   0.3813     0.7830   0.487    0.636
factor(yr)1989   0.3345     0.7785   0.430    0.676
shock                NA         NA      NA       NA

在这里,R在时间效应之后shock按照指定的方式下降。请记住,年份固定效应代表每年影响所有公司的常见冲击,因此对 1988 年的冲击建模的时间虚拟模型与年份固定效应是多余的。

这是订购问题的又一个例子。