如何对两个或多个因子自变量中的每个水平组合进行回归？

Question

我想知道如何对两个或多个因子变量的水平组合执行多个独立线性回归。

假设我们的数据集有一个因连续变量，然后是两个因子自变量和一个连续自变量。

那么假设我们在 r 中的回归公式是这样的：

model <- lm(weight ~ city + diet + height)

或者，要编写伪代码，我正在尝试这样做： lm(weight ~ height) %>% group by city lm(weight ~ height) %>% group by diet lm(weight ~ height) %>% group by city & diet

我知道我们可以对每个城市和饮食逐一进行线性回归，但您知道我们可以创建一个循环以便我们对数据集中的每个城市和饮食进行独立回归吗？

为了更好地说明这一点，我在这张图片中制作了这个假数据集，然后列出了我想要的三种类型的输出。但是，我不想一个一个地手动做，而是宁愿使用循环。

有谁知道如何在 r 中做到这一点？

Answer 1

首先定义一个小的regfun 来计算所需的汇总统计信息。然后，使用by 按组应用它。对于两个组的组合，我们可以paste 列在一起使用交互函数: 作为因子。

regfun <- function(x) summary(lm(w ~ h, x))$coe[2, c(1, 4)]

do.call(rbind, by(d, d$city, regfun))
#     Estimate  Pr(>|t|)
# a -0.1879530 0.4374580
# b -0.2143780 0.4674864
# c -0.2866948 0.5131854

do.call(rbind, by(d, d$diet, regfun))
#     Estimate  Pr(>|t|)
# y -0.1997162 0.3412652
# z -0.3512349 0.4312766

# do.call(rbind, by(d, Reduce(paste, d[1:2]), regfun))
with(d, do.call(rbind, by(d, city:diet, regfun)))  ## credits to @G.Grothendieck
#       Estimate  Pr(>|t|)
# a y -0.2591764 0.5576043
# a z -0.1543536 0.8158689
# b y -0.1966501 0.7485405
# b z -0.4354839 0.7461538
# c y -0.5000000 0.3333333
# c z -1.0671642 0.7221495

编辑

如果我们有一个不平衡的面板，即with(d, city:diet) 给出了实际上不在数据中的“不可能”组合，我们必须对其进行稍微不同的编码。您可以将by 视为第一个split 然后lapply 的组合，所以让我们开始吧。因为我们会得到错误，我们可以使用tryCatch 来提供类似的替代。

s <- with(d2, split(d2, city:diet))
do.call(rbind, lapply(s, function(x) 
  tryCatch(regfun(x), 
           error=function(e) cbind.data.frame(Estimate=NA, `Pr(>|t|)`=NA))))
#       Estimate  Pr(>|t|)
# a:y -0.2591764 0.5576043
# a:z         NA        NA
# b:y  5.2500000       NaN
# b:z         NA        NA
# c:y -0.5000000 0.3333333
# c:z  9.5000000       NaN
# d:y         NA        NA
# d:z  1.4285714       NaN
# e:y         NA        NA
# e:z -7.0000000       NaN
# f:y         NA        NA
# f:z  2.0000000       NaN

---

数据：

d <- structure(list(city = structure(c(1L, 2L, 3L, 1L, 2L, 3L, 1L, 
2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L), .Label = c("a", 
"b", "c"), class = "factor"), diet = structure(c(1L, 1L, 1L, 
2L, 2L, 2L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 2L, 2L, 2L), .Label = c("y", 
"z"), class = "factor"), id = c(1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 
2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L), w = c(66L, 54L, 50L, 
74L, 59L, 53L, 67L, 75L, 66L, 64L, 73L, 56L, 53L, 74L, 54L, 63L, 
69L, 75L), h = c(152L, 190L, 174L, 176L, 185L, 186L, 180L, 194L, 
154L, 169L, 183L, 177L, 189L, 152L, 182L, 191L, 173L, 179L)), out.attrs = list(
    dim = c(city = 3L, diet = 2L, id = 3L), dimnames = list(city = c("city=a", 
    "city=b", "city=c"), diet = c("diet=y", "diet=z"), id = c("id=1", 
    "id=2", "id=3"))), row.names = c(NA, -18L), class = "data.frame")

d2 <- structure(list(city = structure(c(1L, 2L, 3L, 4L, 5L, 6L, 1L, 
2L, 3L, 4L, 5L, 3L, 1L, 6L, 3L, 6L, 2L, 3L), .Label = c("a", 
"b", "c", "d", "e", "f"), class = "factor"), diet = structure(c(1L, 
1L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 2L, 2L, 
2L), .Label = c("y", "z"), class = "factor"), id = c(1L, 1L, 
1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L
), w = c(66L, 54L, 50L, 74L, 59L, 53L, 67L, 75L, 66L, 64L, 73L, 
56L, 53L, 74L, 54L, 63L, 69L, 75L), h = c(152L, 190L, 174L, 176L, 
185L, 186L, 180L, 194L, 154L, 169L, 183L, 177L, 189L, 152L, 182L, 
191L, 173L, 179L)), out.attrs = list(dim = c(city = 3L, diet = 2L, 
id = 3L), dimnames = list(city = c("city=a", "city=b", "city=c"
), diet = c("diet=y", "diet=z"), id = c("id=1", "id=2", "id=3"
))), row.names = c(NA, -18L), class = "data.frame")

Answer 2

我们可以在列表中定义模型规范，然后在所需模型列表上使用lapply()。

代码

models <- list("m1" = c("weight", "height"),
               "m2" = c("weight", "height", "city"),
               "m3" = c("weight", "height", "diet"),
               "m4" = c("weight", "height", "diet", "city"))

lapply(models, function(x){
  lm(weight ~ ., data = df[, x])
})

# $m1
# 
# Call:
# lm(formula = weight ~ ., data = df[, x])
# 
# Coefficients:
# (Intercept)       height  
#     -0.2970       0.1219  
#
#
# $m2
#
# Call:
# lm(formula = weight ~ ., data = df[, x])
#
# Coefficients:
# (Intercept)       height  cityHouston  
#     -0.3705       0.1259       0.1205  
#
#
# $m3
#
# Call:
# lm(formula = weight ~ ., data = df[, x])
# 
# Coefficients:
#    (Intercept)          height       dietVegan  dietVegetarian  
#        -0.1905          0.1270         -0.1288         -0.1757  
#
#
# $m4
#
# Call:
# lm(formula = weight ~ ., data = df[, x])
# 
# Coefficients:
#  (Intercept)          height       dietVegan  dietVegetarian     cityHouston  
#        -0.2615          0.1310         -0.1417         -0.1663          0.1197  

数据

df <- data.frame("weight" = rnorm(100), 
           "height" = rexp(100),
           "diet" = as.factor(sample(c("Vegan", "Vegetarian", "Meat"), replace = TRUE, 100)),
           "city" = as.factor(sample(c("Houston", "Chicago"), replace = TRUE, 100)))

Answer 3

为了让我的 cmets 进一步了解汇总分析相对于亚组分析的效率...

使用starwars 作为（不太理想的）起点：

d <- starwars %>% 
       filter(mass < 1000) %>%   # Exclude Jabba
       mutate(maleOrNot=ifelse(sex=="male", sex, "other")) %>% 
       replace_na(list(maleOrNot="other"))

为了论证，假设我们想仅根据角色是否为男性和身高来回归角色的质量，然后获得平均高度处预测质量的标准误差。

pData <- d %>% 
           group_by(maleOrNot) %>% 
           summarise(height=mean(height), .groups="drop")
pData

# A tibble: 2 x 2
  maleOrNot height
* <chr>      <dbl>
1 male        178.
2 other       162.

按组分析

lapply(
  d %>% pull(maleOrNot) %>% unique(),
  function(x) {
    m <- lm(mass ~ height, d %>% filter(maleOrNot == x))
    predict(m, pData %>% filter(maleOrNot ==  x), se.fit=TRUE)$se.fit
  }
)
[[1]]
[1] 2.656427

[[2]]
[1] 5.855176

现在汇总分析：

m <- lm(mass ~ maleOrNot + height, d)
predict(m, pData, se.fit=TRUE)$se.fit

       1        2 
2.789770 4.945734 

非男性的预测精度稍差 (5%)，但男性的预测精度提高了 15.5%。

但模型不是特别好。也许交互模型会有所改善：

m <- lm(mass ~ maleOrNot:height, d)
predict(m, pData, se.fit=TRUE)$se.fit
       1        2 
2.776478 4.880154 

现在这些数字分别下降了 4.5% 和 16.7%。在模型中包含其他项可能会进一步提高精度。

一般而言（尽管有例外），与拟合多个子组模型相比，拟合合并模型不太可能降低精度，并且可以显着提高精度。这是因为所有组都对（共同）方差的估计做出了贡献。

在计算时间方面：

library(microbenchmark)

byGroup <- function() {
  lapply(
    d %>% pull(maleOrNot) %>% unique(),
    function(x) {
      m <- lm(mass ~ height, d %>% filter(maleOrNot == x))
      predict(m, pData %>% filter(maleOrNot ==  x), se.fit=TRUE)$se.fit
    }
  )
}

pooled <- function() {
  m <- lm(mass ~ maleOrNot + height, d)
  predict(m, pData, se.fit=TRUE)$se.fit
}

microbenchmark(byGroup, pooled, times=100)
Unit: nanoseconds
    expr min   lq  mean median uq  max neval
 byGroup  44 45.5 55.22     47 48  891   100
  pooled  42 44.0 60.27     46 47 1434   100

所以对于这个简单的案例，几乎没有区别。更复杂的例子可能会给出不同的答案。