We conducted a behavior change field experiment using the following variables:
- Two time points (T0, T1)
- Two groups (intervention vs. control)
- Individual ID (n = 62 in group 1, n = 53 in group 2)
- Workshop ID (8 workshops)
We are analyzing the data in lme4 using this model
lmer(DV_T1 ~ 1 + group + DV_T0 + (1|workshop_id), data=df)
Sample size was given by resource constraints which is why we decided to run a sensitivity power analysis to detect the minimum effect size of interest (https://lakens.github.io/statistical_inferences/08-samplesizejustification.html).
We tried doing this using the simr package, but we am not entirely sure whether our approach is appropriate (it throws NAs/0s)? Many thanks for any hints/advice/resources.
library(lme4)
library(simr)
n1 <- 62
n2 <- 53
n_total <- n1 + n2
n_workshops <- 8
simulate_data <- function(effect_size, n_total, n_workshops, n1, n2) {
workshop_ID <- factor(rep(1:n_workshops, length.out = n_total))
group <- factor(rep(c("Group1", "Group2"), times = c(n1, n2)))
DV_T0 <- rnorm(n_total)
DV_T1 <- DV_T0 + effect_size * (group == "Group2") + rnorm(n_total, 0, 1)
data.frame(DV_T1 = DV_T1,
DV_T0 = DV_T0,
group = group,
workshop_ID = workshop_ID)
}
effect_sizes <- seq(0.1, 1.0, by = 0.1)
power_analysis <- function(effect_size) {
simulated_data <- simulate_data(effect_size, n_total, n_workshops, n1, n2)
model <- lmer(DV_T1 ~ DV_T0 + group + (1 | workshop_ID), data = simulated_data)
if (isSingular(model)) {
return(NA)
}
power_simulation <- powerSim(model, nsim = 1000, test = fixed("groupGroup2"))
return(summary(power_simulation)$mean)
value
}
powers <- sapply(effect_sizes, power_analysis)
valid_indices <- !is.na(powers)
plot(effect_sizes[valid_indices], powers[valid_indices], type = "b", xlab = "Effect Size", ylab = "Power",
main = "Sensitivity Analysis")
abline(h = 0.80, col = "red", lty = 2)
power