I’ve been experimenting with the multinma package for multi-level network meta-regression (ML-NMR), and I’ve run into something that’s left me unsure how to interpret LOO-CV results in this setting.
In the multinma package it looks like each individual patient from IPD contributes one observation to the likelihood, while each arm from an aggregate (AgD) study contributes a single observation. This means that if we use LOO on the combined IPD+AgD network, the number of “data points” driving the estimate is the sum of IPD patients and aggregate arms.
This raises two questions I’m not sure how to answer
- Are we really getting a mix of estimates of leave-one-patient-out error and leave-one-study-out error? I tried to explore this a little bit below by comparing an aggregate vs individual model on the aggregate data studies only since when not adjusting it just involves expanding the binomial counts and I do indeed get very different estimates from loo.
- Could this mixed LOO drive different model selection then having everything be leave one study out? My thinking is that for if example If LOO is dominated by the many IPD observations, could this cause model selection to favor a model that predicts well within the IPD study but poorly in the AgD studies?
- I notice that LOO has all bad pareto k estimates when evaluated in aggregate data only. I assume that means we’d have all bad or many bad pareto k estimates if trying to estimate all as leave-one-study out? When expanding the aggregate data to individual patients I still get 2 bad and 1 very bad pareto k but all others are now “good”.
Minimal example
library(multinma)
library(loo)
library(dplyr)
# Load psoriasis data (includes both IPD and AgD)
pso_ipd <- filter(plaque_psoriasis_ipd,
studyc %in% c("UNCOVER-1", "UNCOVER-2", "UNCOVER-3"))
pso_agd <- filter(plaque_psoriasis_agd,
studyc == "FIXTURE")
pso_ipd <- pso_ipd %>%
mutate(# Variable transformations
bsa = bsa / 100,
prevsys = as.numeric(prevsys),
psa = as.numeric(psa),
weight = weight / 10,
durnpso = durnpso / 10,
# Treatment classes
trtclass = case_when(trtn == 1 ~ "Placebo",
trtn %in% c(2, 3, 5, 6) ~ "IL blocker",
trtn == 4 ~ "TNFa blocker"),
# Check complete cases for covariates of interest
complete = complete.cases(durnpso, prevsys, bsa, weight, psa)
) |>
dplyr::filter(complete)
pso_agd <- pso_agd %>%
mutate(
# Variable transformations
bsa_mean = bsa_mean / 100,
bsa_sd = bsa_sd / 100,
prevsys = prevsys / 100,
psa = psa / 100,
weight_mean = weight_mean / 10,
weight_sd = weight_sd / 10,
durnpso_mean = durnpso_mean / 10,
durnpso_sd = durnpso_sd / 10,
# Treatment classes
trtclass = case_when(trtn == 1 ~ "Placebo",
trtn %in% c(2, 3, 5, 6) ~ "IL blocker",
trtn == 4 ~ "TNFa blocker")
)
expand_pasi75_to_ipd <- function(agd) {
do.call(rbind, lapply(seq_len(nrow(agd)), function(i) {
row <- agd[i, ]
# responders: pasi75_r rows of 1
responders <- data.frame(
studyc = row$studyc,
trtc = row$trtc,
pasi75 = rep(1, row$pasi75_r)
)
# non-responders: (pasi75_n - pasi75_r) rows of 0
nonresponders <- data.frame(
studyc = row$studyc,
trtc = row$trtc,
pasi75 = rep(0, row$pasi75_n - row$pasi75_r)
)
rbind(responders, nonresponders)
}))
}
ipd_like <- expand_pasi75_to_ipd(pso_agd)
pso_agd_as_ipd_net <- set_ipd(ipd_like,
study = studyc,
trt = trtc,
r = pasi75)
pso_net <- combine_network(
set_ipd(pso_ipd,
study = studyc,
trt = trtc,
r = pasi75,
trt_class = trtclass),
set_agd_arm(pso_agd,
study = studyc,
trt = trtc,
r = pasi75_r,
n = pasi75_n,
trt_class = trtclass)
)
pso_net <- add_integration(pso_net,
durnpso = distr(qgamma, mean = durnpso_mean, sd = durnpso_sd),
prevsys = distr(qbern, prob = prevsys),
bsa = distr(qlogitnorm, mean = bsa_mean, sd = bsa_sd),
weight = distr(qgamma, mean = weight_mean, sd = weight_sd),
psa = distr(qbern, prob = psa),
n_int = 64
)
pso_net_agd <- set_agd_arm(pso_agd,
study = studyc,
trt = trtc,
r = pasi75_r,
n = pasi75_n,
trt_class = trtclass)
pso_fit_FE <- nma(pso_net,
trt_effects = "fixed",
link = "probit",
likelihood = "bernoulli2",
class_interactions = "common",
prior_intercept = normal(scale = 10),
prior_trt = normal(scale = 10),
prior_reg = normal(scale = 10),
init_r = 0.1,
QR = TRUE)
pso_fit_FE_agd <- nma(pso_net_agd,
trt_effects = "fixed",
link = "probit",
likelihood = "binomial2",
prior_intercept = normal(scale = 10),
prior_trt = normal(scale = 10),
prior_reg = normal(scale = 10),
init_r = 0.1,
QR = TRUE)
pso_fit_FE_agd_as_ipd <- nma(pso_agd_as_ipd_net,
trt_effects = "fixed",
link = "probit",
likelihood = "bernoulli2",
prior_intercept = normal(scale = 10),
prior_trt = normal(scale = 10),
prior_reg = normal(scale = 10),
init_r = 0.1,
QR = TRUE)
# Compute LOO
loo_combined <- loo(pso_fit_FE)
loo_agd_only <- loo(pso_fit_FE_agd)
loo_agd_as_ipd <- loo(pso_fit_FE_agd_as_ipd)
print(loo_combined) # ll matrix is 4000 x 3858 consistent with one obs per IPD patient + one per agd arm
print(loo_agd_only) # ll matrix is 4000 x 4 consistent with one obs per study arm
print(loo_agd_as_ipd) # ll matrix is 4000 x 1297
nrow(pso_ipd) + nrow(pso_agd)
nrow(pso_agd)