Continuing a series of posts discussing the structure of intra-cluster correlations (ICC’s) in the context of a stepped-wedge trial, this latest edition is primarily interested in fitting Bayesian hierarchical models for more complex cases (though I do talk a bit more about the linear mixed effects models). The first two posts in the series focused on generating data to simulate various scenarios; the third post considered linear mixed effects and Bayesian hierarchical models to estimate ICC’s under the simplest scenario of constant between-period ICC’s. Throughout this post, I use code drawn from the previous one; I am not repeating much of it here for brevity’s sake. So, if this is all new, it is probably worth glancing at before continuing on.

Data generation

The data generating model this time around is only subtly different from before, but that difference is quite important. Rather than a single cluster-specific effect \(b_c\), there is now a vector of cluster effects \(\mathbf{b_c} = \left( b_{c1}, b_{c2}, \ldots, b_{cT} \right)\), where \(b_c \sim MVN(\mathbf{0}, \sigma^2 \mathbf{R})\) (see this earlier post for a description of the correlation matrix \(\mathbf{R}\).)

\[ Y_{ict} = \mu + \beta_0t + \beta_1X_{ct} + b_{ct} + e_{ict} \]

By altering the correlation structure of \(\mathbf{b_c}\) (that is \(\mathbf{R}\)), we can the change the structure of the ICC’s. (The data generation was the focus of the first two posts of this series, here and here. The data generating function genDD includes an argument where you can specify the two correlation structures, exchangeable and auto-regressive:

library(simstudy)

defc <- defData(varname = "mu", formula = 0, 
                dist = "nonrandom", id = "cluster")
defc <- defData(defc, "s2", formula = 0.15, dist = "nonrandom")
defc <- defData(defc, "m", formula = 15, dist = "nonrandom")

defa <- defDataAdd(varname = "Y", 
                   formula = "0 + 0.10  * period + 1 * rx + cteffect", 
                   variance = 2, dist = "normal")
genDD <- function(defc, defa, nclust, nperiods, 
                  waves, len, start, rho, corstr) {
  
  dc <- genData(nclust, defc)
  dp <- addPeriods(dc, nperiods, "cluster")
  dp <- trtStepWedge(dp, "cluster", nWaves = waves, lenWaves = len, 
                     startPer = start)
  dp <- addCorGen(dtOld = dp, nvars = nperiods, idvar = "cluster", 
                  rho = rho, corstr = corstr, dist = "normal", 
                  param1 = "mu", param2 = "s2", cnames = "cteffect")
  
  dd <- genCluster(dp, cLevelVar = "timeID", numIndsVar = "m", 
                   level1ID = "id")
  dd <- addColumns(defa, dd)
  dd[]
}

Constant between-period ICC’s

In this first scenario, the assumption is that the within-period ICC’s are larger than the between-period ICC’s and the between-period ICC’s are constant. This can be generated with random effects that have a correlation matrix with compound symmetry (or is exchangeable). In this case, we will have 60 clusters and 7 time periods:

set.seed(4119)
dcs <- genDD(defc, defa, 60, 7, 4, 1, 2, 0.6, "cs")

# correlation of "unobserved" random effects

round(cor(dcast(dcs[, .SD[1], keyby = .(cluster, period)], 
  formula = cluster ~ period, value.var = "cteffect")[, 2:7]), 2)
##      0    1    2    3    4    5
## 0 1.00 0.60 0.49 0.60 0.60 0.51
## 1 0.60 1.00 0.68 0.64 0.62 0.64
## 2 0.49 0.68 1.00 0.58 0.54 0.62
## 3 0.60 0.64 0.58 1.00 0.63 0.66
## 4 0.60 0.62 0.54 0.63 1.00 0.63
## 5 0.51 0.64 0.62 0.66 0.63 1.00


Linear mixed-effects model

It is possible to use lmer to correctly estimate the variance components and other parameters that underlie the data generating process used in this case. The cluster-level period-specific effects are specified in the model as “cluster/period”, which indicates that the period effects are nested within the cluster.

library(lme4)

lmerfit <- lmer(Y ~ period + rx + (1 | cluster/period) , data = dcs)
as.data.table(VarCorr(lmerfit))
##               grp        var1 var2       vcov     sdcor
## 1: period:cluster (Intercept) <NA> 0.05827349 0.2413990
## 2:        cluster (Intercept) <NA> 0.07816476 0.2795796
## 3:       Residual        <NA> <NA> 2.02075355 1.4215321

Reading from the vcov column in the lmer output above, we can extract the period:cluster variance (\(\sigma_w^2\)), the cluster variance (\(\sigma^2_v\)), and the residual (individual level) variance (\(\sigma^2_e\)). Using these three variance components, we can estimate the correlation of the cluster level effects (\(\rho\)), the within-period ICC (\(ICC_{tt}\)), and the between-period ICC (\(ICC_{tt^\prime}\)). (See the addendum below for a more detailed description of the derivations.)

Correlation (\(\rho\)) of cluster-specific effects over time

In this post, don’t confuse \(\rho\) with the ICC. \(\rho\) is the correlation between the cluster-level period-specific random effects. Here I am just showing that it is function of the decomposed variance estimates provided in the lmer output:

\[ \rho = \frac{\sigma^2_v}{\sigma^2_v + \sigma^2_w} \]

vs <- as.data.table(VarCorr(lmerfit))$vcov
vs[2]/sum(vs[1:2])  
## [1] 0.5728948


Within-period ICC

The within-period ICC is the ratio of total cluster variance relative to total variance:

\[ICC_{tt} = \frac{\sigma^2_v + \sigma^2_w}{\sigma^2_v + \sigma^2_w+\sigma^2_e}\]

sum(vs[1:2])/sum(vs)
## [1] 0.06324808


Between-period ICC

The between-period \(ICC_{tt^\prime}\) is really just the within-period \(ICC_{tt}\) adjusted by \(\rho\) (see the addendum):

\[ICC_{tt^\prime} = \frac{\sigma^2_v}{\sigma^2_v + \sigma^2_w+\sigma^2_e}\]

vs[2]/sum(vs)          
## [1] 0.0362345


Bayesian model

Now, I’ll fit a Bayesian hierarchical model, as I did earlier with the simplest constant ICC data generation process. The specification of the model in stan in this instance is slightly more involved as the number of parameters has increased. In the simpler case, I only had to estimate a scalar parameter for \(\sigma_b\) and a single ICC parameter. In this model definition (nested_cor_cs.stan) \(\mathbf{b_c}\) is a vector so there is a need to specify the variance-covariance matrix \(\sigma^2 \mathbf{R}\), which has dimensions \(T \times T\) (defined in the transformed parameters block). There are \(T\) random effects for each cluster, rather than one. And finally, instead of one ICC value, there are two - the within- and between-period ICC’s (defined in the generated quantities block).

data {
  int<lower=0> N;              // number of unique individuals
  int<lower=1> K;              // number of predictors
  int<lower=1> J;              // number of clusters
  int<lower=0> T;              // number of periods
  int<lower=1,upper=J> jj[N];  // group for individual
  int<lower=1> tt[N];          // period for individual
  matrix[N, K] x;              // matrix of predctors
  vector[N] y;                 // matrix of outcomes
}

parameters {
  vector[K] beta;              // model fixed effects
  real<lower=0> sigmalev1;     // cluster variance (sd)
  real<lower=-1,upper=1> rho;  // correlation
  real<lower=0> sigma;         // individual level varianc (sd)
  
  matrix[J, T] ran;            // site level random effects (by period)
}
 
transformed parameters{ 
  
  cov_matrix[T] Sigma;
  vector[N] yhat;
  vector[T] mu0;
  
  for (t in 1:T) 
    mu0[t] = 0;
    
  // Random effects with exchangeable correlation  
  
  for (j in 1:(T-1))
    for (k in (j+1):T) {
      Sigma[j,k] = pow(sigmalev1,2) * rho;  
      Sigma[k,j] = Sigma[j, k];
    }
     
  for (i in 1:T)
      Sigma[i,i] = pow(sigmalev1,2);
  
  for (i in 1:N)  
      yhat[i] = x[i]*beta + ran[jj[i], tt[i]];
}

model {
  
  sigma ~ uniform(0, 10);
  sigmalev1 ~ uniform(0, 10);
  rho ~ uniform(-1, 1);
  
  for (j in 1:J)
    ran[j] ~ multi_normal(mu0, Sigma);
    
  y ~ normal(yhat, sigma);

}

generated quantities {
  
  real sigma2;
  real sigmalev12;
  real iccInPeriod;
  real iccBetPeriod;
  
  sigma2 = pow(sigma, 2);
  sigmalev12 = pow(sigmalev1, 2);
  
  iccInPeriod = sigmalev12/(sigmalev12 + sigma2);
  iccBetPeriod = iccInPeriod * rho;
  
}

Model estimation requires creating the data set (in the form of an R list), compiling the stan model, and then sampling from the posterior to generate distributions of all parameters and generated quantities. I should include conduct a diagnostic review (e.g. to assess convergence), but you’ll have to trust me that everything looked reasonable.

library(rstan)
options(mc.cores = parallel::detectCores())

x <- as.matrix(dcs[ ,.(1, period, rx)])
K <- ncol(x)
N <- dcs[, length(unique(id))]
J <- dcs[, length(unique(cluster))]
T <- dcs[, length(unique(period))]
jj <- dcs[, cluster]
tt <- dcs[, period] + 1
y <- dcs[, Y]

testdat <- list(N=N, K=K, J=J, T=T, jj=jj, tt=tt, x=x, y=y)

rt <- stanc("nested_cor_cs.stan")
sm <- stan_model(stanc_ret = rt, verbose=FALSE)
fit.cs <- sampling(sm, data=testdat, seed = 32748, 
                  iter = 5000, warmup = 1000,
                  control = list(max_treedepth = 15))

Here is a summary of results for \(\rho\), \(ICC_{tt}\), and \(ICC_{tt^\prime}\). I’ve included a comparison of the means of the posterior distributions with the lmer estimates, followed by a more complete (visual) description of the posterior distributions of the Bayesian estimates:

mb <- sapply(
  rstan::extract(fit.cs, pars=c("rho", "iccInPeriod", "iccBetPeriod")),
  function(x) mean(x)
)

cbind(bayesian=round(mb,3), 
      lmer = round(c(vs[2]/sum(vs[1:2]), 
                     sum(vs[1:2])/sum(vs), 
                     vs[2]/sum(vs)),3)
)
##              bayesian  lmer
## rho             0.576 0.573
## iccInPeriod     0.065 0.063
## iccBetPeriod    0.037 0.036

Decaying between-period ICC over time

Now we enter somewhat uncharted territory, since there is no obvious way in R using the lme4 or nlme packages to decompose the variance estimates when the random effects have correlation that decays over time. This is where we might have to rely on a Bayesian approach to do this. (I understand that SAS can accommodate this, but I can’t bring myself to go there.)

We start where we pretty much always do - generating the data. Everything is the same, except that the cluster-random effects are correlated over time; we specify a correlation structure of ar1 (auto-regressive).

set.seed(4119)
dar1 <- genDD(defc, defa, 60, 7, 4, 1, 2, 0.6, "ar1")

# correlation of "unobserved" random effects

round(cor(dcast(dar1[, .SD[1], keyby = .(cluster, period)], 
  formula = cluster ~ period, value.var = "cteffect")[, 2:7]), 2)
##      0    1    2    3    4    5
## 0 1.00 0.60 0.22 0.20 0.18 0.06
## 1 0.60 1.00 0.64 0.45 0.30 0.23
## 2 0.22 0.64 1.00 0.61 0.32 0.30
## 3 0.20 0.45 0.61 1.00 0.61 0.49
## 4 0.18 0.30 0.32 0.61 1.00 0.69
## 5 0.06 0.23 0.30 0.49 0.69 1.00

The model file is similar to nested_cor_cs.stan, except that the specifications of the variance-covariance matrix and ICC’s are now a function of \(\rho^{|t^\prime - t|}\):

transformed parameters{ 
  ⋮
  for (j in 1:T)
    for (k in 1:T)
      Sigma[j,k] = pow(sigmalev1,2) * pow(rho,abs(j-k));
  ⋮
}

generated quantities {
  ⋮
  for (j in 1:T)
      for (k in 1:T)
        icc[j, k] = sigmalev12/(sigmalev12 + sigma2) * pow(rho,abs(j-k));
  ⋮
}

The stan compilation and sampling code is not shown here - they are the same before. The posterior distribution of \(\rho\) is similar to what we saw previously.

print(fit.ar1, pars=c("rho"))
## Inference for Stan model: nested_cor_ar1.
## 4 chains, each with iter=5000; warmup=1000; thin=1; 
## post-warmup draws per chain=4000, total post-warmup draws=16000.
## 
##     mean se_mean   sd 2.5%  25%  50%  75% 97.5% n_eff Rhat
## rho 0.58       0 0.08 0.41 0.53 0.58 0.64  0.73  2302    1
## 
## Samples were drawn using NUTS(diag_e) at Fri Jun 28 14:13:41 2019.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at 
## convergence, Rhat=1).

Now, however, we have to consider a full range of ICC estimates. Here is a plot of the posterior distribution of all ICC’s with the means of each posterior directly below. The diagonal represents the within-period (constant) ICCs, and the off-diagonals are the between-period ICC’s.

An alternative Bayesian model: unstructured correlation

Now, there is no particular reason to expect that the particular decay model (with an AR1 structure) would be the best model. We could try to fit an even more general model, one with minimal structure. For example if we put no restrictions on the correlation matrix \(\mathbf{R}\), but assumed a constant variance of \(\sigma_b^2\), we might achieve a better model fit. (We could go even further and relax the assumption that the variance across time changes as well, but I’ll leave that to you if you want to try it.)

In this case, we need to define \(\mathbf{R}\) and specify a prior distribution (I use the Lewandowski, Kurowicka, and Joe - LKJ prior, as suggested by Stan documentation) and define the ICC’s in terms of \(\mathbf{R}\). Here are the relevant snippets of the stan model (everything else is the same as before):

parameters {
  ⋮
  corr_matrix[T] R;        // correlation matrix
  ⋮
}
 
transformed parameters{ 
  ⋮
  Sigma = pow(sigmalev1,2) * R;
  ⋮
}

model {
  ⋮
  R ~ lkj_corr(1);         // LKJ prior on the correlation matrix 
  ⋮
}

generated quantities {
  ⋮
  for (j in 1:T)
    for (k in 1:T)
      icc[j, k] = sigmalev12/(sigmalev12 + sigma2) * R[j, k];
  ⋮
}

Here are the means of the ICC posterior distributions alongside the means from the previous auto-regressive model.

Looking at the unstructured model estimates on the right, it does appear that a decay model might be reasonable. (No surprise there, because in reality, it is totally reasonable; that’s how we generated the data.) We can use package bridgesampling which estimates marginal log likelihoods (across the prior distributions of the parameters). The marginal likelihoods are used in calculating the Bayes Factor, which is the basis for comparing two competing models. Here, the log-likelihood is reported. If the unstructured model is indeed an improvement (and it could very well be, because it has more parameters), the we would expect the marginal log-likelihood for the second model to be greater (less negative) than the log-likelihood for the auto-regressive model. If fact, the opposite true, suggesting the auto-regressive model is the preferred one (out of these two):

library(bridgesampling)

bridge_sampler(fit.ar1, silent = TRUE)
## Bridge sampling estimate of the log marginal likelihood: -5132.277
## Estimate obtained in 6 iteration(s) via method "normal"
bridge_sampler(fit.ar1.nc, silent = TRUE)
## Bridge sampling estimate of the log marginal likelihood: -5137.081
## Estimate obtained in 269 iteration(s) via method "normal".

 

 

Addendum - interpreting lmer variance estimates

In order to show how the lmer variance estimates relate to the theoretical variances and correlations in the case of a constant between-period ICC, here is a simulation based on 1000 clusters. The key parameters are \(\sigma^2_b = 0.15\), \(\sigma^2_e = 2\), and \(\rho = 0.6\). And based on these values, the theoretical ICC’s are: \(ICC_{within} = 0.15/2.15 = 0.698\), and \(ICC_{bewteen} = 0.698 * 0.6 = 0.042\).

set.seed(4119)
dcs <- genDD(defc, defa, 1000, 7, 4, 1, 2, 0.6, "cs")

The underlying correlation matrix of the cluster-level effects is what we would expect:

round(cor(dcast(dcs[, .SD[1], keyby = .(cluster, period)], 
  formula = cluster ~ period, value.var = "cteffect")[, 2:7]), 2)
##      0    1    2    3    4    5
## 0 1.00 0.59 0.59 0.61 0.61 0.61
## 1 0.59 1.00 0.61 0.60 0.61 0.64
## 2 0.59 0.61 1.00 0.59 0.61 0.61
## 3 0.61 0.60 0.59 1.00 0.59 0.62
## 4 0.61 0.61 0.61 0.59 1.00 0.60
## 5 0.61 0.64 0.61 0.62 0.60 1.00

Here are the variance estimates from the mixed-effects model:

lmerfit <- lmer(Y ~ period + rx + (1 | cluster/period) , data = dcs)
as.data.table(VarCorr(lmerfit))
##               grp        var1 var2       vcov     sdcor
## 1: period:cluster (Intercept) <NA> 0.05779349 0.2404028
## 2:        cluster (Intercept) <NA> 0.09143749 0.3023863
## 3:       Residual        <NA> <NA> 1.98894356 1.4102991

The way lmer implements the nested random effects , the cluster period-specific effect \(b_{ct}\) is decomposed into \(v_c\), a cluster level effect, and \(w_{ct}\), a cluster time-specific effect:

\[ b_{ct} = v_c + w_{ct} \]

Since both \(v_c\) and \(w_{ct}\) are normally distributed (\(v_c \sim N(0,\sigma_v^2)\) and \(w_{ct} \sim N(0,\sigma_w^2)\)), \(var(b_{ct}) = \sigma^2_b = \sigma^2_v + \sigma^2_w\).

Here is the observed estimate of \(\sigma^2_v + \sigma^2_w\):

vs <- as.data.table(VarCorr(lmerfit))$vcov
sum(vs[1:2])
## [1] 0.149231

An estimate of \(\rho\) can be extracted from the lmer model variance estimates:

\[ \begin{aligned} \rho &= cov(b_{ct}, b_{ct^\prime}) \\ &= cov(v_{c} + w_{ct}, v_{c} + w_{ct^\prime}) \\ &= var(v_c) + cov(w_{ct}) \\ &= \sigma^2_v \end{aligned} \]

\[ \begin{aligned} var(b_{ct}) &= var(v_{c}) + var(w_{ct}) \\ &= \sigma^2_v + \sigma^2_w \end{aligned} \]

\[ \begin{aligned} cor(b_{ct}, b_{ct^\prime}) &= \frac{cov(b_{ct}, b_{ct^\prime})}{\sqrt{var(b_{ct}) var(b_{ct^\prime})} } \\ \rho &= \frac{\sigma^2_v}{\sigma^2_v + \sigma^2_w} \end{aligned} \]

vs[2]/sum(vs[1:2])
## [1] 0.6127246

And here are the estimates of within and between-period ICC’s:

\[ICC_{tt} = \frac{\sigma^2_b}{\sigma^2_b+\sigma^2_e} =\frac{\sigma^2_v + \sigma^2_w}{\sigma^2_v + \sigma^2_w+\sigma^2_e}\]

sum(vs[1:2])/sum(vs)
## [1] 0.06979364

\[ \begin{aligned} ICC_{tt^\prime} &= \left( \frac{\sigma^2_b}{\sigma^2_b+\sigma^2_e}\right) \rho \\ \\ &= \left( \frac{\sigma^2_v + \sigma^2_w}{\sigma^2_v + \sigma^2_w+\sigma^2_e}\right) \rho \\ \\ &=\left( \frac{\sigma^2_v + \sigma^2_w}{\sigma^2_v + \sigma^2_w+\sigma^2_e} \right) \left( \frac{\sigma^2_v}{\sigma^2_v + \sigma^2_w} \right) \\ \\ &= \frac{\sigma^2_v}{\sigma^2_v + \sigma^2_w+\sigma^2_e} \end{aligned} \]

vs[2]/sum(vs)
## [1] 0.04276428