Different within- and between-period ICC’s

In a scenario with constant within- and between-period ICC’s, the correlated data can be induced using a single cluster-level effect like \(b_c\) in this model:

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

More complexity can be added if, instead of a single cluster level effect, we have a vector of correlated cluster/time specific effects \(\mathbf{b_c}\). These cluster-specific random effects \((b_{c1}, b_{c2}, \ldots, b_{cT})\) replace \(b_c\), and the slightly modified data generating model is

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

The vector \(\mathbf{b_c}\) has a multivariate normal distribution \(N_T(0, \sigma^2_b \mathbf{R})\). This model assumes a common covariance structure across all clusters, \(\sigma^2_b \mathbf{R}\), where the general version of \(\mathbf{R}\) is

\[ \mathbf{R} = \left( \begin{matrix} 1 & r_{12} & r_{13} & \cdots & r_{1T} \\ r_{21} & 1 & r_{23} & \cdots & r_{2T} \\ r_{31} & r_{32} & 1 & \cdots & r_{3T} \\ \vdots & \vdots & \vdots & \vdots & \vdots \\ r_{T1} & r_{T2} & r_{T3} & \cdots & 1 \end{matrix} \right ) \]

Constant correlation over time

In this first case, the correlation between individuals in the same cluster but different time periods is less than the correlation between individuals in the same cluster and same time period. In other words, \(ICC_{tt} \ne ICC_{tt^\prime}\). However the between-period correlation is constant, or \(ICC_{tt^\prime}\) are constant for all \(t\) and \(t^\prime\). We have these correlations when \(r_0 = \rho\) and \(r = 1\), giving

\[ \mathbf{R} = \mathbf{R}(\rho, 1) = \left( \begin{matrix} 1 & \rho & \rho & \cdots & \rho \\ \rho & 1 & \rho & \cdots & \rho \\ \rho & \rho & 1 & \cdots & \rho \\ \vdots & \vdots & \vdots & \vdots & \vdots \\ \rho & \rho & \rho & \cdots & 1 \end{matrix} \right ) \]

To simulate under this scenario, I am setting \(\sigma_b^2 = 0.15\), \(\sigma_e^2 = 2.0\), and \(\rho = 0.6\). We would expect the following ICC’s:

\[ \begin{aligned} ICC_{tt} &= \frac{0.15}{0.15+2.00} = 0.0698 \\ \\ ICC_{tt^\prime} &= \frac{0.15}{0.15+2.00}\times0.6 = 0.0419 \end{aligned} \]

Here is the code to define and generate the data:

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

defa <- defDataAdd(varname = "Y", 
                   formula = "0 + 0.10  * period + 1 * rx + cteffect", 
                   variance = 2, dist = "normal")

dc <- genData(100, defc)
dp <- addPeriods(dc, 7, "cluster")
dp <- trtStepWedge(dp, "cluster", nWaves = 4, lenWaves = 1, startPer = 2)
dp <- addCorGen(dtOld = dp, nvars = 7, idvar = "cluster", 
                rho = 0.6, corstr = "cs", dist = "normal", 
                param1 = "mu", param2 = "s2", cnames = "cteffect")
  
dd <- genCluster(dp, cLevelVar = "timeID", numIndsVar = 100, 
                 level1ID = "id")
dd <- addColumns(defa, dd)

As I did in my previous post, I’ve generated 200 data sets, estimated the within- and between-period ICC’s for each data set, and computed the average for each. The plot below shows the expected values in gray and the estimated values in purple and green.

Declining correlation over time

In this second case, we make an assumption that the correlation between individuals in the same cluster degrades over time. Here, the correlation between two individuals in adjacent time periods is stronger than the correlation between individuals in periods further apart. That is \(ICC_{tt^\prime} > ICC_{tt^{\prime\prime}}\) if \(|t^\prime - t| < |t^{\prime\prime} - t|\). This structure can be created by setting \(r_0 = 1\) and \(r=\rho\), giving us an auto-regressive correlation matrix \(R\):

\[ \mathbf{R} = \mathbf{R}(1, \rho) = \left( \begin{matrix} 1 & \rho & \rho^2 & \cdots & \rho^{T-1} \\ \rho & 1 & \rho & \cdots & \rho^{T-2} \\ \rho^2 & \rho & 1 & \cdots & \rho^{T-3} \\ \vdots & \vdots & \vdots & \vdots & \vdots \\ \rho^{T-1} & \rho^{T-2} & \rho^{T-3} & \cdots & 1 \end{matrix} \right ) \]

I’ve generated data using the same variance assumptions as above. The only difference in this case is that the corstr argument in the call to addCorGen is “ar1” rather than “cs” (which was used above). Here are a few of the expected correlations:

\[ \begin{aligned} ICC_{t,t} &= \frac{0.15}{0.15+2.00} = 0.0698 \\ \\ ICC_{t,t+1} &= \frac{0.15}{0.15+2.00}\times 0.6^{1} = 0.0419 \\ \\ ICC_{t,t+2} &= \frac{0.15}{0.15+2.00}\times 0.6^{2} = 0.0251 \\ \\ \vdots \\ ICC_{t, t+6} &= \frac{0.15}{0.15+2.00}\times 0.6^{6} = 0.0032 \end{aligned} \]

And here is the code:

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

defa <- defDataAdd(varname = "Y", 
                   formula = "0 + 0.10  * period + 1 * rx + cteffect", 
                   variance = 2, dist = "normal")

dc <- genData(100, defc)
dp <- addPeriods(dc, 7, "cluster")
dp <- trtStepWedge(dp, "cluster", nWaves = 4, lenWaves = 1, startPer = 2)
dp <- addCorGen(dtOld = dp, nvars = 7, idvar = "cluster", 
                rho = 0.6, corstr = "ar1", dist = "normal", 
                param1 = "mu", param2 = "s2", cnames = "cteffect")
  
dd <- genCluster(dp, cLevelVar = "timeID", numIndsVar = 10, 
                 level1ID = "id")
dd <- addColumns(defa, dd)

And here are the observed average estimates (based on 200 datasets) alongside the expected values:

Random slope

In this last case, I am exploring what the ICC’s look like in the context of random effects model that includes a cluster-specific intercept \(b_c\) and a cluster-specific slope \(s_c\):

\[ Y_{ict} = \mu + \beta_0 t + \beta_1 X_{ct} + b_c + s_c t + e_{ict} \]

Both \(b_c\) and \(s_c\) are normally distributed with mean 0, and variances \(\sigma_b^2\) and \(\sigma_s^2\), respectively. (In this example \(\sigma_b^2\) and \(\sigma_s^2\) are uncorrelated, but that may not necessarily be the case.)

Because of the random slopes, the variance of the \(Y\)’s increase over time:

\[ var(Y_{ict}) = \sigma^2_b + t^2 \sigma^2_s + \sigma^2_e \]

The same is true for the within- and between-period covariances:

\[ \begin{aligned} cov(Y_{ict}, Y_{jct}) &= \sigma^2_b + t^2 \sigma^2_s \\ \\ cov(Y_{ict}, Y_{jct^\prime}) &= \sigma^2_b + tt^\prime \sigma^2_s \\ \end{aligned} \]

The ICC’s that follow from these various variances and covariances are:

\[ \begin{aligned} ITT_{tt} &= \frac{\sigma^2_b + t^2 \sigma^2_s}{\sigma^2_b + t^2 \sigma^2_s + \sigma^2_e}\\ \\ ITT_{tt^\prime} & = \frac{\sigma^2_b + tt^\prime \sigma^2_s}{\left[(\sigma^2_b + t^2 \sigma^2_s + \sigma^2_e)(\sigma^2_b + {t^\prime}^2 \sigma^2_s + \sigma^2_e)\right]^\frac{1}{2}} \end{aligned} \]

In this example, \(\sigma^2_s = 0.01\) (and the other variances remain as before), so

\[ ITT_{33} = \frac{0.15 + 3^2 \times 0.01}{0.15 + 3^2 \times 0.01 + 2} =0.1071\] and

\[ ITT_{36} = \frac{0.15 + 3 \times 6 \times 0.01}{\left[(0.15 + 3^2 \times 0.01 + 2)(0.15 + 6^2 \times 0.01 + 2)\right ]^\frac{1}{2}} =0.1392\]

Here’s the data generation:

defc <- defData(varname = "ceffect", formula = 0, variance = 0.15, 
                dist = "normal", id = "cluster")
defc <- defData(defc, "cteffect", formula = 0, variance = 0.01, 
                dist = "normal")

defa <- defDataAdd(varname = "Y", 
  formula = "0 + ceffect + 0.10  * period + cteffect * period + 1 * rx", 
  variance = 2, dist = "normal")

dc <- genData(100, defc)
dp <- addPeriods(dc, 7, "cluster")
dp <- trtStepWedge(dp, "cluster", nWaves = 4, lenWaves = 1, startPer = 2)
  
dd <- genCluster(dp, cLevelVar = "timeID", numIndsVar = 10, 
                 level1ID = "id")
dd <- addColumns(defa, dd)

And here is the comparison between observed and expected ICC’s. The estimates are quite variable, so there appears to be slight bias. However, if I generated more than 200 data sets, the mean would likely converge closer to the expected values.

In the next post (or two), I plan on providing some examples of fitting models to the data I’ve generated here. In some cases, fairly standard linear mixed effects models in R may be adequate, but in others, we may need to look elsewhere.

References: