0

I am trying to create a plot like the following (taken from https://cran.r-project.org/web/packages/dotwhisker/vignettes/dotwhisker-vignette.html), but with no intercept, and with no rescaling:

I have been following the code from https://cran.r-project.org/web/packages/glmmTMB/vignettes/model_evaluation.pdf:

if (requireNamespace("broom.mixed") && requireNamespace("dotwhisker")) {
(t1 <- broom.mixed::tidy(owls_nb1, conf.int = TRUE))
if (packageVersion("dotwhisker")>"0.4.1") {
## to get this version (which fixes various dotwhisker problems)
## use devtools::install_github("bbolker/broom.mixed") or
## wait for pull request acceptance/submission to CRAN/etc.
dwplot(owls_nb1)+geom_vline(xintercept=0,lty=2)
} else {
owls_nb1$coefficients <- TRUE ## hack!
dwplot(owls_nb1,by_2sd=FALSE)+geom_vline(xintercept=0,lty=2)
}
} 

I have tried to apply this code in the following way to achieve a dwplot for multiple models (with no intercepts and without rescaling):

mod1 <- glmmTMB(Species1_Occupancy ~ offset(log(Study_Unit_Area)) + (1|Site) + scale(Var1) + scale(Var2) + scale(Var3) + scale(Var4) + scale(Var5) + scale(Var6), data=data, family=binomial(link="cloglog"))

mod2 <- glmmTMB(Species2_Occupancy ~ offset(log(Study_Unit_Area)) + (1|Site) + scale(Var1) + scale(Var2) + scale(Var3) + scale(Var4) + scale(Var5) + scale(Var6), data=data, family=binomial(link="cloglog")) 

mod3 <- glmmTMB(Species3_Occupancy ~ offset(log(Study_Unit_Area)) + (1|Site) + scale(Var1) + scale(Var2) + scale(Var3) + scale(Var4) + scale(Var5) + scale(Var6), data=data, family=binomial(link="cloglog"))      

dwplot(list(mod1,mod2,mod3), show_intercept = FALSE, by_2sd = FALSE)

But I get this error message:

Error in tidy.glmmTMB(., conf.int = TRUE, ...) : 
  unknown effect type list(par = c(beta = 0, beta = 0, beta = 0, beta = 0, beta = 0, beta = 0, beta = 0, theta = 0), fn = function (x = last.par[-random], ...) 
{
    if (tracepar) {
        cat("par:\n")
        print(x)
    }
    if (!validpar(x)) 
        return(NaN)
    ans <- try({
        if (MCcontrol$doMC) {
        ff(x, order = 0)
        MC(last.par, n = MCcontrol$n, seed = MCcontrol$seed, order = 0)
    }
    else ff(x, order = 0)
}, silent = silent)
if (is.character(ans)) 
    NaN
else ans
}, gr = function (x = last.par[-random], ...) 
{
ans <- try({
    if (MCcontrol$doMC) {
        ff(x, order = 0)
        MC(last.par, n = MCcontrol$n, seed = MCcontrol$seed, order = 1)
        }
        else ff(x, order = 1)
    }, silent = silent)
    if (is.character(ans)) 
        ans <- rep(NaN, length(x))
    if (tracemgc) 
        cat("outer mgc: ", max(abs(ans)), "\n")
    ans
}, he = function (x = last.par[-random], ...) 
{

When I run this without the by_2sd() function, I get no error message, but the intercept is still displayed.

dwplot(list(mod1,mod2,mod3), show_intercept = FALSE)

Any help would be much appreciated.

MoriahT
  • 1
  • 2
  • You could try [`sjPlot::plot_models()`](https://strengejacke.github.io/sjPlot/reference/plot_models.html). Does this work? – Daniel May 11 '20 at 19:35
  • @Daniel I also had no luck with plot_models(). On running: `plot_models(mod1, mod2, mod3, mod4, grid = TRUE)` I got the following error: `Error in is.null(transform) && fitfam$is_multinomial : invalid 'y' type in 'x && y'` Perhaps this is due to the model type I am running...they are binomial GLMMs with the clog-log link, run using the mixed model package glmmTMB. The response variable is made up of 1/0 for presence/absence of species – MoriahT May 20 '20 at 20:28
  • This error sounds like you need to update the _insight_ package. Does this help? – Daniel May 21 '20 at 17:29
  • Yes, the _plot_models()_ function worked perfectly following update of the _insight_ package. Thanks very much – MoriahT May 21 '20 at 21:11
  • @Daniel, a followup question about _plot_models_: My input models use the clog-log link (e.g. `glmmTMB(species_occupancy ~ scale(var1) + scale(var2), family=binomial(link="cloglog"))`) and I have plotted them using the following code: `plot_models(mod1, mod2, mod3, mod4, grid = TRUE, transform = NULL)`, but the output plot refers to the values as log-odds. Should the values not be referred to as the percentage change in the probability[Y = 1] for a one unit change of a predictor variable? – MoriahT May 21 '20 at 22:18
  • That's an issue in the current CRAN version of sjPlot, but should have been fixed on the GitHub version of sjPlot. – Daniel May 22 '20 at 07:35
  • An update will be submitted to CRAN the next days. – Daniel May 22 '20 at 07:36
  • Yes that works, very happy with the plot. _sjPlot_ is more intuitive than _dotwhisker_. Thanks very much. – MoriahT May 22 '20 at 16:18

0 Answers0