不需要重新运行模型，如何从另一个参考水平获取模型估计值？

Question

我想知道是否有一种简单的方法来改变截距中的值，可能是数学上的方式，而不需要重新运行大模型。举个例子：

mtcars$cyl<-as.factor(mtcars$cyl)

summary(
  lm(mpg~cyl+hp,data=mtcars)
)

输出：

    Call:
lm(formula = mpg ~ cyl + hp, data = mtcars)

Residuals:
   Min     1Q Median     3Q    Max 
-4.818 -1.959  0.080  1.627  6.812 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 28.65012    1.58779  18.044  < 2e-16 ***
cyl6        -5.96766    1.63928  -3.640  0.00109 ** 
cyl8        -8.52085    2.32607  -3.663  0.00103 ** 
hp          -0.02404    0.01541  -1.560  0.12995    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.146 on 28 degrees of freedom
Multiple R-squared:  0.7539,    Adjusted R-squared:  0.7275 
F-statistic: 28.59 on 3 and 28 DF,  p-value: 1.14e-08

现在我可以将参考电平更改为6缸，并且可以看到8缸现在与6缸的比较，而不是4缸。

mtcars$cyl<-relevel(mtcars$cyl,"6")

summary(
  lm(mpg~cyl+hp,data=mtcars)
)

输出：

Call:
lm(formula = mpg ~ cyl + hp, data = mtcars)

Residuals:
   Min     1Q Median     3Q    Max 
-4.818 -1.959  0.080  1.627  6.812 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 22.68246    2.22805   10.18 6.48e-11 ***
cyl4         5.96766    1.63928    3.64  0.00109 ** 
cyl8        -2.55320    1.97867   -1.29  0.20748    
hp          -0.02404    0.01541   -1.56  0.12995    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.146 on 28 degrees of freedom
Multiple R-squared:  0.7539,    Adjusted R-squared:  0.7275 
F-statistic: 28.59 on 3 and 28 DF,  p-value: 1.14e-08

我想知道是否有一种方法能够在不重新运行模型的情况下获取这些值？您可以看到，在每个模型中，从4缸到6缸的比较是相同的(-5.96和5.96)，但是如何获得任一模型中“其他”系数的估计值(例如第一个模型中的-2.55)。当然，在这种情况下，运行其他模型只需要几分之一秒钟的时间。但是对于非常大的模型，能够在不重新运行模型的情况下更改参考水平将是方便的。是否有相对简单的方法可以将所有估计值和标准误转换为基于不同参考水平的值，还是这样做过于复杂？

对于lme4、glmmTMB或rstanarm模型的任何解决方案将不胜感激。

- Dylan_Gomes

---

请查看 emmeans 包。 - Ben Bolker

还有multcomp::glht，它可以让您手动构建查询。 - user20650

可能有帮助的链接：https://stats.stackexchange.com/questions/168650/how-to-set-custom-contrasts-with-lmer-in-r - Mike

Answer 1

这里有一个函数，可以为给定因子变量的每个重新排列提供系数，而无需再次运行模型或指定对比方式：

rearrange_model_factors <- function(model, var)
{
  var   <- deparse(substitute(var))
  coefs <- coef(model)
  level_coefs <- grep(paste0("^", var), names(coefs))
  coefs[level_coefs] <- coefs[1] + coefs[level_coefs]
  used_levels <- gsub(var, "", names(coefs[level_coefs]))
  all_levels <- levels(model$model[[var]])
  names(coefs)[1] <- paste0(var, setdiff(all_levels, used_levels))
  level_coefs <- grep(paste0("^", var), names(coefs))
  levs <- coefs[level_coefs]
  perms <- gtools::permutations(length(levs), length(levs))
  perms <- lapply(seq(nrow(perms)), function(i) levs[perms[i,]])

  lapply(perms, function(x) { 
    x[-1] <- x[-1] - x[1]
    coefs[level_coefs] <- x
    names(coefs)[level_coefs] <- names(x)
    names(coefs)[1] <- "(Intercept)"
    coefs
    })
}

假设您有这样的一个模型：

iris_mod <- lm(Sepal.Width ~ Species + Sepal.Length, data = iris)

为了查看如果Species顺序不同，您的系数会如何改变，您只需要执行以下操作：

rearrange_model_factors(iris_mod, Species)
#> [[1]]
#>       (Intercept) Speciesversicolor  Speciesvirginica      Sepal.Length 
#>         1.6765001        -0.9833885        -1.0075104         0.3498801 
#> 
#> [[2]]
#>       (Intercept)  Speciesvirginica Speciesversicolor      Sepal.Length 
#>         1.6765001        -1.0075104        -0.9833885         0.3498801 
#> 
#> [[3]]
#>      (Intercept)    Speciessetosa Speciesvirginica     Sepal.Length 
#>       0.69311160       0.98338851      -0.02412184       0.34988012 
#> 
#> [[4]]
#>      (Intercept) Speciesvirginica    Speciessetosa     Sepal.Length 
#>       0.69311160      -0.02412184       0.98338851       0.34988012 
#> 
#> [[5]]
#>       (Intercept)     Speciessetosa Speciesversicolor      Sepal.Length 
#>        0.66898976        1.00751035        0.02412184        0.34988012 
#> 
#> [[6]]
#>       (Intercept) Speciesversicolor     Speciessetosa      Sepal.Length 
#>        0.66898976        0.02412184        1.00751035        0.34988012 

或者使用您自己的示例：

mtcars$cyl <- as.factor(mtcars$cyl)

rearrange_model_factors(lm(mpg ~ cyl + hp, data = mtcars), cyl)
#> [[1]]
#> (Intercept)        cyl6        cyl8          hp 
#> 28.65011816 -5.96765508 -8.52085075 -0.02403883 
#> 
#> [[2]]
#> (Intercept)        cyl8        cyl6          hp 
#> 28.65011816 -8.52085075 -5.96765508 -0.02403883 
#> 
#> [[3]]
#> (Intercept)        cyl4        cyl8          hp 
#> 22.68246309  5.96765508 -2.55319567 -0.02403883 
#> 
#> [[4]]
#> (Intercept)        cyl8        cyl4          hp 
#> 22.68246309 -2.55319567  5.96765508 -0.02403883 
#> 
#> [[5]]
#> (Intercept)        cyl4        cyl6          hp 
#> 20.12926741  8.52085075  2.55319567 -0.02403883 
#> 
#> [[6]]
#> (Intercept)        cyl6        cyl4          hp 
#> 20.12926741  2.55319567  8.52085075 -0.02403883 

---

我们需要一点解释才能看出这个是如何起作用的。

虽然上面这个函数只运行了一次模型，但让我们从创建一个包含3个版本的mtcars列表开始，其中cyl的基准因子水平都不同。

df_list <- list(mtcars_4 = within(mtcars, cyl <- factor(cyl, c(4, 6, 8))),
                mtcars_6 = within(mtcars, cyl <- factor(cyl, c(6, 4, 8))),
                mtcars_8 = within(mtcars, cyl <- factor(cyl, c(8, 4, 6))))

现在我们可以使用 lapply 一次性提取您模型的所有三个版本的系数。为了清晰起见，我们将删除 hp 系数，因为它在所有三个版本中保持不变：

coefs <- lapply(df_list, function(x) coef(lm(mpg ~ cyl + hp, data = x))[-4])

coefs
#> $mtcars_4
#> (Intercept)        cyl6        cyl8 
#>   28.650118   -5.967655   -8.520851 
#> 
#> $mtcars_6
#> (Intercept)        cyl4        cyl8 
#>   22.682463    5.967655   -2.553196 
#> 
#> $mtcars_8
#> (Intercept)        cyl4        cyl6 
#>   20.129267    8.520851    2.553196

现在，我们需要提醒自己，每个因素水平的系数都是相对于基线水平给出的。这意味着对于非截距系数，我们可以简单地将截距值加上它们的系数以获得它们的绝对值。这意味着这些数字表示当hp等于cyl的三个级别时mpg的期望值为0。

coefs <- lapply(coefs, function(x) c(x[1], x[-1] + x[1]))

coefs
#> $mtcars_4
#> (Intercept)        cyl6        cyl8 
#>    28.65012    22.68246    20.12927 
#> 
#> $mtcars_6
#> (Intercept)        cyl4        cyl8 
#>    22.68246    28.65012    20.12927 
#> 
#> $mtcars_8
#> (Intercept)        cyl4        cyl6 
#>    20.12927    28.65012    22.68246

由于我们现在的所有三个值都是绝对值，让我们将“Intercept”重命名为适当的因子级别：

coefs <- mapply(function(x, y) { names(x)[1] <- y; x}, 
                x = coefs, y = c("cyl4", "cyl6", "cyl8"), SIMPLIFY = FALSE)
coefs
#> $mtcars_4
#>     cyl4     cyl6     cyl8 
#> 28.65012 22.68246 20.12927 
#> 
#> $mtcars_6
#>     cyl6     cyl4     cyl8 
#> 22.68246 28.65012 20.12927 
#> 
#> $mtcars_8
#>     cyl8     cyl4     cyl6 
#> 20.12927 28.65012 22.68246

最后，让我们重新排列顺序，以便比较所有三个因子水平的绝对值：

coefs <- lapply(coefs, function(x) x[order(names(x))])

coefs
#> $mtcars_4
#>     cyl4     cyl6     cyl8 
#> 28.65012 22.68246 20.12927 
#> 
#> $mtcars_6
#>     cyl4     cyl6     cyl8 
#> 28.65012 22.68246 20.12927 
#> 
#> $mtcars_8
#>     cyl4     cyl6     cyl8 
#> 28.65012 22.68246 20.12927

我们可以看到它们都是相同的。这就是为什么在lm中因子的排序是任意的：即使摘要看起来不同，改变因子级别的顺序最终会得出相同的数值预测结果。

---

简而言之

所以回答你的问题，如果你只有第一个模型，那么你从哪里得到-2.55的答案是找到非截距系数之间的差异。 在这种情况下。

(-8.520851) -(-5.967655) 
#> [1] -2.553196

或者，将拦截器添加到非拦截器系数上，您可以看到如果任何水平是基线，则插值将是什么，并且您可以通过简单的减法获取相对于任何其他级别的系数。这就是我的rearrange_model_factors函数的工作原理。

^{由reprex软件包 (v0.3.0)于2020-10-05创建}

- Allan Cameron

---

@BenBolker，我并不是在建议这个。我只是展示了因子不同排序的等价性。OP问道“但是我该如何获得任一模型中'其他'系数的估计值（例如第一个模型中的-2.55）。” - 我正在展示如何做到这一点。 - Allan Cameron

@BenBolker 我的意思是如何在不重新运行模型的情况下检索不同顺序因子的系数。抱歉我的表述不够清晰。 - Allan Cameron

你正在展示如何手动完成它。我认为OP正在寻找机械化的方法... - Ben Bolker

这是如何计算估计值的详细解释 - 感谢您。计算SE/CI是否同样简单？我猜不是。当然，我的实际模型是一个更复杂的rstanarm模型，具有交互项的随机斜率，因此我正在尝试找出一种手动绘制估计值和CI的方法。不幸的是，我很难理解如何从截距中提取一些术语。幸运的是，交互作用是与2级因子相关的，因此我只需要运行两个模型即可获得估计值和不确定性。 - Dylan_Gomes

@Dylan_Gomes 不，SE和CI更难。然而，现在你知道只需要将截距加到因子系数上即可获得每个水平的绝对值，而SE将保持不变 - 你可以使用这个事实为每个水平计算CI。 - Allan Cameron