PLM：无法添加虚拟变量
PLM: Cannot add dummy variable
我目前正在使用 plm() 估计固定效应模型。下面的 table 是我的数据示例（请注意，我在这里使用了任意数字）。我 运行 使用地区和年份固定效应进行回归，正如预期的那样，由于重复的 id-time 而出现错误。因此，我将 district 与 grade 合并在一起以获得回归的唯一 id。
state
district
year
grade
Y
X
id
AK
1001
2009
3
0.1
0.5
1001.3
AK
1001
2010
3
0.8
0.4
1001.3
AK
1001
2011
3
0.5
0.7
1001.3
AK
1001
2009
4
1.5
1.3
1001.4
AK
1001
2010
4
1.1
0.7
1001.4
AK
1001
2011
4
2.1
0.4
1001.4
...
...
...
..
..
..
...
WY
5606
2011
6
4.2
5.3
5606.6
一切都很顺利，直到我尝试在回归中添加年级虚拟变量。我尝试同时使用 factor() 并在等式中添加了虚拟变量。但两者都没有奏效。我没有在结果中看到虚拟变量。请注意，为了简洁起见，我只展示了第一个带有 factor() 的作品。在第二个回归中，我生成了年级虚拟变量，即 g3 和 g4，并将它们放在回归中而不是 factor(grade)。它应该看起来像 plm(formula = Y ~ X + g3 + g4,....
fe <- plm(formula = Y ~ X + factor(grade),
data = df,
index = c("id", "year"),
model = "within",
effect = "twoways")
summary(fe)
Twoways effects Within Model
Call:
plm(formula = Y ~ X + factor(grade), data = df,
effect = "twoways", model = "within", index = c("id",
"year"))
Unbalanced Panel: n = 64302, T = 1-10, N = 499112
Residuals:
Min. 1st Qu. Median 3rd Qu. Max.
-11.35455 -0.34340 0.00000 0.34364 6.42513
Coefficients:
Estimate Std. Error t-value Pr(>|t|)
Y 0.0126717 0.0036019 3.518 0.0004348 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Total Sum of Squares: 173290
Residual Sum of Squares: 173280
R-Squared: 2.8464e-05
Adj. R-Squared: -0.14788
F-statistic: 12.3766 on 1 and 434800 DF, p-value: 0.00043478
问题：为什么会这样？是因为district和id之间的结合吗？如果是这样，我应该如何修复它以获得这些虚拟变量的系数？ plm() 是我应该使用的合适的软件包吗？任何建议，将不胜感激。谢谢！
P.S。绝对不是多重共线性的问题。 This post说跟这个问题有关。我遵循了这个 post，但我的结果得到了 false。
尽管您指出了测试，但这绝对是一个共线性问题。 grade 中没有未被 id 考虑的独立信息。这是一个简单的例子。在这个模型中，唯一的变量是 id 因子——它本质上是估计 grade 的每个值的 y 的平均值，它是截距加上任何虚拟变量的系数特别是 id.
set.seed(123)
dat <- tibble(
dist = sample(LETTERS[1:10], 1000, replace=TRUE),
grade = sample(letters[17:26], 1000, replace=TRUE),
id = paste(dist, grade, sep="-"),
y = rnorm(1000)
)
mod <- lm(y ~ factor(id), data=dat)
现在，假设我们要使用该模型来获取任何 grade 的均值，这里是 y 的 grade 均值。
dat %>%
group_by(grade) %>%
summarise(m = mean(y))
# # A tibble: 10 × 2
# grade m
# <chr> <dbl>
# 1 q -0.0523
# 2 r -0.193
# 3 s -0.0964
# 4 t 0.0647
# 5 u -0.161
# 6 v -0.0273
# 7 w 0.0390
# 8 x 0.109
# 9 y -0.104
# 10 z 0.146
让我们尝试使用模型估计来获得 grade=z 的 y 的 grade 平均值。首先，让我们找出包含 grade=z:
的每个 id 组中的观察百分比
n <- dat %>%
group_by(id) %>%
tally() %>%
filter(str_detect(id, "z$")) %>%
mutate(pct = n/sum(n))
n
# # A tibble: 10 × 3
# id n pct
# <chr> <int> <dbl>
# 1 A-z 11 0.112
# 2 B-z 10 0.102
# 3 C-z 13 0.133
# 4 D-z 6 0.0612
# 5 E-z 12 0.122
# 6 F-z 8 0.0816
# 7 G-z 9 0.0918
# 8 H-z 10 0.102
# 9 I-z 7 0.0714
# 10 J-z 12 0.122
我们现在可以收集截距和包含 grade=z:
的 id 值
ests <- broom::tidy(mod) %>%
filter(str_detect(term, "ntercept|z")) %>%
mutate(term = gsub("factor\(id\)", "", term)) %>%
select(1,2)
ests
# # A tibble: 11 × 2
# term estimate
# <chr> <dbl>
# 1 (Intercept) -0.391
# 2 A-z 0.264
# 3 B-z 0.572
# 4 C-z 0.520
# 5 D-z 0.863
# 6 E-z 0.774
# 7 F-z 0.755
# 8 G-z 0.591
# 9 H-z 0.0538
# 10 I-z -0.0657
# 11 J-z 0.951
然后我们可以将这些数据与上面的百分比结合起来，并将 intercept 项的百分比替换为 1，因为我们想将截距添加到组系数的加权平均值中：
ests <- ests %>%
left_join(n %>% rename(term = id)) %>%
mutate(pct = ifelse(is.na(pct), 1, pct))
ests
# # A tibble: 11 × 4
# term estimate n pct
# <chr> <dbl> <int> <dbl>
# 1 (Intercept) -0.391 NA 1
# 2 A-z 0.264 11 0.112
# 3 B-z 0.572 10 0.102
# 4 C-z 0.520 13 0.133
# 5 D-z 0.863 6 0.0612
# 6 E-z 0.774 12 0.122
# 7 F-z 0.755 8 0.0816
# 8 G-z 0.591 9 0.0918
# 9 H-z 0.0538 10 0.102
# 10 I-z -0.0657 7 0.0714
# 11 J-z 0.951 12 0.122
最后，我们可以将 estimate 列乘以 pct 列求和：
ests %>%
summarise(m = sum(pct*estimate))
# # A tibble: 1 × 1
# m
# <dbl>
# 1 0.146
请注意，这与我们从上面计算的 grade=z 平均值完全相同。这表明我们可以通过使用 id 系数完美地恢复 y 的 grade 均值，这意味着我们无法估计 grade 一次 id 的独立效应] 由于完全共线性已被考虑在内。
我目前正在使用 plm() 估计固定效应模型。下面的 table 是我的数据示例（请注意，我在这里使用了任意数字）。我 运行 使用地区和年份固定效应进行回归，正如预期的那样，由于重复的 id-time 而出现错误。因此，我将 district 与 grade 合并在一起以获得回归的唯一 id。
| state | district | year | grade | Y | X | id |
|---|---|---|---|---|---|---|
| AK | 1001 | 2009 | 3 | 0.1 | 0.5 | 1001.3 |
| AK | 1001 | 2010 | 3 | 0.8 | 0.4 | 1001.3 |
| AK | 1001 | 2011 | 3 | 0.5 | 0.7 | 1001.3 |
| AK | 1001 | 2009 | 4 | 1.5 | 1.3 | 1001.4 |
| AK | 1001 | 2010 | 4 | 1.1 | 0.7 | 1001.4 |
| AK | 1001 | 2011 | 4 | 2.1 | 0.4 | 1001.4 |
| ... | ... | ... | .. | .. | .. | ... |
| WY | 5606 | 2011 | 6 | 4.2 | 5.3 | 5606.6 |
一切都很顺利，直到我尝试在回归中添加年级虚拟变量。我尝试同时使用 factor() 并在等式中添加了虚拟变量。但两者都没有奏效。我没有在结果中看到虚拟变量。请注意，为了简洁起见，我只展示了第一个带有 factor() 的作品。在第二个回归中，我生成了年级虚拟变量，即 g3 和 g4，并将它们放在回归中而不是 factor(grade)。它应该看起来像 plm(formula = Y ~ X + g3 + g4,....
fe <- plm(formula = Y ~ X + factor(grade),
data = df,
index = c("id", "year"),
model = "within",
effect = "twoways")
summary(fe)
Twoways effects Within Model
Call:
plm(formula = Y ~ X + factor(grade), data = df,
effect = "twoways", model = "within", index = c("id",
"year"))
Unbalanced Panel: n = 64302, T = 1-10, N = 499112
Residuals:
Min. 1st Qu. Median 3rd Qu. Max.
-11.35455 -0.34340 0.00000 0.34364 6.42513
Coefficients:
Estimate Std. Error t-value Pr(>|t|)
Y 0.0126717 0.0036019 3.518 0.0004348 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Total Sum of Squares: 173290
Residual Sum of Squares: 173280
R-Squared: 2.8464e-05
Adj. R-Squared: -0.14788
F-statistic: 12.3766 on 1 and 434800 DF, p-value: 0.00043478
问题：为什么会这样？是因为district和id之间的结合吗？如果是这样，我应该如何修复它以获得这些虚拟变量的系数？ plm() 是我应该使用的合适的软件包吗？任何建议，将不胜感激。谢谢！
P.S。绝对不是多重共线性的问题。 This post说跟这个问题有关。我遵循了这个 post，但我的结果得到了 false。
尽管您指出了测试，但这绝对是一个共线性问题。 grade 中没有未被 id 考虑的独立信息。这是一个简单的例子。在这个模型中，唯一的变量是 id 因子——它本质上是估计 grade 的每个值的 y 的平均值，它是截距加上任何虚拟变量的系数特别是 id.
set.seed(123)
dat <- tibble(
dist = sample(LETTERS[1:10], 1000, replace=TRUE),
grade = sample(letters[17:26], 1000, replace=TRUE),
id = paste(dist, grade, sep="-"),
y = rnorm(1000)
)
mod <- lm(y ~ factor(id), data=dat)
现在，假设我们要使用该模型来获取任何 grade 的均值，这里是 y 的 grade 均值。
dat %>%
group_by(grade) %>%
summarise(m = mean(y))
# # A tibble: 10 × 2
# grade m
# <chr> <dbl>
# 1 q -0.0523
# 2 r -0.193
# 3 s -0.0964
# 4 t 0.0647
# 5 u -0.161
# 6 v -0.0273
# 7 w 0.0390
# 8 x 0.109
# 9 y -0.104
# 10 z 0.146
让我们尝试使用模型估计来获得 grade=z 的 y 的 grade 平均值。首先，让我们找出包含 grade=z:
id 组中的观察百分比
n <- dat %>%
group_by(id) %>%
tally() %>%
filter(str_detect(id, "z$")) %>%
mutate(pct = n/sum(n))
n
# # A tibble: 10 × 3
# id n pct
# <chr> <int> <dbl>
# 1 A-z 11 0.112
# 2 B-z 10 0.102
# 3 C-z 13 0.133
# 4 D-z 6 0.0612
# 5 E-z 12 0.122
# 6 F-z 8 0.0816
# 7 G-z 9 0.0918
# 8 H-z 10 0.102
# 9 I-z 7 0.0714
# 10 J-z 12 0.122
我们现在可以收集截距和包含 grade=z:
id 值
ests <- broom::tidy(mod) %>%
filter(str_detect(term, "ntercept|z")) %>%
mutate(term = gsub("factor\(id\)", "", term)) %>%
select(1,2)
ests
# # A tibble: 11 × 2
# term estimate
# <chr> <dbl>
# 1 (Intercept) -0.391
# 2 A-z 0.264
# 3 B-z 0.572
# 4 C-z 0.520
# 5 D-z 0.863
# 6 E-z 0.774
# 7 F-z 0.755
# 8 G-z 0.591
# 9 H-z 0.0538
# 10 I-z -0.0657
# 11 J-z 0.951
然后我们可以将这些数据与上面的百分比结合起来，并将 intercept 项的百分比替换为 1，因为我们想将截距添加到组系数的加权平均值中：
ests <- ests %>%
left_join(n %>% rename(term = id)) %>%
mutate(pct = ifelse(is.na(pct), 1, pct))
ests
# # A tibble: 11 × 4
# term estimate n pct
# <chr> <dbl> <int> <dbl>
# 1 (Intercept) -0.391 NA 1
# 2 A-z 0.264 11 0.112
# 3 B-z 0.572 10 0.102
# 4 C-z 0.520 13 0.133
# 5 D-z 0.863 6 0.0612
# 6 E-z 0.774 12 0.122
# 7 F-z 0.755 8 0.0816
# 8 G-z 0.591 9 0.0918
# 9 H-z 0.0538 10 0.102
# 10 I-z -0.0657 7 0.0714
# 11 J-z 0.951 12 0.122
最后，我们可以将 estimate 列乘以 pct 列求和：
ests %>%
summarise(m = sum(pct*estimate))
# # A tibble: 1 × 1
# m
# <dbl>
# 1 0.146
请注意，这与我们从上面计算的 grade=z 平均值完全相同。这表明我们可以通过使用 id 系数完美地恢复 y 的 grade 均值，这意味着我们无法估计 grade 一次 id 的独立效应] 由于完全共线性已被考虑在内。