# Effect size and power analysis for zero-inflated negative binomial mixed models

I'm working with a zero-inflated dataset (confirmed using vcdExtra::zero.test) of bat calls per unit time (i.e. bat activity), which I'm relating to the occurrence of passenger trains passing the recording site.

To assess the effect of passing trains, I compared the number of bat calls in the 30 seconds before, after and between train passes using a zero-inflated negative binomial mixed model (glmmTMB package) with the structure:

glmmTMB(bat call rate ~ category + (1|site),
ziformula = ~.,
family = nbinom2)


...where category is either before, after or between, and (1|site) accounts for repeated measures at each location.

This showed the number of bat calls was lower in the after category (i.e. there was less activity in the 30 seconds after a train passed).

To examine if the reduction in calls lasts longer than 30 seconds, I re-ran the model comparing bat call rates from time intervals of increasing duration - 60/120/240 seconds - divided into the same before, after and between train categories. And, category was no longer significant in any of these models.

All intervals between trains are ≥ 90 seconds in duration, with frequency decreasing as duration increases. This means sample size reduced sequentially from 3000 in the initial 30 second model, to 1800 in the 60 second model, 1000 in the 120 second model, and 500 in the 240 second model.

I'd like to determine whether the lack of relationship in the 60/120/240 seconds models is because there is only a reduction in bat calls for 30 seconds, or is due to their lower sample sizes.

I assume (perhaps naively) I could do this by comparing the effect size for the 30 seconds model with the power of the 60/120/240 seconds models. But, I'm not sure how to calculate these values (it's my first time using zero-inflated models).

According to Johnson et al. (2015), it is possible, but there's no information in the paper about how to do it, and GLMMmisc::sim.glmm doesn't seem to work with glmmTMB models.

Any advice, or ideas for alternative approaches would be greatly appreciated!

In general, performing a power analysis for your own study is not a good idea; see here from more info. The information you have in your own data is summarized in the corresponding confidence intervals.

If you still like to do a power analysis for a future study using the zero-inflated negative binomial model, then it is indeed the best way to do it via simulation. The following code simulates from a simple version of this model. You could adapt it to your needs:

set.seed(123)
n <- 300 # number of subjects
K <- 8 # number of measurements per subject
t_max <- 5 # maximum follow-up time

# we constuct a data frame with the design:
# everyone has a baseline measurment, and then measurements at random follow-up times
DF <- data.frame(id = rep(seq_len(n), each = K),
time = c(replicate(n, c(0, sort(runif(K - 1, 0, t_max))))),
sex = rep(gl(2, n/2, labels = c("male", "female")), each = K))

# design matrices for the fixed and random effects non-zero part
X <- model.matrix(~ sex * time, data = DF)
Z <- model.matrix(~ 1, data = DF)
# design matrices for the fixed and random effects zero part
X_zi <- model.matrix(~ sex, data = DF)
Z_zi <- model.matrix(~ 1, data = DF)

betas <- c(1.5, 0.05, 0.05, -0.03) # fixed effects coefficients non-zero part
shape <- 2 # shape/size parameter of the negative binomial distribution
gammas <- c(-1.5, 0.5) # fixed effects coefficients zero part
D11 <- 0.5 # variance of random intercepts non-zero part
D22 <- 0.4 # variance of random intercepts zero part

# we simulate random effects
b <- cbind(rnorm(n, sd = sqrt(D11)), rnorm(n, sd = sqrt(D22)))
# linear predictor non-zero part
eta_y <- as.vector(X %*% betas + rowSums(Z * b[DF$$id, 1, drop = FALSE])) # linear predictor zero part eta_zi <- as.vector(X_zi %*% gammas + rowSums(Z_zi * b[DF$$id, 2, drop = FALSE]))
# we simulate negative binomial longitudinal data
DF$$y <- rnbinom(n * K, size = shape, mu = exp(eta_y)) # we set the extra zeros DF$$y[as.logical(rbinom(n * K, size = 1, prob = plogis(eta_zi)))] <- 0