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)
@drfeinberg
Copy link
Copy Markdown

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
Copy Markdown

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
Copy Markdown

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
Copy Markdown

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
Copy Markdown

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
Copy Markdown

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
Copy Markdown

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