Plotting Estimates (Fixed Effects) of Regression Models

Question

I'm using lmer4 package [lmer() function] to estimate several Average Models, which I want to plot their Estimated Coefficients. I found this document, "Plotting Estimates (Fixed Effects) of Regression Models, by Daniel Lüdecke" that explains how to plot Estimates, and it works with Average Models, but uses Conditional Average values instead of Full Average values.

Example of script:

library(lme4)
options(na.action = "na.omit")
PA_model_clima1_Om_ST <- lmer(O.matt ~ mes_N + Temperatura_Ar_PM_ST + RH_PM_ST + Vento_V_PM_ST + Evapotranspiracao_PM_ST + Preci_total_PM_ST + (1|ID), data=Abund)

library(MuMIn)
options(na.action = "na.fail")
PA_clima1_Om_ST<-dredge(PA_model_clima1_Om_ST)

sort.PA_clima1_Om_ST<- PA_clima1_Om_ST[order(PA_clima1_Om_ST$AICc),]
top.models_PA_clima1_Om_ST<-get.models(sort.PA_clima1_Om_ST, subset = delta < 2) 

model.sel(top.models_PA_clima1_Om_ST)
Avg_PA_clima1_Om_ST<-model.avg(top.models_PA_clima1_Om_ST, fit = TRUE) 
summary(Avg_PA_clima1_Om_ST)

Results of this script:

Term codes: 
Evapotranspiracao_PM_ST       Preci_total_PM_ST                RH_PM_ST    Temperatura_Ar_PM_ST 
                      1                       2                       3                       4 
          Vento_V_PM_ST 
                      5 

Model-averaged coefficients:  
(full average) 
                        Estimate Std. Error Adjusted SE z value Pr(>|z|)    
(Intercept)               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    
Evapotranspiracao_PM_ST  -0.1202     0.5174      0.5183   0.232 0.816654    
 
(conditional average) 
                        Estimate Std. Error Adjusted SE z value Pr(>|z|)    
(Intercept)               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    
Evapotranspiracao_PM_ST  -0.3053     0.7897      0.7912   0.386 0.699619    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Plot scrip:

library(sjPlot)
library(sjlabelled)
library(sjmisc)
library(ggplot2)
data(efc)
theme_set(theme_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")

Plot:

As you can see, it uses the values of Conditional Average values instead of Full Average values. How can I plot Estimates of Average Models using Full Average values?

Answer 1

I think it takes the conditional.. so unless you hack the function or maybe contact the author to have such an option, one way is to plot the coefficients yourself:

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

We take it out, pivot and plot:

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

You can also use parameters::model_parameters(), which is internally used by sjPlot::plot_model(). model_parameters() has a component-argument to decide which component to return. However, plot_model() does not yet pass additional arguments down to model_parameters(). I'm going to address this in sjPlot. Meanwhile, using model_parameters() at least offers a quick plot()-method.

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"))

You can do some minor modifications to the plot:

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

^{Created on 2020-06-27 by the reprex package (v0.3.0)}