时间戳

首页 > 解决方案 > 从另一个参考水平获取模型估计，而不运行新模型？

r - 从另一个参考水平获取模型估计，而不运行新模型？

问题描述

我想知道是否有一种简单的方法可以更改截距中的值，也许是数学上的，而无需重新运行大型模型。举个例子：

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 cyl，并且可以查看现在 8 cyl 与 6 cyl 而不是 4 cyl 的比较：

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 cyl 到 6 cyl 的比较在每个模型中都是相同的（-5.96和5.96），但是我将如何获得任一模型中“其他”系数的估计值（例如-2.55来自第一个模型）。当然，在这种情况下，运行另一个模型只需要几分之一秒。但是对于非常大的模型，无需重新运行即可更改参考电平会很方便。是否有相对简单的方法可以将所有估计值和标准误差转换为基于不同的参考水平，还是这样做太复杂了？

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

标签： rlme4coefficientsrstanarmglmmtmb

解决方案

这是一个函数，它将为您提供给定因子变量的每次重新排列的系数，而无需再次运行模型或指定对比：

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

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

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：改变因子水平的顺序最终给出相同的数值预测，即使总结看起来不同。

TL;博士

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

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

或者，将截距添加到非截距系数上，如果任何水平是基线，您可以看到截距是什么，并且您可以通过简单的减法获得相对于任何其他水平的任何水平的系数。这就是我的功能的rearrange_model_factors工作方式。

由reprex 包于 2020-10-05 创建(v0.3.0)

推荐阅读