Dyl*_*mes 5 r lme4 coefficients rstanarm glmmtmb
我想知道是否有一种简单的方法可以在不重新运行大型模型的情况下改变截距中的值，也许是数学上的。举个例子：
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缸到6缸的比较在每种型号（同-5.96和5.96），但我怎么会得到在任何模型中的“其他”系数估计值（例如-2.55从第一个模型）。当然，在这种情况下，运行另一个模型只需要几分之一秒。但是对于非常大的模型，无需重新运行即可更改参考电平会很方便。是否有相对简单的方法可以将所有估计值和标准误差转换为基于不同的参考水平，或者做这样的事情是否太复杂？
对于lme4、glmmTMB或rstanarm模型的任何解决方案将不胜感激。
这是一个函数，可以为您提供给定因子变量的每次重新排列的系数，而无需再次运行模型或指定对比：
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(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
由于我们现在将所有三个值都设为绝对值，因此我们将“Intercept”重命名为适当的因子级别：
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
最后，让我们重新排列顺序，以便我们可以比较所有三个因子水平的绝对值：
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 的问题的答案是找到非截距系数之间的差异。在这种情况下
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
或者，将截距添加到非截距系数上，您可以看到如果任何水平为基线的话截距会是多少，并且您可以通过简单的减法获得任何水平相对于任何其他水平的系数。这就是我的功能的rearrange_model_factors工作原理。
由reprex 包于 2020 年 10 月 5 日创建(v0.3.0)