Closed raoulvanoosten closed 10 months ago
pitakakariki
commented
10 months ago powerCurve works with the order of the data in the data frame (edit: not entirely true, see below). In your example, the binary variable bibukom is "binnen" for the first 40 observations and "buiten" for the next 40.
That means when powerCurve gets down to 40 or fewer rows, it can no longer fit the model. Should work if you alternate bibukom the way you've alternated campagne and fietstype though.
raoulvanoosten
commented
10 months ago Thanks for your reply. Good to know that the order is important. That can be changed so it will work. It does indeed work with the code below. Important is that the letters of block have to alternative (A, C, E, etc for "binnen", and B, D, F for "buiten").
# for-loop
df_sim <- c()
effect_campagne <- .05
for(i in 1:10){
size_nbikes_binnen_voor <- abs(round(rnorm(1, 42.2, 24.3))) # determine length
helm_nbikes_binnen_voor <- abs(rnorm(1, 0.0166, 0.0188)) # create binom
size_nbikes_binnen_na <- abs(round(rnorm(1, 42.2, 24.3)))
helm_nbikes_binnen_na <- abs(rnorm(1, 0.0166 + effect_campagne, 0.0188))
size_ebikes_binnen_voor <- abs(round(rnorm(1, 56.9, 48.8)))
helm_ebikes_binnen_voor <- abs(rnorm(1, 0.0512, 0.0265))
size_ebikes_binnen_na <- abs(round(rnorm(1, 56.9, 48.8)))
helm_ebikes_binnen_na <- abs(rnorm(1, 0.0512 + effect_campagne, 0.0265))
group_binnen <- LETTERS[(i*2)-1]
size_nbikes_buiten_voor <- abs(round(rnorm(1, 10.4, 6.72))) # determine length
helm_nbikes_buiten_voor <- abs(rnorm(1, 0.0510, 0.0633)) # create binom
size_nbikes_buiten_na <- abs(round(rnorm(1, 10.4, 6.72)))
helm_nbikes_buiten_na <- abs(rnorm(1, 0.0510 + effect_campagne, 0.0633))
size_ebikes_buiten_voor <- abs(round(rnorm(1, 19.8, 13.4)))
helm_ebikes_buiten_voor <- abs(rnorm(1, 0.0968, 0.0569))
size_ebikes_buiten_na <- abs(round(rnorm(1, 19.8, 13.4)))
helm_ebikes_buiten_na <- abs(rnorm(1, 0.0968 + effect_campagne, 0.0569))
group_buiten <- LETTERS[i * 2]
df_sim <- rbind(df_sim,
data.frame(block = group_binnen,
campagne = "voor",
fietstype = "normale fiets",
bibukom = "binnen",
aantal_fietsers = size_nbikes_binnen_voor,
proportie_helm = helm_nbikes_binnen_voor))
df_sim <- rbind(df_sim,
data.frame(block = group_binnen,
campagne = "na",
fietstype = "normale fiets",
bibukom = "binnen",
aantal_fietsers = size_nbikes_binnen_na,
proportie_helm = helm_nbikes_binnen_na))
df_sim <- rbind(df_sim,
data.frame(block = group_binnen,
campagne = "voor",
fietstype = "elektrische fiets",
bibukom = "binnen",
aantal_fietsers = size_ebikes_binnen_voor,
proportie_helm = helm_ebikes_binnen_voor))
df_sim <- rbind(df_sim,
data.frame(block = group_binnen,
campagne = "na",
fietstype = "elektrische fiets",
bibukom = "binnen",
aantal_fietsers = size_ebikes_binnen_na,
proportie_helm = helm_ebikes_binnen_na))
df_sim <- rbind(df_sim,
data.frame(block = group_buiten,
campagne = "voor",
fietstype = "normale fiets",
bibukom = "buiten",
aantal_fietsers = size_nbikes_buiten_voor,
proportie_helm = helm_nbikes_buiten_voor))
df_sim <- rbind(df_sim,
data.frame(block = group_buiten,
campagne = "na",
fietstype = "normale fiets",
bibukom = "buiten",
aantal_fietsers = size_nbikes_buiten_na,
proportie_helm = helm_nbikes_buiten_na))
df_sim <- rbind(df_sim,
data.frame(block = group_buiten,
campagne = "voor",
fietstype = "elektrische fiets",
bibukom = "buiten",
aantal_fietsers = size_ebikes_buiten_voor,
proportie_helm = helm_ebikes_buiten_voor))
df_sim <- rbind(df_sim,
data.frame(block = group_buiten,
campagne = "na",
fietstype = "elektrische fiets",
bibukom = "buiten",
aantal_fietsers = size_ebikes_buiten_na,
proportie_helm = helm_ebikes_buiten_na))
}
df_sim$block <- factor(df_sim$block)
df_sim$campagne <- factor(df_sim$campagne)
df_sim$fietstype <- factor(df_sim$fietstype)
df_sim$bibukom <- factor(df_sim$bibukom)
# modelleren
model_campagne <- lme4::glmer(proportie_helm ~ campagne + fietstype + bibukom + (1 | block), data = df_sim,
family = binomial, weight = aantal_fietsers)
simr::powerCurve(model_campagne, test = simr::fixed("campagne", method = "chisq"), along = "block", nsim = 10)
I have data of the proportion of cylists wearing helmets, and simulated new data from these: one set before a campaign with the same settings, and one set after the campaign with .03 more cyclists wearing helmets. I wanted to use
simr::powerSim()to calculate how many groups (blocks/locations) I need to sample for 90% power. I did that by running a generalized linear model (glmer). However, the simulations yield 0 power for 1 to 9 blocks, and 100% power for 10 blocks. Here is the output:This does not seem right to me. What is going on? I included the code below, which can be run with packages simr, lme4 and optionally emmeans.
The parameters in the data:
These parameters were relevant in the model of the actual data (high AIC)
I tried the
simr::powerSim()function with manual editing of the number of groups for 100 simulations each. This works. I think I should run/randomize the for-loop several times and then run many simulations.