You have been asked to design an (unmatched) cluster-randomised trial to measure the impact of immunisation with a newly developed pneumococcal vaccine on all-cause child mortality in a rural African population. Villages are to be randomly allocated to the intervention or control arms. In the intervention arm, all infants will receive the pneumococcal vaccine in a 3 dose regimen before the age of 6 months, and will then be followed up for two years (to age 30 months) to record mortality. In the control arm, infants will receive a placebo vaccine and will be followed up in the same way.
Current child mortality in the age-range 6-29 months is estimated to be 35 per 1000 child-years, and you estimate that this is likely to vary between 25 and 45 per 1000 child-years in most villages. You wish to have 80% power of detecting a 20% reduction in all-cause mortality due to the pneumococcal vaccine. Each year, around 100 infants reach the age of six months in each village, and you plan to recruit infants to the trial for three years.
Estimate roughly how many villages you will need to recruit to each arm of the trial. What is the design effect?
library(haven)
# create table
library(gtsummary)
# data wrangling and other functions
library(tidyverse)
# if you want to have a output similar to stata
library(tidylog)
#cox
library(survival)
# Limit significant digits to 3, remove scientific notation
options(digits = 3, scipen = 9)t <- 3
m <- 100
g <- 2
lc <- 0.035
le <- 0.035*0.8
alpha <- 0.05
power <- 0.8
CV <- 0.14
T <- t * m * g
IFt <- ((lc +le)/ T + (CV^2) * (lc^2 + le^2))
f <- (qnorm(1 - alpha/2) + qnorm(power))^2
sample <- (IFt * f)/(lc - le)^2 +1This question uses data from the Gambian National Impregnated Bednet Programme, datasets GAMINDIV.DTA and GAMVILL.DTA. The data are from a cross-sectional survey of children aged 1-4 years conducted at the end of the transmission season in a sample of villages from each treatment arm. See course manual for further details.
dt_cluster <- read_dta("datasets/GAMVILL.DTA") |>
mutate(across(where(is.labelled),as_factor))
dt_ind <- read_dta("datasets/GAMINDIV.DTA") |>
mutate(across(where(is.labelled),as_factor))Using GAMVILL.DTA, look at histograms of village parasite prevalence by treatment arm. Is there any apparent difference?
Variables that we will use:
**Outcome*: - in individual dataset: para (malaria status: 1 = positive, 2 = negative) - in summary dayaset: parapos (number of children positive), paraneg (number of children negative), rate (village parasite prevalence: parapos/(parapos+paraneg)) **Exposure*: group (treatment arm: 1 = bednets, 2 = control) **Cluster*: vid (village)
dt_cluster |>
mutate(group = if_else(group==1,"Bednets","Control")) |>
ggplot() +
geom_histogram(aes(rate,fill=group), binwidth = 0.1, color="white")+
scale_y_continuous(breaks = seq(0:10))
Estimate the overall parasite prevalence for each group: (i) by summing over all children (GAMINDIV.DTA); and (ii) by averaging village prevalences (GAMVILL.DTA). Is there any difference between these two estimates?
dt_ind |> count(para) |>
mutate(perc = prop.table(n))| para | n | perc |
|---|---|---|
| 1 | 385 | 0.33 |
| 2 | 781 | 0.67 |
dt_cluster |> summarise(prev = mean(rate))| prev |
|---|
| 0.335 |
Test the hypothesis that net impregnation reduces the prevalence of malaria parasites using the standard chi-squared test, ignoring the randomisation by village (use GAMINDIV.DTA).
tab1 <- table(dt_ind$para, dt_ind$group)
chisq.test(tab1)
Pearson's Chi-squared test with Yates' continuity correction
data: tab1
X-squared = 6, df = 1, p-value = 0.01Test the same null hypothesis, but allowing for randomisation by village, by using (i) a t-test on village prevalences; (ii) a Wilcoxon rank sum test on village prevalences.
t.test(rate~group, data=dt_cluster,
var.equal = T)
Two Sample t-test
data: rate by group
t = -2, df = 32, p-value = 0.06
alternative hypothesis: true difference in means between group 1 and group 2 is not equal to 0
95 percent confidence interval:
-0.23729 0.00673
sample estimates:
mean in group 1 mean in group 2
0.277 0.392 wilcox.test(rate~group, data=dt_cluster,
correct=F)
Wilcoxon rank sum test
data: rate by group
W = 96, p-value = 0.09
alternative hypothesis: true location shift is not equal to 0The default value for t.test in R is considering unequal variances (while in stata is considering equal variances), that’s why we need to include “var.equal=T” to get the same answers. Similarly, for the Wilcoxon, by default it applies a continuity correction in the normal approximation for the p-value (not by default in stata)
Compare the results you obtain from (d) with your unadjusted result from (c). How does the adjustment for clustering change your interpretation of the effect of bednet impregnation on parasite status?
The standard \(\chi^2\) test shows a highly significant effect of treatment, and is invalid! The other tests, which allow for the clustered design, give larger P-values indicating only weak evidence that the treatment has an effect on parasite status.
#get the risks
# dt_ind |> group_by(group) |>
# count(para) |>
# mutate(risk=n/sum(n)) |>
# filter(para==1) |>
# ungroup() |>
# mutate(rr = risk[group=="1"]/risk[group=="2"])
# Risk in each arm
r0 <- 0.367
r1 <- 0.298
# Risk ratio
RR <- r1/r0
logRR <- log(RR)
# Standard deviations of rates in each group (by cluster )
#dt_cluster |> group_by(group) |> summarise(sd(rate)) - get the SD
s1 <- 0.15
s0 <- 0.196
# Numbers of clusters
c1 <- 17
c0 <- 17
# Variances of risks
var_r1 <- (s1^2)/c1
var_r0 <- (s0^2)/c0
# Variance of the log RR
var_logRR <- var_r1/(r1^2) + var_r0/(r0^2)
# Confidence intervals
lower_CI <- exp(logRR - 1.96 * sqrt(var_logRR))
upper_CI <- exp(logRR + 1.96 * sqrt(var_logRR))
data.frame(rr= RR, lcl = lower_CI, ucl = upper_CI)| rr | lcl | ucl |
|---|---|---|
| 0.812 | 0.573 | 1.15 |
library(geepack)
library(lme4)
dt_ind <- dt_ind |>
mutate(para_bin = (if_else(para==1,1,0)),
vid_fac = as.factor(vid),
group_bin = as.factor(if_else(group==2,0,1)))
dt_ind_order <- dt_ind |>
arrange(vid_fac) # geeglm (geepack) assumes that the observations from the same cluster appears contiguous
mod_unorder <- geeglm(para_bin~group_bin,
id = vid_fac,
data=dt_ind,
corstr = "exchangeable",
family = binomial(link="logit"))
mod_order <- geeglm(para_bin~group_bin,
id = vid_fac,
data=(dt_ind_order),
corstr = "exchangeable",
family = binomial(link="logit"))
#broom::tidy(mod_unorder, exponentiate=T) #check that it gives the wrong answer
broom::tidy(mod_order, exponentiate = T) |> gt::gt()| term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|
| (Intercept) | 0.635 | 0.196 | 5.37 | 0.0205 |
| group_bin1 | 0.610 | 0.264 | 3.50 | 0.0614 |
# geeM doesnt require ordering
# mod <- geeM::geem(para_bin~group_bin,
# id = vid_fac,
# data=as.data.frame(dt_ind),
# corstr = "exchangeable",
# family = binomial(link="logit"))
# summary(mod)
mod <- glmer(para_bin~group_bin + (1|vid),
data=dt_ind,
family = "binomial")
broom.mixed::tidy(mod, exponentiate=T, effects="fixed") |> gt::gt()| effect | term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|---|
| fixed | (Intercept) | 0.593 | 0.124 | -2.50 | 0.0124 |
| fixed | group_bin1 | 0.584 | 0.174 | -1.81 | 0.0705 |
Dataset MZTRIAL.DTA contains individual-level data from the Mwanza trial of the impact of improved STD treatment services on HIV incidence. There were 12 communities in 6 matched pairs, and a random cohort of around 1000 adults aged 15-54 years was followed up in each community. The dataset contains data on HIV seroconversion among those individuals who were seronegative at baseline and who were successfully followed up after two years. We are going to analyse the results of the trial with and without adjustment for confounding variables (see Section 6).
dt_mzt <- read_dta("datasets/MZTRIAL.DTA") |>
mutate(across(where(is.labelled),as_factor))
dt_mzt <- dt_mzt |> mutate(arm = factor(arm,
levels = c(0, 1),
labels = c("control", "intervention")),
comp = factor(comp),
pair = factor(pair),
agegp = factor(agegp,
levels = 1:4,
labels = c("15-24", "25-34", "35-44", "45-54")),
sex = factor(sex,
levels = c(1, 2),
labels = c("male", "female")))model = glm(hiv ~ agegp + sex + pair + hivbase) predict(model, type="response")
mod_hiv_prob <- glm(hiv ~ agegp + sex + pair + hivbase,
family = binomial,
data = dt_mzt)
dt_mzt$fitted_hiv_prob <- predict(mod_hiv_prob, type = "response")dt_mzt |>
ggplot(aes(fitted_hiv_prob))+
geom_histogram(aes(fill=arm))
The predict takes a model and give the prediction using a new dataset. In the case of no “newdata” being declared, it will use the dataset used for fitting the model. The type="response" is related to give the prediction in the response scale, i.e., in the probability scale
dt_mzt|> group_by(pair, comp, arm) %>%
summarise("observed" = sum(hiv),
"at_risk" = n(),
"expected" = sum(fitted_hiv_prob))| pair | comp | arm | observed | at_risk | expected |
|---|---|---|---|---|---|
| 1 | 1 | intervention | 5 | 568 | 7.16 |
| 1 | 2 | control | 10 | 702 | 7.84 |
| 2 | 3 | intervention | 4 | 766 | 5.34 |
| 2 | 4 | control | 7 | 833 | 5.66 |
| 3 | 5 | intervention | 17 | 650 | 16.66 |
| 3 | 6 | control | 20 | 630 | 20.34 |
| 4 | 7 | intervention | 13 | 734 | 18.77 |
| 4 | 8 | control | 23 | 760 | 17.23 |
| 5 | 9 | intervention | 4 | 732 | 6.89 |
| 5 | 10 | control | 12 | 782 | 9.11 |
| 6 | 11 | intervention | 5 | 699 | 6.46 |
| 6 | 12 | control | 10 | 693 | 8.54 |
dd_unaj <- dt_mzt|> group_by(pair, comp, arm) %>%
summarise("observed" = sum(hiv),
"at_risk" = n(),
"expected" = sum(fitted_hiv_prob)) |>
group_by(pair) |>
mutate(rr = observed[arm=="intervention"]/at_risk[arm=="intervention"]/(observed[arm=="control"]/at_risk[arm=="control"])) |>
filter(arm=="intervention") |>
select(pair, rr)
dd_adj <- dt_mzt|> group_by(pair, comp, arm) %>%
summarise("observed" = sum(hiv),
"at_risk" = n(),
"expected" = sum(fitted_hiv_prob)) |>
group_by(pair) |>
mutate(adjrr = observed[arm=="intervention"]/expected[arm=="intervention"]/(observed[arm=="control"]/expected[arm=="control"])) |>
filter(arm=="intervention") |>
select(pair, adjrr)
dd_unaj |>
left_join(dd_adj, by = "pair")| pair | rr | adjrr |
|---|---|---|
| 1 | 0.618 | 0.548 |
| 2 | 0.621 | 0.606 |
| 3 | 0.824 | 1.038 |
| 4 | 0.585 | 0.519 |
| 5 | 0.356 | 0.441 |
| 6 | 0.496 | 0.661 |
Calculated the paired t-tests in R
Unadjusted
dd_unadj_test <- dt_mzt|> group_by(pair, comp, arm) %>%
summarise("observed" = sum(hiv),
"at_risk" = n(),
"expected" = sum(fitted_hiv_prob)) |>
group_by(pair) |>
mutate(logrr = log(observed/at_risk))
t_test_values<- t.test(dd_unadj_test[dd_unadj_test$arm=="intervention",]$logrr,dd_unadj_test[dd_unadj_test$arm=="control",]$logrr,paired=T)
t_test_values
Paired t-test
data: dd_unadj_test[dd_unadj_test$arm == "intervention", ]$logrr and dd_unadj_test[dd_unadj_test$arm == "control", ]$logrr
t = -5, df = 5, p-value = 0.004
alternative hypothesis: true mean difference is not equal to 0
95 percent confidence interval:
-0.864 -0.277
sample estimates:
mean difference
-0.57 data.frame(estimate = exp(t_test_values$estimate), #need to back transform
lcl = exp(t_test_values$conf.int[1]),
ucl = exp(t_test_values$conf.int[2]))| estimate | lcl | ucl | |
|---|---|---|---|
| mean difference | 0.565 | 0.422 | 0.758 |
Adjusted
dd_adj_test <- dt_mzt|> group_by(pair, comp, arm) %>%
summarise("observed" = sum(hiv),
"at_risk" = n(),
"expected" = sum(fitted_hiv_prob)) |>
group_by(pair) |>
mutate(logadjrr = log(observed/expected))
t_test_values<- t.test(dd_adj_test[dd_adj_test$arm=="intervention",]$logadjrr,dd_adj_test[dd_adj_test$arm=="control",]$logadjrr,paired=T)
t_test_values
Paired t-test
data: dd_adj_test[dd_adj_test$arm == "intervention", ]$logadjrr and dd_adj_test[dd_adj_test$arm == "control", ]$logadjrr
t = -4, df = 5, p-value = 0.009
alternative hypothesis: true mean difference is not equal to 0
95 percent confidence interval:
-0.801 -0.184
sample estimates:
mean difference
-0.492 data.frame(estimate = exp(t_test_values$estimate),
lcl = exp(t_test_values$conf.int[1]),
ucl = exp(t_test_values$conf.int[2]))| estimate | lcl | ucl | |
|---|---|---|---|
| mean difference | 0.611 | 0.449 | 0.832 |