分类多元回归 R 的截距不是平均值吗？

Question

假设我有 2 个（分类）变量和一个连续变量：

library(tidyverse)
set.seed(123)
ds <- data.frame(
  depression=rnorm(90,10,2),
  schooling_dummy=c(0,1,2),
  sex_dummy=c(0,1)
)

当我回归对性的抑郁（0或1）时，截距为10.0436，sex = 0的平均值是多少。好的！

ds %>% group_by(sex_dummy) %>% 
+   summarise(formatC(mean(depression),format="f", digits=4))
# A tibble: 2 x 2
  sex_dummy `formatC(mean(depression), format = "f", digits = 4)`
      <dbl> <chr>                                                
1      0    10.0436                                              
2      1.00 10.1640

当我回归对学校教育的抑郁时，也会发生同样的事情。截距值为10.4398。学校教育的平均值 = 0 是一样的。

ds %>% group_by(schooling_dummy) %>% 
+   summarise(formatC(mean(depression),format="f", digits=4))
# A tibble: 3 x 2
  schooling_dummy `formatC(mean(depression), format = "f", digits = 4)`
            <dbl> <chr>                                                
1            0    10.4398                                              
2            1.00 9.7122                                               
3            2.00 10.1593    

现在，当我计算一个包含两个变量的模型时，为什么当两个组 = 0 时截距不是均值？回归**截距为 10.3796，但当性别 = 0，学校教育 = 0 时的平均值为 10.32548：

ds %>% group_by(schooling_dummy,sex_dummy) %>% 
+   summarise(formatC(mean(depression),format="f", digits=5))
# A tibble: 6 x 3
# Groups: schooling_dummy [?]
  schooling_dummy sex_dummy `formatC(mean(depression), format = "f", digits = 5)`
            <dbl>     <dbl> <chr>                                                
1            0         0    10.32548                                             
2            0         1.00 10.55404                                             
3            1.00      0    9.59305                                              
4            1.00      1.00 9.83139                                              
5            2.00      0    10.21218                                             
6            2.00      1.00 10.10648   

当我预测模型都为0时：

predict(mod3, data.frame(sex_dummy=0, schooling_dummy=0))
       1 
10.37956 

这个结果与抑郁症有关（当然......）但仍然不是我所期待的，因为： （参考：https://www.theanalysisfactor.com/interpret-the-intercept/）

这个previous forum post有什么相同之处 我知道我的变量是分类的，我正在调整我的脚本，因为您可以使用下面的代码重现： 谢谢

library(tidyverse)
set.seed(123)
ds <- data.frame(
  depression=rnorm(90,10,2),
  schooling_dummy=c(0,1,2),
  sex_dummy=c(0,1)
)
mod <- lm(data=ds, depression ~ relevel(factor(sex_dummy), ref = "0"))
summary(mod)
ds %>% group_by(sex_dummy) %>% 
  summarise(formatC(mean(depression),format="f", digits=4))

mod2 <- lm(data=ds, depression ~ relevel(factor(schooling_dummy), ref = "0"))
summary(mod2)
ds %>% group_by(schooling_dummy) %>% 
  summarise(formatC(mean(depression),format="f", digits=4))

mod3 <- lm(data=ds, depression ~ relevel(factor(sex_dummy), ref = "0") + 
             relevel(factor(schooling_dummy), ref = "0"))
summary(mod3)
ds %>% group_by(schooling_dummy,sex_dummy) %>% 
  summarise(formatC(mean(depression),format="f", digits=5))

predict(mod3, data.frame(sex_dummy=0, schooling_dummy=0))

Answer 1

您的想法中有两个错误（尽管您的 R 代码有效，所以这不是编程错误。

首先，你违反了你自己的声明，你有 not 虚拟编码 schooling 它不仅有零，而且它有 0,1 和 2。

二是你忘了lm建模中的交互效果...

试试这个...

library(tidyverse)
set.seed(123)
ds <- data.frame(
  depression=rnorm(90,10,2),
  schooling_dummy=c(0,1,2),
  sex_dummy=c(0,1)
)
# if you explicitly make these variables factors not integers R will do the right thing with them
ds$schooling_dummy<-factor(ds$schooling_dummy)
ds$sex_dummy<-factor(ds$sex_dummy)
ds %>% group_by(schooling_dummy,sex_dummy) %>%
    summarise(formatC(mean(depression),format="f", digits=5))
# you need an asterisk in your lm model to include the interaction term
lm(depression ~ schooling_dummy * sex_dummy, data = ds)

结果为您提供了预期的平均值...

Call:
lm(formula = depression ~ schooling_dummy * sex_dummy, data = ds)

Coefficients:
                (Intercept)             schooling_dummy1             schooling_dummy2  
                  10.325482                    -0.732433                    -0.113305  
                 sex_dummy1  schooling_dummy1:sex_dummy1  schooling_dummy2:sex_dummy1  
                   0.228561                     0.009778                    -0.334254  

如果您的数据以字符开头，那么您可以避免这种对分类变量的意外误用......因此，如果您的数据以这种方式编码：

ds <- data.frame(
  depression=rnorm(90,10,2),
  schooling=c("A","B","C"),
  sex=c("Male","Female")
) 

您犯同样错误的可能性较小，而且结果更易于阅读...