A correlation structure for the analysis of Gaussian and non-Gaussian responses in crossover experimental designs with repeated measures

Experimental designs are a very useful tool to analyze the effects of treatments applied to a set of experimental units (Hinkelmann and Kempthorne, 2005). In this kind of studies, it is frequent that experimental units are observed at a unique time point, this is known as a transversal experiment; notwithstanding, sometimes experimental units are observed several times during the study, keeping the same treatment, giving rise to longitudinal studies (Davis, 2002). There are situations in which an experimental unit receives all treatments, each one in a different period, this induces a setup where a sequence of treatments is applied to each unit. This kind of designs is known as crossover design (Jones and Kenward, 2015). In the scope of crossover designs, published results focus on the case of a normally distributed (Jones and Kenward, 2015) or binary (two possible outputs, namely success or failure) response variable. The latter case has been treated by means of generalized linear models for binary data (Ratkowsky et al. (1992), Curtin (2017) and Li et al. (2018)).
Jones and Kenward (2015, pag 204) described a crossover experiment with three treatments to control arterial pressure: treatment A is a placebo, treatments B and C are 20 and 40 mg doses of a test drug. Thus, there were six three-period sequences: ABC, ACB, BCA, BAC, CAB, and CBA, each one of them was applied to two individuals. In each period, 10 consecutive measurements of diastolic arterial pressure were taken: 30 and 15 minutes before, and 15, 30, 45, 60, 75, 90, 120 and 240 minutes after the administration of the treatment, as shown in Table 1.

. A second experiment carried out by researchers of the Department of Animal Sciences of Universidad Nacional de Colombia, was focused on inferring the effects of two diets, A (grass) and B (a mixture of grass and paper recycling waste) on dairy cattle performance. Eight cows split into two groups of four received the diets in such a way that the first group was fed diet A from day 1 to day 42, and diet B from day 43 to day 84, while the second group was fed diet B during the first period and diet A during the second one. Measurements of milk yield and quality, live weight and body condition score were taken at 1, 14, 28, 42, 56, 70 and 84 days (Jaime, 2019); the design structure is shown in Table2.

To analyze this sort of designs, Basu and Santra (2010), Josephy et al. (2015), Hao et al. (2015), Lui (2015), Rosenkranz (2015), Grayling et al. (2018), Madeyski and Kitchenham (2018) and Kitchenham et al. (2018) used mixed models for crossover designs with Gaussian response and a single observation per period including additive carryover effects. Biabani et al. (2018) presented a review on crossover designs in neuroscience, all papers considered Gaussian responses and did not account for carryover effects due to the presence of a washout period, even when it lasted a few days. On the other hand, Oh et al. (2003), Curtin (2017) and Li et al. (2018) used generalized linear models for crossover designs of two periods and two sequences of two treatments with a continuous (normal or gamma) or binary response variable and each experimental unit observed once per period. They used generalized estimating equations (GEE) to estimate model parameters. A Bayesian formulation of generalized linear models where the response was the survival time with a single observation per period was presented in Shkedy et al. (2005) and Liu and Li (2016).

Moreover, Dubois et al. (2011), Diaz et al. (2013) and Forbes et al. (2015) used Gaussian mixed models to analyze records from crossover designs with repeated measures using the area under curve as a strategy to obtain a single observation per period and one experimental unit, and they did not account for carryover effects.

In all aforementioned approaches, it is assumed that the crossover design considers a washout period between treatments and that the response variable of each individual is observed once per period. These assumptions are not fulfilled in the experiments described above because of two reasons: i) in patients with arterial hypertension the treatment cannot be stopped, while in the case of dairy cattle, cows must be fed every day; moreover, in both studies, the placebo is part of the treatment design, so considering it as a washout period would modify the experiment, ii) in each treatment, experimental units were observed several times per period. Therefore, there is a necessity for developing a consistent methodology to analyze this sort of experiments. In this study, we develop a methodology to analyze data from crossover designs with repeated measures using GEE and considering two correlation structures, between and within periods, which are combined via a Kronecker product to yield the overall correlation matrix. These approach was found to improve the estimation of parameters of interest with respect to methods based on GEE that consider a single correlation structure for all the observations of a subject; in addition, this method accounts for carryover effects, hence, it does not need a washout period.

In the second section, we present some background and define the methodology, in the third section we propose a method to estimate the correlation matrix, and provide theoretical results that support our estimation and modelling approaches. In the fourth section, we present a simulation study showing some advantages of our model as compared to the model with a single correlation structure for each period. Lastly, in the fifth section, the methodology is applied to real data from the aforementioned arterial pressure and dairy cattle experiments.

A crossover design has the following components (Patterson, 1951): i) Sequences which are each one of the distinct treatment combinations to be sequentially and randomly applied to the experimental units, ii) Treatments, which are randomly applied to each experimental unit within each sequence, iii) Periods, the time intervals in which the treatments that make up a sequence are applied, usually, each period is the same for all sequences and, consequently, the number of periods is equal to the number of elements of each sequence, iv) Experimental unit, the subjects or elements that receive the treatments, in each sequence there are $n_{l}$ experimental units, so the total number of experimental units in the study is denoted by $n=\sum_{l=1}^{S}n_{l}$.

Another important feature of the crossover design is the existence of carryover effects, defined in Vegas et al. (2016) as follows: persistency of the effect of a treatment over those that will be applied later, i.e., the treatment is applied in a given period, but its effect still affects the response in later periods when other treatments are applied, this residual effect is known as carryover effect. When the effect persists after one period, it is known as first order carryover effect, and when it persists after two periods, it is called a second order effect and so on (Patterson, 1951).

In a crossover design with $S$ sequences of length (the number of elements comprising each sequence) $P$, $Y_{ijk}$ is defined as the $k$th respose of the $i$th experimental unit at $j$th period. Let $n_{ij}$ be the number of observations of the $i$th experimental unit during the $j$th period. Then, vector $\boldsymbol{Y}_{ij}$ is defined as:

Also, vector $\boldsymbol{Y}_{i}$ is defined as:

which has dimension $\sum_{j=1}^{P}n_{ij}$, where $P$ is the number of periods.
According to these definitions, the ideas presented when discussing the arterial pressure and dairy cattle experiments, and assuming that $Y_{ijk}$ has a distribution in the exponential family, we propose the following model based on GEE:

where $g(\cdot)$ is the link function related to the exponential family, $\boldsymbol{x}_{ijk}$ is the vector of the design matrix corresponding to $k$th response from the $i$th experimental unit at the $j$th period, $\boldsymbol{\beta}$ is the vectors of fixed effects, $\mu$ is the overall mean, $\alpha_{i}$ is the effect of the $i$th sequence, $\gamma_{j}$ is the effect of the $j$th period, $\tau_{d[i,j]}$ is the effect of treatment $d$ applied in the period $j$ to $i$th experimental unit ($d=1,\ldots,q$), $\theta_{(1)[i,j-1]}$ is the first order carryover effect of $d$ treatment, $\theta_{d[i,j-u]}$ is the carry over effect of order $u$ of $d$ treatment. The carryover effects are considered in the two experiments discussed above because there was not a washout period. Due to the fact that observations of the same experimental unit are correlated, parameter estimation is carried out using GEE (Liang and Zeger, 1986). To this end, the following system of $q=dim(\boldsymbol{\beta})$ equations has to be solved:

where $\boldsymbol{\mu}_{i}=(\mu_{i11},\ldots,\mu_{iPn_{iP}})^{T}$, $\boldsymbol{y}_{i}=(y_{i11},\ldots,y_{iPn_{iP}})^{T}$, $\boldsymbol{x}_{mi}$ is the $m$th column of the design matrix of the ith experimental unit, and the covariance component $\boldsymbol{V}(\boldsymbol{\mu}_{i})$ is defined as:

where $\boldsymbol{D}(\cdot)$ is a diagonal matrix, $V(\mu_{ijk})$ is the variance function corresponding to the exponential family, and $R(\boldsymbol{\alpha})$ is the correlation matrix related to the covariance matrix $\boldsymbol{\Sigma}_{i}=Var(\boldsymbol{Y}_{i})$ as follows:

Due to the repeated measures structure in the crossover design, we propose a correlation structure of the form:

where $\boldsymbol{R}_{1}(\boldsymbol{\alpha}_{1})$ is the within-period correlation matrix, $\boldsymbol{\Psi}$ is the between-period correlation matrix, and $\otimes$ represents the Kronecker product (Harville, 1997). To estimate this matrix, we propose the following modified GEE to estimate $\boldsymbol{\beta}$ is:

where $\boldsymbol{V}(\boldsymbol{\mu}_{i})$ is defined as in Equation (5), $\boldsymbol{\mu}_{i}=\{\mu_{i11},\ldots,\mu_{iPL}\}$, $\boldsymbol{x}_{mi}$ is the $m$th column of the design matrix of the $i$th experimental unit.
The estimating equation for $\boldsymbol{\alpha}_{1}$ is:

where $\boldsymbol{H}_{i}=\boldsymbol{D}(V(r_{ijk}))_{q\times q}$ is a diagonal matrix, $\boldsymbol{\varepsilon}_{i}=E(\boldsymbol{W}_{i})_{q\times 1}$,
$\boldsymbol{W}_{i}=(r_{i11}r_{i12},$ $r_{i11}r_{i13},\ldots,r_{iP(L-1)}r_{iPL})^{T}_{q\times 1}$, and $r_{ijk}$ is the $ijk$th observed Pearson residual defined as:

and $q={P\choose 2}$. On the other hand, we propose the following estimator for $\boldsymbol{\Psi}=\left(\psi_{jj^{\prime}}\right)_{P\times P}$:

where $\boldsymbol{r}_{(j)i}$ is the vector of Pearson residuals for the $i$th experminetal unit at the $j$th period, and $\bar{\boldsymbol{r}}_{(j)}$ is the average (over i), $j\neq j^{\prime}$, with $j=1,\ldots,P$

The form of $\boldsymbol{R}_{1}(\boldsymbol{\alpha}_{1})$ in (6) is given by correlation structures such as: independence, autoregressive, exchangeable, etc (Davis, 2002).

The correlation structures are compared via the quasi-likelihood information criterion ($QIC$, Pan (2001)), which is defined as:

where $\hat{\mu}_{ijk}=\hat{\eta}_{ijk}=g^{-1}(\boldsymbol{x}_{ijk}\hat{\beta})$ is the estimated expected value of observation $Y_{ijk}$ under the model assuming the correlation matrix $R$, $\hat{\Omega}_{I}$ is the estimated covariance matrix of $\boldsymbol{\beta}$ under independence, and $\hat{V}_{R}$ is the covariance matrix of $\boldsymbol{\beta}$ under the model with correlation matrix $\boldsymbol{R}$ defined as in (6).

A simulation study was carried out to evaluate the proposed model. The setup was a crossover model with three periods $(P=3)$, five repetitions $(L=5)$ and two sequences $(S=2)$ within each period. The number of replicates per sequence $(n)$ varied from 2 to 100, and for each replicate, the simulation was ran 100 times. The following parameters are used in the simulation:

Four GEE models were fit per replicate, all had the same linear predictor and differed in the correlation structure: 1) independence, 2) first order autoregressive, 3) exchangeable, and 4) the Kronecker matrix proposed in this paper. Both, the simulations and fitted model were carried out using the libraries gee (Carey., 2019) and geeM (McDaniel et al., 2013) of R (R Core Team, 2022)

Figure 1 shows the proportion of times that each model had the smallest QIC relative to the 100 simulations for each number of replications per sequence (the mean and 95% confidence interval are shown). Moreover, an analysis of parameter estimation under each one of these models using the simulated data was performed; Figure 2 presents 95% confidence intervals for the location parameters under each model. In the simulation study it is observed that most of the times the lowest $QIC$ is obtained for the model with the true correlation structure. This suggets that the proposed estimation method along with theorem 3.1, which will allow to detect kronecker correlation structures in the data of the crossover design. In addition, this behavior is maintained even when the number of replicas within the design is low.
Regarding the coverage of confidence intervals for each parameter, in Figure 2 it can be noticed that the model with independence structure has a very low coverage for period and sequence effects; this will not allow a correct inference about the parameters of interest in the design. The other three models had more appropriate coverages, but the model with the within-period xchangeable matrix (that is, the true structure) exhibited values closer to 95%. As to the time effects, corresponding coverages were close to 95% in all four models; however, a subcoverage was observed for the model with a within-period exchangeable structure and an overcoverage for the independence model. The true model showed coverages corresponding to 95%.

The systolic blood pressure data is presented in Figure 3. The model used to analyze this experiment had a linear predictor considering fixed effects of treatment, period, baseline (the two measurements taken before applying the treatment), first and second order linear and quadratic carryover effects as a function of time. On the other hand, the model use to analyze data from the second experiment had a linear predictor considering fixed effects of baseline, treatment, period, and first order linear carryover effects as a function of time (quadratic effects were not considered because there were three observations within each period). The structures of the working correlation matrices were: i) independence, ii) first order autoregressive, and iii) fitted individually and were compared via the $QIC$ to determine the one exhibiting the best fit to each dataset. Table 3 presents the $QIC$ for each correlation structure in both datasets.

According to 3 the $\boldsymbol{\Psi}\otimes Exch_{10}$ and $\boldsymbol{\Psi}\otimes Exch_{3}$ correlation matrices had the smallest QIC values for the arterial pressure and dairy cattle experiments, respectively. Hence, correlation matrices having the Kronecker structure had the best fit in the two datasets. The correlation matrices estimated using equations (7) and (9) for the arterial pressure data are, respectively:

and for the dairy cattle dataset:

Figure 4 shows the matrix computed as the Kronecker product of matrices in equations (12) and (13). Notice the positive correlation between periods 1 and 3, a small positive correlation between periods 1 and 2. On the other hand, the matrix obtained as the Kronecker product of matrices in equations (14) and (15) is shown in Figure 5 where a positive but small correlation between periods 1 and 2 can be seen.

The matrices with the Kronecker structure are used to estimate the location parameters of the linear model (since these had the smallest QIC) yielding Table 4 4 for the arterial pressure data and Table 5 for the dairy cattle data. Each table shows the estimates, their standard errors (computed using the “sandwich” variance estimator (Hardin and Hilbe, 2003), the z statistic and the p-value corresponding to the null hypothesis $\beta=0$

Table 4 shows significant effects of baseline, treatment C (with respect to A), time (linear and quadratic), period, and carryover effect of treatment C. The interaction between carryover effects and time was not significant, which means that the effect of the previous treatment remains the same during the rest of the experiment. Moreover, Table 5 shows significant effects of baseline (milk yield before the beginning of the experiment), linear (negative) and quadratic (positive) regression coefficients of time, which means that milk yield shows a convex behavior during the experiment. There was no significant treatment effect, but there was a significant carryover effect of the paper-based diet over the grass-based diet, i.e., the cows being fed paper took a time to recover after changing to the grass-based diet.

Defining the correlation matrix structure is a highly relevant decision when using GEE, a proper specification of correlation structures in longitudinal data analysis improves estimation efficiency, leading to more reliable statistical inferences (Hin and Wang, 2009) . Hence, in this paper we develop a family of correlation structures that combine some of the classical structures used in longitudinal data analysis with an additional matrix that allows more flexibility in the overall correlation structure by adding a few parameters. Moreover, we provide an explicit estimation method of these parameters which features some sound statistical properties.

The theoretical results supporting some asymptotic properties of the proposed estimators are illustrated through the simulation study where a gain in goodness of fit was observed across the different simulation scenarios. The QIC was able to select the correct model, which had the proposed correlation structure. In addition, the confidence intervals built from the GEE showed a coverage that matched the nominal value. The estimation by intervals for each of the parameters presents coverage close to 95%, which shows a correct theoretical specification of each univariate interval.

As to the real data analysis, the results for the arterial pressure data showed the importance of accounting for carryover effects as they are useful for correctly estimating and interpreting the treatment affects across the time. If not included in the model, these residual effects may induce confusion problems. On the other hand, in the dairy cattle data, significant carryover effects were detected as well. Regarding the estimates of correlation matrices in both datasets, the QIC selected the proposed correlation matrix. Thus, the theoretical and empirical results from real and simulated data analyses suggest that the proposed methodology is promising and may be applied to perform better inferences from data obtained under crossover designs without washout periods and a repeated measures structure.