Power for linear models of longitudinal data with applications to Alzheimer's Disease Phase II study design
In longpower: Sample Size Calculations for Longitudinal Data

library(knitr)
opts_chunk$set(echo = TRUE, message = FALSE, warning = FALSE, tidy = FALSE, comment = NA)

\newcommand{\var}{\mathrm{var}} \newcommand{\corr}{\mathrm{corr}} \newcommand{\cov}{\mathrm{cov}} \newcommand{\x}{\mathbf{x}} \newcommand{\y}{\mathbf{y}} \newcommand{\ttt}{\mathbf{t}} \newcommand{\uu}{\mathbf{u}} \newcommand{\vv}{\mathbf{v}} \newcommand{\Ipl}{I_{\psi\lambda}} \newcommand{\Ill}{I_{\lambda\lambda}}

Introduction

We will discuss power and sample size estimation for randomized placebo controlled studies in which the primary inference is based on the interaction of treatment and time in a linear mixed effects model [@laird]. We will demonstrate how the sample size formulas of @liu for marginal or model fit by generalized estimating equation (GEE) [@zeger] can be adapted for mixed effects models. Finally, using mixed effects model estimates based on data from the Alzheimer's Disease Neuroimaging Initiative (ADNI), we will give examples of sample size calculations for models with and without baseline covariates which may help explain heterogeneity in cognitive decline and improve power.

Power calculations

Exchangeable correlation and random intercept models

Suppose we wish to estimate the required sample size for inference regarding the interaction of treatment and time in a longitudinal, placebo controlled study. Such calculations are relatively straightforward when the inference is based on a GEE model in which the correlation structure is assumed to be "exchangeable." An exchangeable correlation structure specifies that all observations from within the same cluster, or repeated measures on the same subject, are equally correlated. This is exactly equivalent to a random effects model which includes a random intercept for each cluster of correlated observations. Sample sizes for study designs using these models can be calculated using a simple formula such as that in [@diggle], page 29. The formula requires the number visits, the interval between visits, the estimated model variance ($\sigma^2$), the within subject correlation ($\rho$), and of course the usual sample size calculation inputs (power, significance level, and effect size).

To translate the formula of \cite{diggle} to the random effects setting, let us first consider the details of the assumed error structure of the GEE framework. The GEE model assumes that the response for subject $i$ at time $t_{ij}$, denoted $Y_{ij}$, is the group mean, dependent on time and treatment, plus an error term $\varepsilon_{ij}$. Or, borrowing notation from \cite{diggle}, for group A: [ Y_{ij} = \beta_{0A} + \beta_{1A}t_{ij} + \varepsilon_{ij},\quad i=1,\ldots,m; j=1,\ldots,n. ] and similarly for Group B. The null hypothesis is $H_0: d=\beta_{1A} - \beta_{1B} = 0$. Under an exchangeable correlation structure $\var(Y_{ij})=\var(\varepsilon_{ij})=\sigma^2$ and $\corr(Y_{ij},Y_{ik})=\corr(\varepsilon_{ij}, \varepsilon_{ik})=\rho$, for all subjects, $i$, and time points $j, k$.

In the mixed effects framework we can assume a random intercept model which is equivalent to the GEE model with exchangeable correlation structure. In this case we believe $\varepsilon_{ij} = \alpha_i + \varepsilon_{ij}^$, where $\alpha_i$ is the random intercept term shared by all observations and $\varepsilon_{ij}^$ are independent and identically distributed (iid) error terms. We see that $\var(Y_{ij})=\var(\varepsilon_{ij})=\var(\alpha_i) + \var(\varepsilon_{ij}^)$ and $\corr(Y_{ij},Y_{ik})=E[(\alpha_i + \varepsilon_{ij}^)(\alpha_i + \varepsilon_{ik}^)]/\sigma^2=\var(\alpha_i)/\sigma^2$. The variance of the random intercept, $\var(\alpha_i)$, and the residual variance, $\var(\varepsilon_{ij})$, are easily obtainable from the output of mixed effects fitting software so that one might fit a random effects model to pilot data to educate a power calculation using the GEE formula of \cite{diggle}. Assuming equal numbers in the placebo and active groups, a common visit schedule for all subjects ($t_{ij} = t_{kj}$ for all $i,j,k$), and a random intercept model; the number of subjects per group is: [ m = \frac{2(z_\alpha + z_Q)^2(\var(\alpha_i) + \var(\varepsilon_{ij}^))^2(1-\var(\alpha_i)/\sigma^2)}{ns_x^2d^2} ] where $z_p$ is the $p$th standard normal quantile, $Q$ is $1-P$, $P$ is the specified power, and $s_x^2=n^{-1}\sum_j(t_{j}-\bar x)^2$.

General correlation and random slope models

The random intercept model is not equipped to handle variations in the rate of change from subject to subject. In many diseases, such as Alzheimer's disease, the rate of improvement or decline will vary greatly within the treatment group, regardless of treatment. This variation can be modeled with a random slope term. That is, we assume: [ Y_{ij} = \beta_{0A} + \beta_{1A}t_{ij} + \alpha_{0i} + \alpha_{1i}t_{ij} + \varepsilon_{ij}^, ] where we use $\varepsilon_{ij}^$ again to denote iid error and reserve $\varepsilon_{ij}$ for possibly correlated error. If we derive the correlation structure of $\varepsilon_{ij}=\alpha_{0i} + \alpha_{1i}t_{ij} + \varepsilon_{ij}^$, which is necessary in order to use GEE-based sample size formulas, we find that we no longer have an exchangeable correlation structure. In fact $\var(Y_{ij})=\var(\varepsilon_{ij})=\var(\alpha_{0i})+t_{ij}^2\var(\alpha_{1i}) + 2t_{ij}\cov(\alpha_{0i},\alpha_{1i}) +\var(\varepsilon_{ij}^)$ and $\cov(Y_{ij},Y_{ik})=\cov(\varepsilon_{ij},\varepsilon_{ik})=\var(\alpha_{0i})+t_{ij}t_{ik}\var(\alpha_{1i}) + (t_{ij}+t_{ik})\cov(\alpha_{0i},\alpha_{1i})$. For the common visit schedule case, the covariance matrix for the vector of correlated errors, $\mathbf{\varepsilon_i}=(\varepsilon_{i1},\ldots,\varepsilon_{in})'$, is of the form: [ \Sigma = [(\var(\alpha_{0})+t_{j}t_{k}\var(\alpha_{1}) + (t_{j}+t_{k})\cov(\alpha_{0},\alpha_{1}))]_{jk}+{\rm diag}(\var(\varepsilon^*_j)) ] With this specification of the covariance matrix, one can use the sample size formula of \cite{liu} for linear GEE models (page 941). (Warning: The formula given on the bottom page 29 of \cite{diggle} for general correlation matrices, $R$, is wrong).

The formula for linear models provided by \cite{liu} is useful for testing $H_0: \mathbf{\psi=0}$ for any linear model of the form: [ Y_{ij} = \x_{ij}'\mathbf{\psi} + \mathbf{z}{ij}'\mathbf{\lambda} + \varepsilon{ij} ] where $\mathbf{\varepsilon_i}\sim N(\mathbf{0},\sigma^2R)$ and the covariates for individual $i$, $\x_{i}=(\x_{i1}', \ldots, \x_{i1}')'{n\times p}$ and $\mathbf{z}{i}=(\mathbf{z}{i1}', \ldots, \mathbf{z}{i1}')'{n\times q}$, arise from a known discrete distribution. For our placebo controlled longitudinal study, the fully specified model is of the form: [ Y{ij}=\beta_0 + \beta_1{{\rm Group}{i}=A} + \beta_2t{ij} + \beta_3t_{ij}{{\rm Group}{i}=A}. ] That is, the parameter of interest for the interaction of treatment and time is $\psi = \beta_3$ and nuisance parameter is $\mathbf{\lambda} = (\beta_0,\beta_1,\beta_2)'$. The covariates are distributed as $\x{i} = \ttt = (t_1, \ldots, t_n)'$ and $\mathbf{z}j = [\mathbf{1}\, \mathbf{1}\, \ttt]{n\times3}$ with probability 1/2 (Group A); and $\x_{i} = \mathbf{0}$ and $\mathbf{z}j = [\mathbf{1}\, \mathbf{0}\, \ttt]{n\times3}$ with probability 1/2 (Group B).

The Liu and Liang's formula for linear models can be coded In R as:

library(longpower)

liu.liang.linear.power

The parameters include \texttt{d}, the effect size (possibly vector); \texttt{u}, the list of covariate vectors or matrices associated with the parameter of interest; \texttt{v}, the respective list of covariate vectors or matrices associated with the nuisance parameter; \texttt{sigma2}, the error variance; \texttt{R}, the correlation structure; and \texttt{Pi} the proportion of covariates of each type (\texttt{u}, \texttt{v}, and \texttt{Pi} are expected to be the same length and sorted with respect to each other).

For example, we can reproduce the table of sample sizes (per group) on page 29 of \cite{diggle} for the given exchangeable correlations with $\mathbf{t} = (0,2,5)'$, $\alpha=0.05$, power=0.80, and $d=0.5$ via the \texttt{\hlkwd{diggle.linear.power}} function:

n = 3 # visits
t = c(0,2,5)
rho = c(0.2, 0.5, 0.8)
sigma2 = c(100, 200, 300)
tab.diggle = outer(rho, sigma2, 
      Vectorize(function(rho, sigma2){
        ceiling(diggle.linear.power(
          d=0.5,
          t=t,
          sigma2=sigma2,
          R=rho,
          alternative="one.sided",
          power=0.80)$n[1])}))
colnames(tab.diggle) = paste("sigma2 =", sigma2)
rownames(tab.diggle) = paste("rho =", rho)
tab.diggle

or via the \texttt{\hlkwd{liu.liang.linear.power}} function:

u <- list(u1 = t, u2 = rep(0,n))
v <- list(v1 = cbind(1,1,t),
         v2 = cbind(1,0,t))         
tab.ll <- outer(rho, sigma2, 
      Vectorize(function(rho, sigma2){
        ceiling(liu.liang.linear.power(
          delta=0.5, u=u, v=v,
          sigma2=sigma2,
          R=rho, alternative="one.sided",
          power=0.80)$n[1])}))
colnames(tab.ll) <- paste("sigma2 =", sigma2)
rownames(tab.ll) <- paste("rho =", rho)
tab.ll

As a second example, consider an Alzheimer's disease trial in which assessments are taken every three months for 18 months (7 visits). We assume an smallest detectable effect size of 1.5 points on the cognitive portion of the Alzheimer's Disease Assessment Scale (ADAS-Cog). This is a 70 point scale with great variability among sick individuals. We assume the random intercept to have a variance of 55, the random slope to have a variance of 24, and a residual variance of 10. The correlation between random slope term and random intercept term is 0.8. We can estimate the necessary sample size by first generating the correlation structure. Since $\varepsilon = \var(Y_{ij})$ is not constant over time in this model, we fix \texttt{sigma2=1} and set \texttt{R} equal to the covariance matrix for $\mathbf{\varepsilon_i}$:

# var of random intercept
sig2.i <- 55
# var of random slope
sig2.s <- 24
# residual var
sig2.e <- 10
# covariance of slope and intercep
cov.s.i <- 0.8*sqrt(sig2.i)*sqrt(sig2.s)

cov.t <- function(t1, t2, sig2.i, sig2.s, cov.s.i){
        sig2.i + t1*t2*sig2.s + (t1+t2)*cov.s.i 
}

t <- seq(0,1.5,0.25)
n <- length(t)
R <- outer(t, t, function(x,y){cov.t(x,y, sig2.i, sig2.s, cov.s.i)})
R <- R + diag(sig2.e, n, n)
u <- list(u1 = t, u2 = rep(0,n))
v <- list(v1 = cbind(1,1,t),
         v2 = cbind(1,0,t))         

liu.liang.linear.power(d=1.5, u=u, v=v, R=R, sig.level=0.05, power=0.80)

So the study would require about 207 subjects per arm to achieve 80\% power, with a two-tailed $\alpha=0.05$.

The simple formula provided in \cite{diggle} suggests the required number of subjects can be found by $2(z_\alpha+z_Q)^2\xi/d^2$, where [ \xi_\textrm{WRONG}= \left(\begin{array}{cc} 0 & 1\ \end{array}\right) \left(\begin{array}{ccc} 1 & \ldots & 1 \ t_1 & \ldots & t_n \ \end{array}\right)R^{-1} \left(\begin{array}{ccc} 1 & t_1 \ \vdots & \vdots\ 1 & t_n \end{array}\right) \left(\begin{array}{c} 0 \ 1 \end{array}\right). ] Executing this for our Alzheimer's example, we get a sample size of:

x <- (rbind(1,t)%*%solve(R)%*%cbind(1,t))[2,2]
x*2*(qnorm(1-0.05/2) + qnorm(0.80))^2/1.5^2

which is clearly wrong. In fact, there is a typo in \cite{diggle}. The correct formula for $\xi$ is: \begin{equation}\label{eq:diggle3} \xi = \left(\begin{array}{cc} 0 & 1\ \end{array}\right) \left[ \left(\begin{array}{ccc} 1 & \cdots & 1 \ t_1 & \cdots & t_2\end{array}\right) (\sigma^2R)^{-1} \left(\begin{array}{cc} 1 & t_1 \ \vdots & \vdots \ 1 & t_m\end{array}\right)\right]^{-1} \left(\begin{array}{c} 0 \ 1 \end{array}\right). \end{equation}

Applying the correct formula, we get

x <- solve(rbind(1,t)%*%solve(R)%*%cbind(1,t))[2,2]
x*2*(qnorm(1-0.05/2) + qnorm(0.80))^2/1.5^2

Similarly, using \cite{liu}, we attempt to derive the correct closed form formula for this specific linear model. The required sample size per group is given as [ m = \nu/(\psi_1'\tilde\Sigma_1\psi_1) ] where [ \tilde\Sigma_1 =\sigma^{-2}\sum_{l=1}^m\pi_l (\uu_l'-\Ipl\Ill^-1\vv_l')R^{-1} (\uu_l'-\vv_l\Ill^-1\Ipl'), ] [ \Ipl=\sigma^{-2}\sum_{i=1}^m\pi_l \uu_l'R^{-1}\vv_l, ] and [ \Ill=\sigma^{-2}\sum_{i=1}^m\pi_l \vv_l'R^{-1}\vv_l. ] Again, in our case the probability of each of the two covariate values is $\pi_1=\pi_2=1/2$; and $\uu_1 = (t_1, \ldots, t_n)'$, $\vv_1 = [\mathbf{1}\, \mathbf{0}\, \x_i]{n\times3}$, $\uu_2 = \mathbf{0}$, and $\vv_2 = [\mathbf{1}\, \mathbf{0}\, \x_i]{n\times3}$. We have

\begin{eqnarray} \Ipl & = & \sigma^{-2}/2\uu_1'R^{-1}\vv_1\ \Ill & = & \sigma^{-2}/2[\vv_1'R^{-1}\vv_1 + \vv_2'R^{-1}\vv_2]=1/2X]\ \Ipl\Ill^{-1} & = & \uu_1'R^{-1}\vv_1X^{-1}\ \Ill^{-1}\Ipl' & = & X^{-1}\vv_1'R^{-1}\uu_1 \end{eqnarray}

\begin{eqnarray} \tilde\Sigma_1 & = & \sigma^{-2}/2 [(\uu_1-\uu_1'R^{-1}\vv_1X^{-1}\vv_1') R^{-1} (\uu_1-\vv_1X^{-1}\vv_1'R^{-1}\uu_1)\ & & + \uu_1'R^{-1}\vv_1X^{-1}\vv_2'R^{-1}\vv_2X^{-1}\vv_1R^{-1}\uu_1\ & = & \sigma^{-2}/2[\uu_1R^{-1}\uu - \uu_1'R^{-1}\vv_1X^{-1}\vv_1'R^{-1}\uu_1] \end{eqnarray} Applying this to our working example:

X <- t(v[[1]])%*%solve(R)%*%v[[1]] + 
    t(v[[2]])%*%solve(R)%*%v[[2]]

Sigma1 <- ((t(u[[1]])%*%solve(R)%*%t - 
           t(u[[1]])%*%solve(R)%*%v[[1]]%*%solve(X)%*%t(v[[1]])%*%solve(R)%*%t)/2)

(qnorm(1-0.05/2) + qnorm(0.80))^2/(Sigma1*(1.5)^2)/2