Overview
In a previous blogpost, Comprehending complex designs: Cluster randomized trials, I walked through the nuances and challenges of cluster randomized trials (CRTs). Cluster randomized trials randomize groups of individuals, such as families or clinics, rather than individuals themselves. Cluster randomized trials are used for a variety of reasons, including evaluating the spread of infectious disease within a household or evaluating whether a new intervention is effective or feasible in real-world settings. Participants within the same cluster may share the same environment or care provider, for example, leading to correlated responses. If this intracluster correlation is not accounted for, variances will be underestimated and inference methods will not have the operating characteristics (i.e., type I error) we expect. Linear mixed models represent one approach for obtaining cluster-adjusted estimates, and their application was demonstrated using data from the SHARE cluster trial evaluating different sex ed curriculums (interventions) in schools (clusters).
Individually randomized group treatment trials (IRGTs) are closely related to CRTs, but can require slightly more complex analytic strategies. IRGT designs arise naturally when individuals do not initially belong to a group or cluster, but are individually randomized to receive a group-based intervention or receive treatment through a shared agent. As a result, individuals are independent at baseline, but intracluster correlation can increase with follow-up as individuals interact within their respective group or with their shared agent. IRGTs can be “fully-nested,” meaning that both the control and experimental conditions feature a group-based intervention, or “partially-nested,” meaning that the experimental condition is group-based while the control arm is not. A fully-nested IRGT may be used to compare structured group therapy versus group discussion for mental health outcomes, for example. If both arms feature groups and the same intracluster correlation, analysis of fully-nested IRGTs is practically identical to that of CRTs. In comparison, a partially-nested IRGT may be used to compare group therapy versus individual standard of care or a waitlist control, for example. Analysis of partially-nested IRGTs is more complex because intracluster correlation is only present in one arm, and methods must be adapted to handle heterogeneous covariance or correlation matrices. If fully-nested but arms do not share the same correlation, similar considerations are required.
To provide insight into data generating mechanisms and inference, this blog post demonstrates how to simulate normally distributed outcomes from (1) a two-arm cluster randomized trial and (2) a two-arm, partially-nested individually randomized group treatment trial. I only use base R for data generation, so these approaches can be widely implemented. Simulation of complex trial designs is helpful for sample size calculation and understanding operating characteristics of inference methods in different scenarios, such as small samples. Analysis of the simulated data proceeds using linear mixed models fit by the nlme library. Visualization uses ggplot2.
Cluster randomized trial
Consider a two-arm cluster randomized trial which randomly allocates clusters to a treatment arm and a control arm. There are 
treated clusters and 
control clusters, corresponding to a total of 
clusters randomized. Each cluster features the same number of participants, 
, corresponding to a total of 
participants. We will generate normally distributed responses according to a specified intracluster correlation and effect size.
Within each cluster, responses are generated according to 
where 
represents the expected response in the 
arm (i=1 if treatment cluster, i=0 if control) and 
represents the total variance, or the sum of the between-cluster variance 
and within-cluster variance 
. Here, we assume that arms have the same variance and intracluster correlation so 
.The intracluster correlation is the proportion of total variation explained by differences between clusters, 
. If we know the ICC and the total variance, we can recover the values of 
and 
Cohen’s effect size is defined as, 
, i.e., the standardized mean difference. Cohen proposed that 
corresponds to a “small effect,” 
a “moderate effect,” and 
a “large effect.” If we fix the treatment mean, 
, and the total variance , we must set the control mean 
to obtain desired effect size 
. We reference to the unstandardized difference in means as the average treatment effect, 
.
In the following code block, we first specify key parameters and calculate all others. More details on data structures and parameters can be found in my previous blog post on cluster randomized trials.
# STEP 0: Specify important design and response parameters ---------------------
## Seed for reproducibility
set.seed(1243)
## Helper function
r2 <- function(x) round(x,2) # function to clean up reported numbers
## Design
r = 0.5 # randomization ratio (here, r=0.5 implies 1:1 trt:ctrl)
C = 1000 # expected number of total clusters
C1 = r*C # expected number of treatment clusters
C0 = (1-r)*C # expected number of control clusters
n = 10 # expected number of participants per cluster
N = C*n # expected number of total participants
## Intracluster correlation and variance
icc = 0.1 # intracluster correlation = between / (between + within)
var = 16 # total variance = between + within cluster variance
var_a = icc * var # between-cluster variance (random intercept)
var_e = var - var_a # within-cluster variance (residual)
## Effects and means
ES = 0.5 # cohen's effect size (0.2=small, 0.5=moderate, 0.8=large)
mu1 = 25 # mean in treatment arm (fixed)
mu0 = mu1 - ES*sqrt(var) # mean in control arm (mu1 - ES*sd)
ATE = mu1 - mu0 # average treatment effect
## Summarize parameter values
param_text <- paste0("You are generating a cluster randomized trial with ",
C, " clusters. ",
r*100,"% of clusters are expected in the treatment arm and ",
(1-r)*100,"% in the control arm. ",
n, " participants are expected in each cluster, for a total of ",
N, " participants. ",
"The intracluster correlation is ", icc,
" and the total variance is ", var,
", corresponding to a random intercept variance of ", var_a,
" and residual variance ", var_e, ". The effect size is ", ES,
" and the treatment mean is ", mu1,
", corresponding to a control mean of ", r2(mu0),
". The average treatment effect (ATE) is ", r2(ATE),".")
print(param_text)
## [1] "You are generating a cluster randomized trial with 1000 clusters. 50% of clusters are expected in the treatment arm and 50% in the control arm. 10 participants are expected in each cluster, for a total of 10000 participants. The intracluster correlation is 0.1 and the total variance is 16, corresponding to a random intercept variance of 1.6 and residual variance 14.4. The effect size is 0.5 and the treatment mean is 25, corresponding to a control mean of 23. The average treatment effect (ATE) is 2."
Now that we have specified that we have participants in each of clusters, we need to assign a unique ID to each cluster (cid) and to each participant (pid) within a given cid.
# STEP 1: Initialize dataframe structure to store simulated cluster responses --
df <- data.frame(cid = rep(1:C, each=n), # cluster id assigned to each participant
pid = rep(1:n, times=C)) # participant id repeated within each cluster
Each participant has an associated residual, which is generated 
according to within-cluster variance 
. Because all residuals share the same distribution, we can simply generate and attach 
residuals to our existing dataframe.
# STEP 2: Generate residual for each participant, epsilon ~ N(0, var_e) --------
df$epsilon <- rnorm(n=N, mean=0, sd=sqrt(var_e))
Generating random intercepts for each cluster follows a similar process. We generate random intercepts according to 
where 
is the between-cluster variance. At the same time, we also randomly assign each cluster to treatment or control according to 
where 
is the desired proportion of treated clusters. As a result, the observed number of treatment and control clusters may not match our desired 
and 
, but should be very close.
# STEP 3: Generate cluster intercepts, alpha ~ N(0, var_a), and randomized arm -
df_a <- data.frame(cid = 1:C,
arm = rbinom(n=C, size=1, prob=r),
alpha = rnorm(n=C, mean=0, sd=sqrt(var_a)))
cluster_text <- paste0("C1=", sum(df_a$arm),
" clusters were randomized to the treated arm. ",
"C0=", C-sum(df_a$arm),
" clusters were randomized to the control arm.")
cluster_text
## [1] "C1=502 clusters were randomized to the treated arm. C0=498 clusters were randomized to the control arm."
Now, we need to attach the random cluster intercept to each participant within cluster cid.
# STEP 4: Attach random intercept to each participant within cid ---------------
df <- merge(df, df_a, by='cid')
We generate responses according to a random intercept model. That is, for the 
participant 
in the 
cluster 
in the arm 
,
![\[Y_{ijk} = \beta_0 + \beta_1 \mathbf{1}(i=1) + \alpha_{ij} + \epsilon_{ijk}\]](https://statisticelle.com/wp-content/ql-cache/quicklatex.com-a06c6aff440b7844710bf23a43324134_l3.png "Rendered by QuickLaTeX.com")
where 
equals 1 if a cluster belongs to the treatment and 0 otherwise. It follows that 
corresponds to the control mean, 
the average treatment effect, 
the random intercept for the 
cluster, and 
the residual for the 
participant.
We’ve done the hard work of simulating the random components of this equation, now we just need to add them all together! Voila, generated responses!
# STEP 5: Generate responses according to random intercept structure -----------
df$Y <- mu0 + ATE * df$arm + df$alpha + df$epsilon
## Preview dataframe of generated responses:
head(df)
## cid pid epsilon arm alpha Y ## 1 1 1 -1.4359086 0 -0.75642 20.80767 ## 2 1 2 -3.8790426 0 -0.75642 18.36454 ## 3 1 3 -2.6940894 0 -0.75642 19.54949 ## 4 1 4 3.1403289 0 -0.75642 25.38391 ## 5 1 5 -0.8084272 0 -0.75642 21.43515 ## 6 1 6 -2.2389995 0 -0.75642 20.00458
Once I’ve generated responses, I always visualize and summarize the response distribution in each arm to ensure it matches my expectations. The distributions suggest a moderate effect size (enough shift to tell a difference, but still mostly overlapping) and estimated means align with our known parameter values.
# STEP 6: Visualize and summarize distribution to ensure correctness -----------
library(ggplot2) |> suppressWarnings() # for data viz
df$label <- ifelse(df$arm==1, "Treatment", "Control") # for a tidy legend
df |> # visualize arm-specific distributions using density
ggplot(aes(x=Y, col=label, fill=label)) +
geom_density(alpha=0.5) +
labs(x="Generated Response (Y)", y="Density", col="Arm", fill="Arm") +
theme_bw() +
theme(legend.position="bottom")
# Summary of responses in the treatment arm:
summary(df$Y[df$arm==1])
## Min. 1st Qu. Median Mean 3rd Qu. Max. ## 9.329 22.475 25.248 25.190 28.038 38.810
# Summary of responses in the control arm:
summary(df$Y[df$arm==0])
## Min. 1st Qu. Median Mean 3rd Qu. Max. ## 8.043 20.158 22.895 22.878 25.606 35.913
Now we can analyze our generated data as if it arised from a cluster randomized trial! We fit a linear mixed model with a random intercept for each cid using the lme function within the nlme library. Our analysis model aligns perfectly with our data generating model, so we should see good performance here. The summary of our model provides model fit statistics, between (intercept) and within-cluster (residual) variance components, and fixed effects (regression coefficients). The coefficient of arm should correspond to the ATE.
Note that the degrees of freedom for the treatment effect, i.e., coefficient of arm, equals 
. There is no unique way of specifying the degrees of freedom within a linear mixed model. may be fine in large samples, but other methods perform better in small samples, like Satterthwaite degrees of freedom.
# STEP 7: Fit linear mixed model to estimate ATE and variance components
library(nlme) # to fit linear mixed models
LMM <- lme(Y~arm, random=~1|cid, data=df) # fixed effects, random effect, data
# Summary of LMM:
summary(LMM)
## Linear mixed-effects model fit by REML ## Data: df ## AIC BIC logLik ## 55819.35 55848.19 -27905.67 ## ## Random effects: ## Formula: ~1 | cid ## (Intercept) Residual ## StdDev: 1.285042 3.793098 ## ## Fixed effects: Y ~ arm ## Value Std.Error DF t-value p-value ## (Intercept) 22.877611 0.07877184 9000 290.42881 0 ## arm 2.312229 0.11117807 998 20.79753 0 ## Correlation: ## (Intr) ## arm -0.709 ## ## Standardized Within-Group Residuals: ## Min Q1 Med Q3 Max ## -4.050530724 -0.640817269 -0.001490679 0.669057102 3.439998952 ## ## Number of Observations: 10000 ## Number of Groups: 1000
Finally, we can extract values from our LMM object to calculate confidence intervals for the average treatment effect, fixed and random variance components, and the intracluster correlation.
# STEP 8: Extract and compare estimates to parameters
# View variance-covariance matrix for LMM coefficients:
vcov(LMM)
## (Intercept) arm ## (Intercept) 0.006205002 -0.006205002 ## arm -0.006205002 0.012360562
# View random intercept and residual variance estimates:
VarCorr(LMM)
## cid = pdLogChol(1) ## Variance StdDev ## (Intercept) 1.651332 1.285042 ## Residual 14.387590 3.793098
# Fixed effects:
est_ATE <- LMM$coefficients$fixed['arm']
est_ATE_df <- LMM$fixDF$X['arm']
est_ATE_se <- sqrt(vcov(LMM)['arm', 'arm'])
est_ATE_ci <- est_ATE + c(-1, 1) * qt(0.975, df=est_ATE_df) * est_ATE_se
fixed_text <- paste0("The estimated ATE is ", r2(est_ATE),
" with SE ", r2(est_ATE_se), " on ", est_ATE_df, " df.",
" The corresponding 95% CI for the ATE is: ",
r2(est_ATE_ci[1]), " to ", r2(est_ATE_ci[2]), ".")
print(fixed_text)
## [1] "The estimated ATE is 2.31 with SE 0.11 on 998 df. The corresponding 95% CI for the ATE is: 2.09 to 2.53."
# Random effects:
est_var_a <- as.numeric(VarCorr(LMM)[1,'Variance'])
est_var_e <- as.numeric(VarCorr(LMM)[2,'Variance'])
est_var <- est_var_a + est_var_e
est_icc <- est_var_a / est_var
random_text <- paste0("The estimated random intercept variance is ", r2(est_var_a),
" and the estimated residual variance is ", r2(est_var_e),
". The estimated ICC is ", r2(est_icc) , ".")
print(random_text)
## [1] "The estimated random intercept variance is 1.65 and the estimated residual variance is 14.39. The estimated ICC is 0.1."
Huzzah!!! Everything aligns with our parameter values: ATE = 2, random intercept variance = 1.6, random residual variance = 14.4, and ICC = 0.1.
Partially-nested individually randomized group treatment trial
Consider a two-arm, partially-nested cluster randomized trial featuring a group-based intervention in the treatment arm and individual intervention, or no groups, in the control arm.
Data generation for a partially-nested IRGT is similar to that of a cluster trial, with two twists. First, the intracluster correlation should be non-zero in the treatment arm and zero in the individual control arm. This means we now need to specify two values of the ICC and two sets of variance components. Second, the number of clusters in the treatment arm is equal to the number of groups (size 
). The number of clusters in the control arm is equal to the number of participants, or group size 
. That is, each control participant is considered their own group or cluster.
Suppose there are groups in the treatment arm, each featuring 
participants, and 
participants in the control arm. In total, there are 
participants in the study. We will generate normally distributed responses. We assume the total variance in each arm is the same, even if the intracluster correlation differs.
Assume we want 25 groups of 10 individuals in the group treatment arm and 250 individuals in the control arm.
# STEP 0: Specify important design and response parameters ---------------------
## Seed for reproducibility
set.seed(1234)
## Helper function
r2 <- function(x) round(x,2) # function to clean up reported numbers
## Design
N = 500 # total number of individuals to randomize
C1 = 25 # number of treatment groups
n1 = 10 # expected number of participants per group in trt arm
n0 = 1 # expected number of participants per group in ctrl arm
C0 = N - C1*n1 # expected number of control clusters
C = C1 + C0 # expected number of total clusters
N1 = C1*n1 # expected number of total participants in trt arm
N0 = C0*n0 # expected number of total participants in ctrl arm
## GROUP TREATMENT ARM: Intracluster correlation and variance
icc = 0.1 # intracluster correlation = between / (between + within)
var = 16 # total variance = between + within cluster variance
var_a = icc * var # between-cluster variance (random intercept)
var_e = var - var_a # within-cluster variance (residual)
## INDIVIDUAL CONTROL ARM:
var0 = 16 # total variance = between + within cluster variance
## Effects and means
ES = 0.5 # cohen's effect size (0.2=small, 0.5=moderate, 0.8=large)
mu1 = 25 # mean in treatment arm (fixed)
mu0 = mu1 - ES*sqrt(var) # mean in control arm (mu1 - ES*sd)
ATE = mu1 - mu0 # average treatment effect
## Summarize parameter values
param_text <- paste0("You are generating a partially-nested IRGT with ",
C1, " groups of ", n1, " participants in the group treatment arm ",
"and ", N0, " participants in the individual control arm.",
" In the group treatment arm, the intracluster correlation is ", icc,
" and the total variance is ", var,
", corresponding to a random intercept variance of ", var_a,
" and residual variance ", var_e, ". The effect size is ", ES,
" and the treatment mean is ", mu1,
", corresponding to a control mean of ", r2(mu0),
". The average treatment effect (ATE) is ", r2(ATE),".")
print(param_text)
## [1] "You are generating a partially-nested IRGT with 25 groups of 10 participants in the group treatment arm and 250 participants in the individual control arm. In the group treatment arm, the intracluster correlation is 0.1 and the total variance is 16, corresponding to a random intercept variance of 1.6 and residual variance 14.4. The effect size is 0.5 and the treatment mean is 25, corresponding to a control mean of 23. The average treatment effect (ATE) is 2."
We need to assign a unique ID to each cluster (cid) and to each participant (pid) within a given cid. Notice that we need to do this separately by arm now, with each control participant being assigned their own cluster ID and participant ID.
# STEP 1: Initialize dataframe structure to store simulated cluster responses --
df <- data.frame(cid = c(rep(1:C1, each=n1), rep((C1+1):C, each=1)), # trt, ctrl
pid = c(rep(1:n1, times=C1), rep(1:n0, times=C0))) # trt, ctrl
Each participant has an associated residual, but the residual variance now differs by arm. So, we generate residuals separately by arm, and then combine.
Note that we have assigned IDs so that all treatment responses are stacked on top of all control responses. i.e., the first 
rows correspond to the group treatment arm and last 
rows correspond to the individual control arm.
# STEP 2: Generate residual for each participant, epsilon ~ N(0, var_e) --------
df$epsilon <- c(rnorm(n=N1, mean=0, sd=sqrt(var_e)), # group treatment
rnorm(n=N0, mean=0, sd=sqrt(var0))) # individual control
In the group treatment arm, we generate random intercepts according to where is the between-cluster variance. In the individual control arm, we set all random intercept values to 0.
The way we’ve set this up, all groups must be assigned to the treatment arm and all singletons to the control arm. In reality, you might accumulate a pool of eligible individuals at a given site and then individually randomize them to be in either arm and form groups accordingly. e.g., If a site has a pool of 20 individuals, truly random 1:1 allocation may randomize 11 participants to treatment and 9 to control. Then, there may be one treatment group with 11 individuals and 9 control individuals.
The way we have assigned groups here should approximate behavior well, unless large variation in group sizes is permitted. In my experience, group sizes must fall within a specified range. For example, if our desired group size is 10 individuals, we may allow sites to form groups of 8 to 12 individuals to permit drop-out or unequal numbers of individuals allocated to each arm, for example.
# STEP 3: Generate cluster intercepts, alpha ~ N(0, var_a), and randomized arm -
df_a <- data.frame(cid = 1:C,
arm = c(rep(1, C1), rep(0, C0)),
alpha = c(rnorm(n=C1, mean=0, sd=sqrt(var_a)),
rep(0, C0)))
cluster_text <- paste0("C1=", sum(df_a$arm),
" treatment groups. ",
"C0=", C-sum(df_a$arm),
" control participants.")
cluster_text
## [1] "C1=25 treatment groups. C0=250 control participants."
Now, we need to attach the random cluster intercept to each participant within cluster cid.
# STEP 4: Attach random intercept to each participant within cid ---------------
df <- merge(df, df_a, by='cid')
We assume a random intercept structure for our generated data, as before. We’ve done the hard work of simulating the random components of this equation, now we just need to add them all together! Voila, generated responses!
# STEP 5: Generate responses according to random intercept structure -----------
df$Y <- mu0 + ATE * df$arm + df$alpha + df$epsilon
## Preview dataframe of generated responses:
head(df)
## cid pid epsilon arm alpha Y ## 1 1 1 -4.580492 1 1.245659 21.66517 ## 2 1 2 1.052770 1 1.245659 27.29843 ## 3 1 3 4.115165 1 1.245659 30.36082 ## 4 1 4 -8.901297 1 1.245659 17.34436 ## 5 1 5 1.628414 1 1.245659 27.87407 ## 6 1 6 1.920347 1 1.245659 28.16601
Once I’ve generated responses, I always visualize and summarize the response distribution in each arm to ensure it matches my expectations. The distributions suggest a moderate effect size (enough shift to tell a difference, but still mostly overlapping) and estimated means align with our known parameter values.
# STEP 6: Visualize and summarize distribution to ensure correctness -----------
library(ggplot2) |> suppressWarnings() # for data viz
df$label <- ifelse(df$arm==1, "Treatment", "Control") # for a tidy legend
df |> # visualize arm-specific distributions using density
ggplot(aes(x=Y, col=label, fill=label)) +
geom_density(alpha=0.5) +
labs(x="Generated Response (Y)", y="Density", col="Arm", fill="Arm") +
theme_bw() +
theme(legend.position="bottom")
# Summary of responses in the treatment arm:
summary(df$Y[df$arm==1])
## Min. 1st Qu. Median Mean 3rd Qu. Max. ## 12.73 22.51 24.50 25.01 27.55 37.72
# Summary of responses in the control arm:
summary(df$Y[df$arm==0])
## Min. 1st Qu. Median Mean 3rd Qu. Max. ## 9.416 20.218 23.076 23.031 25.876 35.784
Now we can fit a linear mixed model! However, our random effects structure is now much more complicated. Although we have specified the same total variance 
for both arms, their within and between-group variances differ. We need our model to allow for heterogeneous between-group (random intercept) variances and heterogeneous within-group (residual) variances.
To specify heterogeneous between-group variances, we need to specify a covariance matrix with diagonal entries equal to 
for the treatment arm and 
for the control arm. Off-diagonal entries of this matrix are 0, because clusters are independent due to randomization. We do this with an expanded random statement that specifies a diagonal structure (pdDiag) with entries by arm.
To specify heterogeneous within-group variances, we need to let the residual variance equal 
for the treatment arm and for the control arm. We do this through the weights statement. Note: Specifying and interpreting these flexible structures in lme is not particularly straight-forward. I found this blog post by The Change Lab very helpful in this regard.
The summary of LMM appears similar to before, but now we have an additional variance function section. Values here represent multipliers of the residual standard deviation for each arm. That is, the treatment residual SD would be (residual SD x 1) whereas for control it would be (residual SD x 1.09), demonstrated shortly.
# STEP 7: Fit linear mixed model to estimate ATE and variance components
library(nlme) # to fit linear mixed models
LMM <- lme(Y~arm,
random=list(cid = pdDiag(form = ~arm)), # between-cluster var by arm
weights=varIdent(form = ~1 | arm), # within-cluster var by arm
data=df)
# Summary of LMM:
summary(LMM)
## Linear mixed-effects model fit by REML ## Data: df ## AIC BIC logLik ## 2830.109 2855.373 -1409.055 ## ## Random effects: ## Formula: ~arm | cid ## Structure: Diagonal ## (Intercept) arm Residual ## StdDev: 1.092631 0.8838134 3.733902 ## ## Variance function: ## Structure: Different standard deviations per stratum ## Formula: ~1 | arm ## Parameter estimates: ## 1 0 ## 1.000000 1.090465 ## Fixed effects: Y ~ arm ## Value Std.Error DF t-value p-value ## (Intercept) 23.03103 0.2666270 273 86.37920 0 ## arm 1.98229 0.4537145 273 4.36902 0 ## Correlation: ## (Intr) ## arm -0.588 ## ## Standardized Within-Group Residuals: ## Min Q1 Med Q3 Max ## -3.11926915 -0.62385663 -0.03420332 0.63052705 3.02485189 ## ## Number of Observations: 500 ## Number of Groups: 275
Finally, we can extract values from our LMM object to calculate confidence intervals for the average treatment effect, fixed and random variance components, and the intracluster correlation. There’s more to extract here, now that we have heterogeneous variance components.
# STEP 8: Extract and compare estimates to parameters
# View variance-covariance matrix for LMM coefficients:
vcov(LMM)
## (Intercept) arm ## (Intercept) 0.07108996 -0.07108996 ## arm -0.07108996 0.20585681
# View random intercept and residual variance estimates:
VarCorr(LMM)
## cid = pdDiag(arm) ## Variance StdDev ## (Intercept) 1.1938431 1.0926313 ## arm 0.7811261 0.8838134 ## Residual 13.9420215 3.7339016
# Fixed effects:
est_ATE <- LMM$coefficients$fixed['arm']
est_ATE_df <- LMM$fixDF$X['arm']
est_ATE_se <- sqrt(vcov(LMM)['arm', 'arm'])
est_ATE_ci <- est_ATE + c(-1, 1) * qt(0.975, df=est_ATE_df) * est_ATE_se
fixed_text <- paste0("The estimated ATE is ", r2(est_ATE),
" with SE ", r2(est_ATE_se), " on ", est_ATE_df, " df.",
" The corresponding 95% CI for the ATE is: ",
r2(est_ATE_ci[1]), " to ", r2(est_ATE_ci[2]), ".")
print(fixed_text)
## [1] "The estimated ATE is 1.98 with SE 0.45 on 273 df. The corresponding 95% CI for the ATE is: 1.09 to 2.88."
# Random effects:
resid_wt0 <- coef(LMM$modelStruct$varStruct, unconstrained=F)
est_var_a <- as.numeric(VarCorr(LMM)['(Intercept)','Variance'])
est_var_e <- as.numeric(VarCorr(LMM)['Residual','Variance'])
est_var0 <- est_var_e * resid_wt0^2
est_var <- est_var_a + est_var_e
est_icc <- est_var_a / est_var
random_text <- paste0("Within the group treatment arm, ",
"the estimated random intercept variance is ", r2(est_var_a),
" and the estimated residual variance is ", r2(est_var_e),
". The estimated ICC is ", r2(est_icc) ,
". Within the individual control arm, the estimated ",
"residual variance is ", r2(est_var0),".")
print(random_text)
## [1] "Within the group treatment arm, the estimated random intercept variance is 1.19 and the estimated residual variance is 13.94. The estimated ICC is 0.08. Within the individual control arm, the estimated residual variance is 16.58."
Huzzah!!! Everything aligns with our parameter values. Everything aligns with our parameter values: ATE = 2, treatment random intercept variance = 1.6, treatment residual variance = 14.4, treatment ICC = 0.1, and control residual variance 16.
Notice that there are many more parameters to estimate for a partially-nested individually randomized group treatment trial, so it is important to ensure you have enough individuals to precisely capture these quantities! Simulation can help you figure this out. 🙂 If this feels like a nightmare to you, the NIH also provides lovely, easy-to-use sample size calculators for fully and partially-nested IRGTs.
In closing…
So, admittedly, I’ve kind of cheated here… I used huge, unrealistic sample sizes to ensure that I could get close to the true parameter value with a single replicate, rather than performing multiple replicates to assess sampling distributions.
Cluster randomized trials and individually randomized group treatment trials can be very logistically challenging. In the context of cluster randomized trials, clusters like hospitals can be limited in number and difficult to recruit, making few large clusters a common reality. In the context of individually randomized group treatment trials, it can be hard to accumulate large pools of eligible participants for randomization. If it takes too long to accumulate a pool, drop-out is a real threat to study feasibility. As a result, the numbers of clusters and groups used here are not realistic.
With smaller sample sizes, there will be greater variability in effect estimates. Especially because clustered designs will always be less efficient than traditional individually randomized trials due to correlation. A single replicate with realistic sample sizes can yield estimates that are not particularly close to the truth, even if our model is correctly specified.
The data generation processes described here can be turned into a function, for example, and used to allow for multiple replicates to be constructed and multiple estimates obtained. This allows you to vary sample sizes and parameters to investigate how estimates vary, and the distribution of results you may reasonably expect if you were to conduct such a study. This is very helpful for sample size or power calculation, for example.
Thanks for reading!!!
