如何解决在不运行新模型的情况下从另一个参考水平获取模型估计值？
我想知道是否有一种简单的方法可以改变截距中的值，也许是数学上的方法，而无需重新运行大型模型。例如：
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,p-value: 1.14e-08
我想知道有没有一种方法可以在不重新运行模型的情况下获得这些值？您可以看到，每个模型（-5.96和5.96）中从4缸到6缸的比较是相同的，但是在两个模型中，我如何获得“其他”系数的估算值（例如{em> first 模型中的-2.55）。当然，在这种情况下，运行另一个模型只需要几分之一秒。但是对于非常大的模型，无需重新运行就可以更改参考级别。是否有相对简单的方法可以将所有估计值和标准误差转换为基于不同的参考水平，或者这样做太复杂了？
对于lme4，glmmTMB或rstanarm模型的任何解决方案都将受到赞赏。
解决方法
这是一个函数，它将为您提供给定因子变量的每次重排的系数，而无需再次运行模型或指定对比：
countPositives :: [Int] -> Int
countPositives xs = length [ x | x <- xs,x > 0]假设您有一个这样的模型：
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("^",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
我们需要一些说明，以了解其工作原理。
尽管上面的函数只运行一次模型，但让我们首先创建一个包含3个版本的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
的列表，其中mtcars的基线因子水平都不同。
cyl现在，我们可以使用df_list <- list(mtcars_4 = within(mtcars,cyl <- factor(cyl,c(4,6,8))),mtcars_6 = within(mtcars,c(6,4,mtcars_8 = within(mtcars,c(8,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等于0时，这些数字代表hp的期望值。
cyl由于我们现在拥有所有三个值的绝对值，因此将“ 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 <- 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
中的因子顺序是任意的：更改因子级别的顺序最终会给出相同的数值预测，即使摘要看起来不同。
TL; DR
因此，如果您只有第一个模型，那么从哪里得到-2.55的问题的答案是找出非截距系数之间的差异。在这种情况下
lm或者，将截距添加到非截距系数上，您可以看到如果任何水平是基线 ，截距将是什么，并且您可以获得相对于任何水平的任何水平的系数其他通过简单的减法。这就是我的功能(-8.520851) -(-5.967655)
#> [1] -2.553196
的工作方式。
由reprex package（v0.3.0）于2020-10-05创建
版权声明：本文内容由互联网用户自发贡献，该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容， 请发送邮件至 dio@foxmail.com 举报，一经查实，本站将立刻删除。
