绘制回归模型的估计值（固定效应）

Question

我正在使用 lmer4 包 [lmer() 函数] 来估计几个平均模型，我想绘制它们的估计系数。我找到了这篇文档“Plotting Estimates (Fixed Effects) of Regression Models, by Daniel Lüdecke”，它解释了如何绘制 Estimates，它适用于平均模型，但使用条件平均值而不是完全平均值。

脚本示例：

库（lme4）
选项（na.action = “na.omit”）
PA_model_clima1_Om_ST 
此脚本的结果：
术语代码：
Evapotranspiracao_PM_ST Preci_total_PM_ST RH_PM_ST Temperatura_Ar_PM_ST
                      1 2 3 4
          Vento_V_PM_ST
                      5

模型平均系数：
（全平均）
                        估计标准。误差调整 SE z 值 Pr(>|z|)
(截取) 5.4199 1.4094 1.4124 3.837 0.000124 ***
Preci_total_PM_ST -0.8679 1.0300 1.0313 0.842 0.400045
RH_PM_ST 0.6116 0.8184 0.8193 0.746 0.455397
Temperatura_Ar_PM_ST -1.9635 0.7710 0.7725 2.542 0.011026 *
Vento_V_PM_ST -0.6214 0.7043 0.7052 0.881 0.378289
蒸散量_PM_ST -0.1202 0.5174 0.5183 0.232 0.816654
 
（条件平均）
                        估计标准。误差调整 SE z 值 Pr(>|z|)
(截取) 5.4199 1.4094 1.4124 3.837 0.000124 ***
Preci_total_PM_ST -1.2200 1.0304 1.0322 1.182 0.237249
RH_PM_ST 1.0067 0.8396 0.8410 1.197 0.231317
Temperatura_Ar_PM_ST -1.9635 0.7710 0.7725 2.542 0.011026 *
Vento_V_PM_ST -0.8607 0.6936 0.6949 1.238 0.215546
蒸散量_PM_ST -0.3053 0.7897 0.7912 0.386 0.699619
---
意义。代码：0‘***’0.001‘**’0.01‘*’0.05‘.’0.1‘’1

剧情脚本：
库（sjPlot）
图书馆（sjlabeled）
图书馆（sjmisc）
图书馆（ggplot2）
数据（EFC）
主题集（主题_sjplot（））

plot_model(Avg_PA_clima1_Om_ST, type="est", vline.color="black", sort.est = TRUE, show.values = TRUE, value.offset = .3, title="O. mattogrossae")

剧情：

如您所见，它使用条件平均值而不是完全平均值的值。
如何使用完全平均值绘制平均模型的估计值？
    

Answer 1

我认为它需要条件..所以除非你破解函数或者联系作者以获得这样的选择，一种方法是自己绘制系数：

library(lme4)
library(MuMIn)
options(na.action = "na.fail")
set.seed(888)
dat= data.frame(y = rnorm(100),
var1 = rnorm(100),var2 = rnorm(100),
var3=rnorm(100),rvar = sample(1:2,replace=TRUE,100))

lme_mod <- lmer(y ~ var1+ var2+ var3 + (1|rvar), dat)
dre_mod <- dredge(lme_mod)
avg_mod = model.avg(dre_mod,fit=TRUE) 

summary(avg_mod)

Model-averaged coefficients:  
(full average) 
            Estimate Std. Error Adjusted SE z value Pr(>|z|)
(Intercept) -0.02988    0.18699     0.18936   0.158    0.875
var2        -0.03791    0.08817     0.08858   0.428    0.669
var1        -0.02999    0.07740     0.07778   0.386    0.700
var3         0.01521    0.05371     0.05404   0.281    0.778
 
(conditional average) 
            Estimate Std. Error Adjusted SE z value Pr(>|z|)
(Intercept) -0.02988    0.18699     0.18936   0.158    0.875
var2        -0.16862    0.11197     0.11339   1.487    0.137
var1        -0.15293    0.10841     0.10978   1.393    0.164
var3         0.11227    0.10200     0.10327   1.087    0.277

The matrix is under:

    summary(avg_mod)$coefmat.full
                    Estimate Std. Error Adjusted SE   z value  Pr(>|z|)
(Intercept) -0.02988418 0.18698720  0.18935677 0.1578194 0.8745991
var2        -0.03791016 0.08816936  0.08857788 0.4279867 0.6686608
var1        -0.02998709 0.07740247  0.07778028 0.3855360 0.6998404
var3         0.01520633 0.05371407  0.05404100 0.2813850 0.7784151

我们将其取出、旋转并绘制：

library(ggplot2)

df = data.frame(summary(avg_mod)$coefmat.full)
df$variable = rownames(df)
colnames(df)[2] = "std_error"
df = df[df$variable !="(Intercept)",]
df$type = ifelse(df$Estimate>0,"pos","neg")

ggplot(df,aes(x=variable,y=Estimate))+
geom_point(aes(col=type),size=3) + 
geom_errorbar(aes(col=type,ymin=Estimate-1.96*std_error,ymax=Estimate+1.96*std_error),width=0,size=1) + 
geom_text(aes(label=round(Estimate,digits=2)),nudge_x =0.1) +
geom_hline(yintercept=0,col="black")+ theme_bw()+coord_flip()+
scale_color_manual(values=c("#c70039","#111d5e")) +
theme(legend.position="none")

Answer 2

也可以使用parameters::model_parameters()，sjPlot::plot_model() 内部使用。 model_parameters() 有一个 component 参数来决定返回哪个组件。但是，plot_model() 尚未将其他参数传递给model_parameters()。我将在 sjPlot 中解决这个问题。同时，使用model_parameters() 至少提供了一个快速的plot() 方法。

library(lme4)
library(MuMIn)
options(na.action = "na.fail")
set.seed(888)
dat= data.frame(y = rnorm(100),
                var1 = rnorm(100),var2 = rnorm(100),
                var3=rnorm(100),rvar = sample(1:2,replace=TRUE,100))

lme_mod <- lmer(y ~ var1+ var2+ var3 + (1|rvar), dat)
dre_mod <- dredge(lme_mod)
xavg_mod = model.avg(dre_mod,fit=TRUE) 

library(parameters)
model_parameters(avg_mod)
#> Parameter   | Coefficient |   SE |        95% CI |    z | df |     p
#> --------------------------------------------------------------------
#> (Intercept) |       -0.03 | 0.19 | [-0.40, 0.34] | 0.16 | 96 | 0.875
#> var2        |       -0.17 | 0.11 | [-0.39, 0.05] | 1.49 | 96 | 0.137
#> var1        |       -0.15 | 0.11 | [-0.37, 0.06] | 1.39 | 96 | 0.164
#> var3        |        0.11 | 0.10 | [-0.09, 0.31] | 1.09 | 96 | 0.277

model_parameters(avg_mod, component = "full")
#> Parameter   | Coefficient |   SE |        95% CI |    z | df |     p
#> --------------------------------------------------------------------
#> (Intercept) |       -0.03 | 0.19 | [-0.40, 0.34] | 0.16 | 96 | 0.875
#> var2        |       -0.04 | 0.09 | [-0.21, 0.14] | 0.43 | 96 | 0.669
#> var1        |       -0.03 | 0.08 | [-0.18, 0.12] | 0.39 | 96 | 0.700
#> var3        |        0.02 | 0.05 | [-0.09, 0.12] | 0.28 | 96 | 0.778

plot(model_parameters(avg_mod, component = "full"))

你可以对情节做一些小的修改：

library(ggplot2)
plot(model_parameters(avg_mod, component = "full")) +
  geom_text(aes(label = round(Coefficient, 2)), nudge_x = .2)

^{由reprex package (v0.3.0) 于 2020 年 6 月 27 日创建}