I am estimating a mixed logit model and I have run the code for the power analysis (without interaction terms) and I would like to know how to run a similar power analysis once I have interaction terms. For ML1 (without interaction terms) I have used the code below based on Bekker-Grob et al. (2015) and then extended the code for ML2 when I have interactions (see second part of the code) but I find that I need less observations for some interaction terms which does not make sense so there must be something wrong in my code for ML2. What approach should I take when running a power analysis for a mixed logit with interaction terms? .
Thank you so much for the help!
Code: Select all
## When using this code, please cite this article: de Bekker-Grob, E.W., Donkers, B., Jonker, M.F. et al. Patient (2015) 8: 373. https://doi.org/10.1007/s40271-015-0118-z##
############################################################################################################################################################################
############################################################################################################################################################################
library(tidyr)
library(dplyr)
## ================================================================
## 1. Reshape to long format
## ================================================================
df_long <- database %>%
pivot_longer(
cols = starts_with("alt"),
names_to = c("Alternative", ".value"),
names_pattern = "alt(\\d+)_(.*)"
)
df_long <- df_long %>%
mutate(constant = ifelse(Alternative == 3, 1, 0))
df_long <- df_long %>%
mutate(supportindiv = ifelse(support == 2, 1, 0))
df_long <- df_long %>%
mutate(supportcoop = ifelse(support == 3, 1, 0))
df_long <- df_long %>%
mutate(plantindp = ifelse(plantation == 2, 1, 0))
df_long <- df_long %>%
mutate(plantcoop = ifelse(plantation == 3, 1, 0))
df_long <- df_long %>%
mutate(sprayindgp = ifelse(spraying == 2, 1, 0))
df_long <- df_long %>%
mutate(spraycoop = ifelse(spraying == 3, 1, 0))
df_long <- df_long %>%
mutate(techcoll = ifelse(suivi == 2, 1, 0))
df_long <- df_long[ , -c(4:5, 7:8)]
df_long <- df_long %>%
select(
1:3,
5:9,
4,
everything()
)
## ================================================================
## 2. Base (fixed) attribute design matrix, from the Ngene design
## ================================================================
design_base <- as.matrix(df_long[, -c(1:3)])
colnames(design_base) <- c("ASC_SQ", "supportindiv", "supportcoop",
"plantindp", "plantcoop", "investcost",
"sprayindgp", "spraycoop", "techcoll")
nalts <- 3
nchoices <- 24
test_alpha <- 0.05
z_one_minus_alpha <- qnorm(1 - test_alpha)
test_beta <- 0.20
z_one_minus_beta <- qnorm(1 - test_beta)
## ================================================================
## 3. ML1 power analysis
## ================================================================
parameters_ML1 <- c(-11.925, 0.741, 0.932, -0.017, -0.172, -1.418, -0.953, -1.354, -0.448)
names(parameters_ML1) <- colnames(design_base)
ncoefficients_ML1 <- length(parameters_ML1)
info_mat_ML1 <- matrix(0, ncoefficients_ML1, ncoefficients_ML1)
exputilities_ML1 <- exp(design_base %*% parameters_ML1)
for (k_set in 1:nchoices) {
alternatives <- ((k_set - 1) * nalts + 1):(k_set * nalts)
p_set <- exputilities_ML1[alternatives] / sum(exputilities_ML1[alternatives])
p_diag <- diag(p_set)
middle_term <- p_diag - p_set %o% p_set
full_term <- t(design_base[alternatives, ]) %*% middle_term %*% design_base[alternatives, ]
info_mat_ML1 <- info_mat_ML1 + full_term
}
sigma_beta_ML1 <- solve(info_mat_ML1, diag(ncoefficients_ML1))
N_ML1 <- ((z_one_minus_beta + z_one_minus_alpha) * sqrt(diag(sigma_beta_ML1)) / abs(parameters_ML1))^2
names(N_ML1) <- names(parameters_ML1)
cat("=== Required sample size per coefficient - ML1 ===\n")
print(N_ML1)
## ================================================================
## 4. ML2 power analysis, including interaction terms
## ================================================================
set.seed(123)
parameters_ML2 <- c(
ASC_SQ = -8.019,
ASC_SQ_Age45 = -0.857,
ASC_SQ_Female = -1.178,
ASC_SQ_HigherEdu = 0.431,
ASC_SQ_FarmsizeNSZ = 0.005,
ASC_SQ_KnowledgeFRG = -3.783,
ASC_SQ_ExternalSpray = 1.778,
ASC_SQ_NoEquip = 0.012,
ASC_SQ_Renew = -0.129,
ASC_SQ_NoSellWine = 0.633,
ASC_SQ_Coop = -0.881,
ASC_SQ_Organic = -0.293,
ASC_SQ_Mildew = 0.973,
ASC_SQ_Turnover = 0.0007,
ASC_SQ_Tradition = 0.396,
ASC_SQ_BeliefsKnowledge = -1.380,
ASC_SQ_BeliefsEco = 0.134,
supportindiv = 0.689,
supportcoop = 1.963,
plantindp = -1.301,
plantcoop = -0.757,
investcost = -1.328,
sprayindgp = -2.273,
spraycoop = -2.415,
techcoll = -0.363,
plantindp_BeliefsEco = 1.722,
plantcoop_BeliefsEco = 0.857,
sprayindgp_BeliefsEco = 1.673,
spraycoop_BeliefsEco = 1.311,
supportcoop_BeliefsKnowledge = -1.109,
techcoll_BeliefsKnowledge = -0.070
)
ncoefficients_ML2 <- length(parameters_ML2) # 31
n_draws <- 2000
covariate_draws <- data.frame(
Age45 = rbinom(n_draws, 1, 0.48),
Female = rbinom(n_draws, 1, 0.167),
HigherEdu = rbinom(n_draws, 1, 0.76),
FarmsizeNSZ = rnorm(n_draws, mean = 2, sd = X),
KnowledgeFRG = rbinom(n_draws, 1, 0.7),
ExternalSpray = rbinom(n_draws, 1, X),
NoEquip = rbinom(n_draws, 1, 0.29),
Renew = rbinom(n_draws, 1, 0.43),
NoSellWine = rbinom(n_draws, 1, 0.146),
Coop = rbinom(n_draws, 1, 0.85),
Organic = rbinom(n_draws, 1, 0.17),
Mildew = rbinom(n_draws, 1, 0.55),
Turnover = rnorm(n_draws, mean = 506,993, sd = x),
Tradition = rbinom(n_draws, 1, 0.5),
BeliefsKnowledge = rnorm(n_draws, mean = 0, sd = 1),
BeliefsEco = rnorm(n_draws, mean = 0, sd = 1)
)
## -----------------------------------------------------------------------
ASC <- design_base[, "ASC_SQ"]
plantindp_d <- design_base[, "plantindp"]
plantcoop_d <- design_base[, "plantcoop"]
sprayindgp_d <- design_base[, "sprayindgp"]
spraycoop_d <- design_base[, "spraycoop"]
supportcoop_d <- design_base[, "supportcoop"]
techcoll_d <- design_base[, "techcoll"]
info_mat_ML2 <- matrix(0, ncoefficients_ML2, ncoefficients_ML2)
for (r in 1:n_draws) {
cov_r <- covariate_draws[r, ]
design_r <- cbind(
ASC,
ASC * cov_r$Age45,
ASC * cov_r$Female,
ASC * cov_r$HigherEdu,
ASC * cov_r$FarmsizeNSZ,
ASC * cov_r$KnowledgeFRG,
ASC * cov_r$ExternalSpray,
ASC * cov_r$NoEquip,
ASC * cov_r$Renew,
ASC * cov_r$NoSellWine,
ASC * cov_r$Coop,
ASC * cov_r$Organic,
ASC * cov_r$Mildew,
ASC * cov_r$Turnover,
ASC * cov_r$Tradition,
ASC * cov_r$BeliefsKnowledge,
ASC * cov_r$BeliefsEco,
design_base[, "supportindiv"],
design_base[, "supportcoop"],
design_base[, "plantindp"],
design_base[, "plantcoop"],
design_base[, "investcost"],
design_base[, "sprayindgp"],
design_base[, "spraycoop"],
design_base[, "techcoll"],
plantindp_d * cov_r$BeliefsEco,
plantcoop_d * cov_r$BeliefsEco,
sprayindgp_d * cov_r$BeliefsEco,
spraycoop_d * cov_r$BeliefsEco,
supportcoop_d * cov_r$BeliefsKnowledge,
techcoll_d * cov_r$BeliefsKnowledge
)
exputilities_r <- exp(design_r %*% parameters_ML2)
for (k_set in 1:nchoices) {
alternatives <- ((k_set - 1) * nalts + 1):(k_set * nalts)
p_set <- exputilities_r[alternatives] / sum(exputilities_r[alternatives])
p_diag <- diag(p_set)
middle_term <- p_diag - p_set %o% p_set
full_term <- t(design_r[alternatives, ]) %*% middle_term %*% design_r[alternatives, ]
info_mat_ML2 <- info_mat_ML2 + full_term
}
}
## population-average information matrix per respondent
info_mat_ML2 <- info_mat_ML2 / n_draws
sigma_beta_ML2 <- solve(info_mat_ML2, diag(ncoefficients_ML2))
N_ML2 <- ((z_one_minus_beta + z_one_minus_alpha) * sqrt(diag(sigma_beta_ML2)) / abs(parameters_ML2))^2
names(N_ML2) <- names(parameters_ML2)
cat("\n=== Required sample size per coefficient - ML2 ===\n")
print(N_ML2)
cat("\nOverall required sample size (max across coefficients):\n")
cat("ML1:", ceiling(max(N_ML1)), "\n")
cat("ML2:", ceiling(max(N_ML2)), "\n")