Skip to content

Instantly share code, notes, and snippets.

@bbolker
Created December 13, 2022 22:32
Show Gist options
  • Select an option

  • Save bbolker/e6694fef9cdc1e509379faed5b5dc39d to your computer and use it in GitHub Desktop.

Select an option

Save bbolker/e6694fef9cdc1e509379faed5b5dc39d to your computer and use it in GitHub Desktop.
plotting sjPlots in a grid
library(lme4)
library(sjPlot)
set.seed(101)
n <- 1e3 ## obs per group
ns <- 4 ## groups
dd <- data.frame(s = rep(1:ns, each= n),
f = factor(sample(1:20, size = n*ns, replace = TRUE)),
x = rnorm(n*ns))
sims <- lapply(1:ns,
function(i) simulate(~x + (1|f),
family = poisson,
newdata = subset(dd, s==i),
newparams = list(beta = c(1,1), theta = 1))[[1]])
dd$y <- unlist(sims)
fits <- lapply(1:ns,
function(i) glmer(y~x+(1|f),
family = poisson,
data = subset(dd, s==i)))
plots <- lapply(fits, function(f) plot_model(f, type = "pred")$x)
cowplot::plot_grid(plotlist = plots)
@strengejacke
Copy link

Not sure whether predictions or forest plots were requested? plot_models() usually plots forest plots only.

library(lme4)
#> Loading required package: Matrix
library(sjPlot)
#> #refugeeswelcome
library(parameters)

set.seed(101)
n <- 1e3 ## obs per group
ns <- 4  ## groups
dd <- data.frame(s = rep(1:ns, each= n),
                 f = factor(sample(1:20, size = n*ns, replace = TRUE)),
                 x = rnorm(n*ns))
sims <- lapply(1:ns,
               function(i) simulate(~x + (1|f),
                                    family = poisson,
                                    newdata = subset(dd, s==i),
                                    newparams = list(beta = c(1,1), theta = 1))[[1]])
dd$y <- unlist(sims)

fits <- lapply(1:ns,
               function(i) glmer(y~x+(1|f),
                                 family = poisson,
                                 data = subset(dd, s==i)))


# sjPlot
plot_models(fits)

# parameters
compare_parameters(fits, exponentiate = TRUE) |> plot()

# cowplot, but prediction
plots <- lapply(fits, function(f) plot_model(f, type = "pred")$x)
cowplot::plot_grid(plotlist = plots)

Created on 2022-12-15 with reprex v2.0.2

@drfeinberg
Copy link

I don't know how to really simulate this. I'm so sorry.
I have run plot_model on each subset and can make separate graphs, but I get lots of errors when trying this as plot_models to get 1 figure. If you need more info, please let me know.

emotions = c("ANG", "DIS", "FEA", "HAP", "NEU", "SAD")
sex_effects = c(-0.5, 0.5)

# Sample model.  The subsets change.. 6 emotions, 2 sexes.
sad_model_male_voices = glmer(data=subset(
    df, emotion=='SAD' & sex_effect==0.5
), 
  Trusted ~ zpitch + zcpp + zvtl_delta_f + zsubharmonic.to.harmonic.ratio + zspeechrate + zLTAS_slope +  zHNR +
    (1 | ActorID) + 
    (1 | participant), 
  family=binomial(link="logit"))



var_labels = c(
  zpitch="Pitch (mel)", 
  zcpp="Cepstral Peak Prominence", 
  zvtl_delta_f="Predicted Vocal Tract Length", 
  zsubharmonic.to.harmonic.ratio="Subharmonic to Harmonic Ratio",
  zspeechrate="Speech Rate", 
  zLTAS_slope = "Slope of Long-term Average Spectrum", 
  zHNR = "Harmonics to Noise Ratio"
  )

predictors <-c("zpitch", "zcpp", "zvtl_delta_f", "zsubharmonic.to.harmonic.ratio", "zspeechrate", "zLTAS_slope", "zHNR")



models <- c(anger_model_female_voices, disgust_model_female_voices, fear_model_female_voices, happy_model_female_voices, neutral_model_female_voices, sad_model_female_voices, anger_model_male_voices, disgust_model_male_voices, fear_model_male_voices, happy_model_male_voices, neutral_model_male_voices, sad_model_male_voices)



plot_models(
  models,
  grid=TRUE,
  vline.color = "red",
  axis.labels=var_labels,
  show.values = TRUE,
  value.offset = .3,
  sort.est = FALSE,
  transform = NULL,
  title = "Predictors of voice trust"
  )```

@strengejacke
Copy link

It's difficult to say without a reproducible example, but based on your code here and the error message you mentioned, can you try changing

models <- c(anger_model_female_voices, disgust_model_female_voices, fear_model_female_voices, happy_model_female_voices, neutral_model_female_voices, sad_model_female_voices, anger_model_male_voices, disgust_model_male_voices, fear_model_male_voices, happy_model_male_voices, neutral_model_male_voices, sad_model_male_voices)

into

models <- list(anger_model_female_voices, disgust_model_female_voices, fear_model_female_voices, happy_model_female_voices, neutral_model_female_voices, sad_model_female_voices, anger_model_male_voices, disgust_model_male_voices, fear_model_male_voices, happy_model_male_voices, neutral_model_male_voices, sad_model_male_voices)

I.e. replacing c() by list()?

@drfeinberg
Copy link

Thanks. Changing it to a list and commenting out all of the options that I set produced a plot. I can try to figure the errors out from here. Thanks so much!!!

@drfeinberg
Copy link

drfeinberg commented Dec 16, 2022

sort.est = TRUE or sort.est = FALSE causes the error "Warning: Could not access model information.Error in if (fam.info$is_linear) { : argument is of length zero"

@strengejacke
Copy link

Ok, I found the bug. sort.est is no valid argument for plot_models(), so it was captured by the ... and thereby considered as regression model object. I added a check for valid arguments and plot_models() now returns an informative message when this happens.

In short: plot_models(fits, sort.est = TRUE) is not intended to work.

There is a reason why there is no sort options, because estimates are grouped by models, and I'm assuming a case where you can also have more than 1 estimate. And I think you can't sort both by models/groups and within group by estimate in ggplot2.

library(lme4)
library(sjPlot)
library(parameters)

set.seed(101)
n <- 1e2 ## obs per group
ns <- 4  ## groups
dd <- data.frame(y = rpois(n*ns, 3),
                 f = factor(sample(1:20, size = n*ns, replace = TRUE)),
                 x1 = rnorm(n*ns),
                 x2 = runif(n*ns))
fits <- lapply(1:ns,
               function(i) glmer(y~x1+x2+(1|f),
                                 family = poisson,
                                 data = dplyr::sample_frac(dd, .6)))

plot_models(fits)

Created on 2022-12-17 with reprex v2.0.2

@strengejacke
Copy link

You can still use parameters::compare_parameters(fits) to get the raw data, but then probably need to build the plot from scratch.

@drfeinberg
Copy link

Thank you so much for all your help here!!!

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment