pitakakariki / simr

Power Analysis of Generalised Linear Mixed Models by Simulation
70 stars 19 forks source link

powerSim fails with neg.binom model created with makeGlmer #241

Open Ben-Cox opened 2 years ago

Ben-Cox commented 2 years ago

Running powerSim() with a glmer.nb created with makeGlmer throws errors in each simulation. The error is: object 'Negative Binomial()' of mode 'function' was not found.

Reprex:

library(tidyverse)
library(simr)

n <- 500

x <- sample(c("a","b","c"), n, replace=TRUE)
y <- sample(c("w","x","y","z"), n, replace=TRUE)

d <- tibble(x, y)

mod <- makeGlmer(count ~ x + (1|y), 
                 family = negative.binomial(theta=0.1), 
                 fixef = c(-1.2,0.5,1.2), 
                 VarCorr = 0.75^2, 
                 data=d)

summary(mod)

sim <- powerSim(mod, fixed("x", "z"), nsim=10)  

sim$errors
pitakakariki commented 2 years ago

Hmmm makeGlmer isn't supposed to work with negative binomials so it probably shouldn't have let you get this far without throwing a helpful error message.

Ben-Cox commented 2 years ago

Ahh I did not realize it isn't supposed to work with negative binomials. Thanks for the heads up.

pitakakariki commented 2 years ago

You could probably hack something together - simulate from a Poisson, fit a negative binomial to that data, then use the barely-tested negative binomial features in powerSim.

Ben-Cox commented 2 years ago

I've got a pretty hacky workaround using lme4 to simulate, fitting to the simulated data, then fixing the betas and RE variances on the model and passing that to powerSim. When I come at it that way power sim works...is this a bad idea for statistical reasons? I got similar results brute forcing the simulations, fitting, and checking the Pr(>|z|) value.

pitakakariki commented 2 years ago

Should be fine - just be aware that glmer.nb support is experimental so if the results don't make sense you might have found a bug.

Statistically, you probably want to adjust the overdispersion parameter as well as the usual beta and RE parameters - although I can't remember if that's possible. Might be able to approximate that by simulating from a "NegBin w/ known Overdispersion" model.

Ben-Cox commented 2 years ago

Ok sounds good. I forgot to mention I fixed the neg.binomal theta to the value estimated from pilot data as well.

Thanks again.