# 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 年的冲击建模的时间虚拟模型与年份固定效应是多余的。
这是订购问题的又一个例子。