A scientific review on advances in statistical methods for crossover design

Cox (1958) analyzed the crossover design, assuming that there is no carryover effect, and pointed out that it is nothing but a split-plot design. Grizzle (1965) pioneered testing for treatment effect in the presence of carryover assuming model (3) and calculated the analysis of variance (ANOVA) table. A correct table version can be found in the appendix of Hills & Armitage (1979). The author concluded that if there is a significant carryover, it is better to use a randomized parallel-group trial than a crossover design. But, when the presence of residual is in doubt, he proposed a two-stage (TS) procedure that first tests for carryover effect based on $T_{\lambda}$ in (5) and then tests for treatment effect using the statistic $T_{\tau,W}$ if the test for carryover is not significant, otherwise tests using $T_{\tau,B}$. It is worth noting that the preliminary carryover test lacks power as the standard error of the statistic $T_{\lambda}$ is influenced by between-subject variation. Hence, the author recommends testing at a higher significance level ($10\%$). Under the normality assumption, the test-statistics for testing $H_{0}:\tau=0$ is given by

	$\displaystyle F_{co}=nT_{\tau,w}^{2}/\hat{\sigma}^{2}\sim F_{1,2n-2}\left(\phi\left(\left(\tau-\lambda/2\right)^{2}\right)\right)$		(6)
	$\displaystyle F_{par}=nT_{\tau,b}^{2}/2\left[\hat{\sigma}^{2}+\hat{\sigma}^{2}_{s}\right]\sim F_{1,2n-2}\left(\phi\left(\tau^{2}\right)\right)$		(7)
	$\displaystyle F_{TS}=F_{co}\;I\left(\lvert T_{\lambda}\rvert\leq c\right)+F_{par}\;I\left(\lvert T_{\lambda}\rvert>c\right)$		(8)

Where $c$ is an appropriate cutoff value and $\hat{\sigma},\hat{\sigma_{s}}$ be the maximum likelihood estimators of $\sigma,\sigma_{s}$ respectively. The subscripts of the F-statistics in (6)-(8) are meant to highlight that the tests are based on crossover (co) and parallel group (par) trials and two-stage procedure (TS). The non-centrality parameters vanish when $\tau-\lambda/2=0$ and $\tau=0$, respectively. Confidence intervals (CI) for the treatment effect can be obtained by inverting the above $F$ tests. This can be done using SAS procedure proc mixed (SAS Institute Inc, 2008) or glm() function in R (R Core Team, 2018)

Brown Jr (1980) analyzed the cost efficiency of the crossover design by incorporating the sample size that would be required to get the same precision and concluded that if there is insignificant carryover, crossover design is super economical than a parallel group design, but then it can be costlier otherwise.

Under the null hypothesis $\tau=0$, $F_{par}$ has a central F-distribution. If $\lambda\neq 0$, the non-centrality parameter in $F_{co}$ does not vanish, and the test becomes biased. Willan & Pater (1986) pointed out that if there is no carryover effect without the treatment effect, then under $H_{0}$ $F_{co}$ has a central F-distribution. Consequently, the test becomes unbiased. But, the problem occurs otherwise. They mentioned that for a placebo-controlled trial, or when a new treatment is compared against a standard treatment, it is expected not to have a carryover when there are no treatment effects. They also differentiated between physical and psychological carryover and asserted that significant carryover can exist without treatment effect for the former but not for the latter. The authors compared the power of the crossover and parallel group tests and derived

$\displaystyle P_{co}(\tau,\lambda)>P_{par}(\tau,\lambda)\quad\textrm{if}\quad\left\lvert\lambda/\tau\right\rvert<2-\sqrt{2(1-\rho)}$

(9)

where $P$ is the power function of the crossover and parallel group tests and $\rho=\sigma_{s}^{2}/(\sigma_{s}^{2}+\sigma^{2})$ is the intra-group correlation. This means the power of F-tests in (6) and (7) depends on the value of intra-group correlation and the relative magnitude of carryover and treatment effect. Given $\sigma_{s}^{2}>\sigma^{2}$, if the carryover effect is small relative to the direct treatment effect, the crossover test of $\tau$ can outperform the parallel group test. Later on, Willan (1988) proposed to use the maximum test-statistic $F_{max}=\max\;(F_{co},F_{par})$ for testing the treatment effect under the assumption that no differential carryover exists in the absence of treatment effects.

Because the test-statistic $T_{\lambda}$ and $T_{\tau,B}$ are highly correlated, the TS procedure can be biased even if $\lambda$ is small. Freeman (1989) examined the performance of $F_{co}$, $F_{par}$ and $F_{TS}$ for testing treatment effect but they allowed the possibility that $\lambda$ can be significant even if $\tau=0$. Simulation results showed that the nominal significance level of the TS procedure is higher than the stipulated type-I error ($\alpha$) even when $\lambda=0$; it goes far above the actual as $\lambda$ deviates from zero. In the absence of carryover, $F_{TS}$ and $F_{co}$ are more powerful than $F_{par}$, but $F_{TS}$ achieves superiority at the expense of increased significance level. Senn (1996) proposed a correction to the TS procedure to maintain the type-I error rate ($\alpha$). If the type-I error of the carryover test is $\alpha_{1}$, he proposed to conduct the $F_{par}$ test at a reduced significance level of $\alpha\alpha_{1}$. This procedure is conservative except when within-subject variability $\sigma^{2}=0$. Wang & Hung (1997) proposes another two-stage procedure that controls the type-error rate like the above and offers a detailed comparison between maximum test statistics and two-stage procedures. The conclusion is similar: if the relationship between carryover and treatment is unknown, then the two-stage procedures can be very misleading. But, under the assumption that carryover does not persist in the absence of treatment effect, $\lambda/\tau$ is negative and large in magnitude, two-stage procedures are more robust to the possible overestimation that maximum test-statistic or the crossover test can lead to.

Another disturbing feature of the TS procedure mentioned by (Senn (2002)) is if there is a significant effect of sequence (which is aliased with the carryover in AB/BA design), the estimator $T_{\tau,b}$ is not unbiased for $\tau$. The test of treatment effect based on only first-period observations will be biased.

The performances of all the existing frequentist approaches rely on the information on the significance of carryover. Selwyn et al. (1981) pioneered the Bayesian analysis for the crossover trials in the Bio-equivalence study by incorporating the uncertainty as an informative prior to the carryover effect. Arguing that such information is rarely available in practice, Grieve (1985) assumed a completely uninformative prior for the parameters. The posterior distribution of $\tau$ is centered at the difference of the first-period observation, and the second-period observation comes through the variance. This means if we do not assume that the carryover effect is insignificant, then we have to do our analysis based on first-period observations, a similar point made by Grizzle (1965). To incorporate a prior belief on the carryover effect, Grieve (1985) proposed a Bayes factor approach that quantifies the model’s probability with carryover and no carryover. This is more practical than Selwyn et al. (1981) approach because the variance $\sigma_{R}^{2}$ in the prior distribution of $\lambda$ is not available in practice. He mentioned that if we have prior knowledge about carryover, then we should have prior information on the treatment effect too, so why not incorporate that information in the prior of $\tau$? The conclusion is that the inference for the treatment effect is heavily based on our prior belief on the carryover effect, which all the frequentist approaches conclude.

Often, the experimenter does not have an adequate extent of certainty about the carryover effects. Therefore, a complete specification of the prior distribution becomes difficult. David & Seaman Jr (2002) introduced a hierarchical approach for the Bayesian analysis of two-period crossover design to ensure robustness to prior choice. Li & Sivaganesan (2017) proposed an objective Bayesian estimation approach for AB/BA crossover trials. The authors assume a mixture distribution of the model with carryover and no carryover similar to Grieve (1985), in addition to a Zellner–Snow prior distribution (Zellner & Siow, 1980) for $\lambda$ when carryover is present. The simulation study shows that the treatment effect’s nominal coverage probability is maintained close to $95\%$ level, irrespective of whether $\lambda$ is small or large. In terms of MSE, this method significantly improves the TS procedure when $\lambda$ is moderate.

Contemporary to Grizzle (1965)’s parametric inference, Koch (1972) proposed a non-parametric approach for analyzing $2\times 2\times 2$ crossover design using Wilcoxon rank sum test or equivalently by Mann-Whitney test. Let $R_{ik},\tilde{R}_{ik}$ be respectively the rank of subject totals $T_{ik}$ and period differences $P_{ik}$ for the kth subject in the ith group based on the entire sample. Let, $R_{i}=\sum_{k}R_{ik}$ and $\tilde{R}_{i}=\sum_{k}\tilde{R}_{ik},i=1,2$. The Wilcoxon rank-sum test-statistic for testing $H_{01}:\lambda=0$ and $H_{02}:\tau=0$ under $\lambda=0$, is given by,

$\displaystyle Z_{\lambda}=\frac{R_{1}-E(R_{1})}{\sqrt{\textrm{Var}(R_{1})}}\quad;\quad Z_{\tau}=\frac{\tilde{R}_{1}-E(\tilde{R}_{1})}{\sqrt{\textrm{Var}(\tilde{R}_{1})}}$

(10)

Where the expectation and variance of $R_{1}$ and $\tilde{R}_{1}$ under the null hypothesis can be obtained from Mann & Whitney (1947). If the sample size is large, under null, the test statistics follow an approximately standard normal distribution. The exact CI for $\tau$ can be obtained using the methods of Hodges Jr & Lehmann (1963) by computing the median of all possible response differences in the two groups. SAS procedure proc nparway and R function wilcox.test() offers CI, an exact and asymptotic significance test. It should be noted that the test of treatment effect is valid only if $\lambda=0$. Otherwise, the test should be done based on period-1 responses. A similar two-stage procedure can also be adopted here to test for the effects of the treatment. Koch (1972) also proposed non-parametric test for testing $H_{03}:\lambda=0,\tau=0$. Under $H_{03}$, the responses for each subject are an identically distributed bivariate random variable, and the Wilcoxon rank-sum test of location following Chatterjee & Sen (1964) can be constructed. Let $\gamma_{ijk}$ be the rank of the kth subject in group $i$ on period $j$ based on the entire sample. Define ${\bm{\mathbf{{\gamma}}}}_{ik}=(\gamma_{i1k},\gamma_{i2k})^{\top}$ as the rank vector and $\bar{{\bm{\mathbf{{\gamma}}}}}_{i}=(\bar{\gamma}_{i1},\bar{\gamma}_{i2})^{\top}$, where $\bar{\gamma}_{ij}=\sum_{k=1}^{n_{i}}\gamma_{ijk}/n_{i}\;i=1,2$ and $\bar{{\bm{\mathbf{{\gamma}}}}}=\left(n_{1}\bar{{\bm{\mathbf{{\gamma}}}}}_{1}+n_{2}\bar{{\bm{\mathbf{{\gamma}}}}}_{2}\right)/(n_{1}+n_{2})$ Then, the test statistic is given by

		$\displaystyle W_{\tau}=\left(\bar{{\bm{\mathbf{{\gamma}}}}}_{1}-\bar{{\bm{\mathbf{{\gamma}}}}}_{2}\right)^{\top}{\bm{\mathbf{{\Sigma}}}}_{\gamma}^{-1}\left(\bar{{\bm{\mathbf{{\gamma}}}}}_{1}-\bar{{\bm{\mathbf{{\gamma}}}}}_{2}\right)\sim\chi^{2}_{2}\quad(\textrm{under }H_{03}\textrm{ asymptotically})$		(11)
	where	$\displaystyle{\bm{\mathbf{{\Sigma}}}}_{\gamma}=\frac{(n_{1}+n_{2})}{n_{1}n_{2}(n_{1}+n_{2}-1)}\sum_{i=1}^{2}\sum_{k=1}^{n_{i}}\left({\bm{\mathbf{{\gamma}}}}_{ik}-\bar{{\bm{\mathbf{{\gamma}}}}}\right)\left({\bm{\mathbf{{\gamma}}}}_{ik}-\bar{{\bm{\mathbf{{\gamma}}}}}\right)^{\top}$		(12)

The non-parametric tests described in (11) are based on the assumption that the variations of the samples are the same for both groups. Cornell (1991) proposed two non-parametric tests for equality of the dispersion of the responses in two groups. One is based on the rank correlation coefficient of $T_{ik}$ and $V_{ik}$, and the other is based on their concordance. The test is based on the fact that under model (3), $\textrm{Cov}(T_{ik},P_{ik})=\sigma_{1}^{2}-\sigma_{2}^{2},i=1,2$, which vanishes under the null hypothesis.

Tudor & Koch (1994) proposed a more general non-parametric approach to test for significant total treatment (combined of direct and carryover) effect. To accomplish that, they tested $H_{0}:2\tau-\lambda=0$. The test statistics are based on the period differences of two groups $P_{ik},i=1,2$, which is given by

$\displaystyle Q=\frac{\left(\bar{P}_{1\bullet}-\bar{P}_{2\bullet}\right)^{2}}{S_{p}^{2}}\quad;\quad\textrm{where }S_{p}^{2}=\frac{(n_{1}+n_{2})}{n_{1}n_{2}(n_{1}+n_{2}-1)}\sum_{i=1}^{2}\sum_{k=1}^{n_{i}}(P_{ik}-\bar{P})^{2}$

(13)

Where $\bar{P}=(n_{1}\bar{P}_{1\bullet}+n_{2}\bar{P}_{2\bullet})/(n_{1}+n_{2})$. Under $H_{0}$, if $n_{i}$ are large ($>15$), $Q$ follows approximately $\chi^{2}$ distribution with $1$ degrees of freedom. The authors mentioned the above sample size requirement for each group to ensure that the total number of subjects in the study is large enough ($>30$) to hold asymptotic results. As they suggested, the above test has a higher power of detecting any significant total treatment effect because it is entirely based on within-subject differences. If the test is rejected, further tests must be done to detect whether the difference is due to direct treatment and/or carryover effects. In the second stage, a test for carryover based on $Z_{\lambda}$ needs to be done. If the test is rejected, a test for $\tau$ based on first-period responses can be done; otherwise, a test using $Z_{\tau}$ as in (10) can be carried out. Sometimes, knowledge of the length of the washout period can be incorporated to judge the absence of carryover effects. This involves a three-stage procedure. The authors recommended setting the significance levels for each test before the start of the test to ensure that the nominal significance level is not inflated.

Jung & Koch (1999) proposed a rank measure association-based non-parametric test under crossover design applicable to an ordinal response. Consider a randomized parallel group design with two treatments. Let $\mu_{i}$ be the mean response for the ith group. The difference $\mu_{1}-\mu_{2}$ represents a measure of disparity between the two groups. In the parametric approach, we model this disparity as $\mu_{1}-\mu_{2}=E(X_{1k}-X_{2k})=\tau_{A}-\tau_{B}$. For the Wilcoxon rank-sum test, we model the median of differences between the responses through $\tau_{A}-\tau_{B}$. Similarly, Jung proposed to model the disparity between the two groups $\psi=P(X_{1k}-X_{2k}>0)$ non-parametrically by $\tau=\tau_{A}-\tau_{B}$. In an AB-BA crossover design, we have $4$ cell means. So, there will be four different rank measures, ${\bm{\mathbf{{\psi}}}}=(\psi_{1},\psi_{2},\psi_{3},\psi_{4})^{\top}$.

If $\hat{{\bm{\mathbf{{\psi}}}}}$ be the estimate of ${\bm{\mathbf{{\psi}}}}=(\psi_{1},\psi_{2},\psi_{3},\psi_{4})^{\top}$ and ${\bm{\mathbf{{\beta}}}}=(\tau,\pi,\lambda)^{\top}$ be our parameter of interest, then, they proposed the model $\hat{{\bm{\mathbf{{f}}}}}=\textrm{logit}(\hat{{\bm{\mathbf{{\psi}}}}})={\bm{\mathbf{{X}}}}{\bm{\mathbf{{\beta}}}}$ where the design matrix ${\bm{\mathbf{{X}}}}$ can easily be obtained from the relation between ${\bm{\mathbf{{\psi}}}}$ and ${\bm{\mathbf{{\beta}}}}$. The estimates $\hat{{\bm{\mathbf{{\psi}}}}}_{m}$ is constructed non-parametrically using the U-statistics of the responses related in $m$th comparison $m=1\dots 4$, e.g., for the first comparison, we consider only the responses in the first period. The author also derived the asymptotic multivariate normal distribution of the estimator $\hat{{\bm{\mathbf{{\beta}}}}}$ with a suitably chosen covariance matrix, which can be used to test any hypothesis of the form $H_{0}:{\bm{\mathbf{{C}}}}{\bm{\mathbf{{\beta}}}}=0$ and to derive approximate $100(1-\alpha)\%$ CI. The performance of this method is not shown in the paper through simulation but is applied to real data, and the results were similar to the standard procedures. For an extensive application of the method, see Chapter $2$ of Jones & Kenward (2014).

For higher order designs with more than two periods, the model in (3) assumes that the carryover effect of a treatment in the period $j$ depends only on the treatment in period $j-1$. This is called first-order carryover. Moreover, the carryover effect of a treatment on itself (self-carryover) is assumed to be identical to the carryover of the treatment on another treatment (mixed-carryover), e.g. in the design $d(2,3,4)$ shown in table 1(a), although the second group receives treatment A and the group $4$ receives the same treatment (B) again at period $2$, carryover effect of B from the first period in group $2$ assumed to be the same as the carryover of B in group $4$. However, in many practical applications, this simple carryover model is inadequate. For example, in customer survey experiments for the evaluation of certain products, if a person strongly likes or dislikes a product, its lingering effect influences the evaluation of the products that follow it. If the product is intensely disliked, the subsequent product will likely get a below-average evaluation. Also, this lingering effect diminishes proportionally as the person evaluates other products. Similarly, in pharmaceutical trials, while comparing more than two drugs under a crossover design, the carryover effect of a drug in a certain period can be significantly different if the same drug is administered in the next period compared to the carryover effect if a different drug is applied in the next period. This is because the impact of a drug is significantly dependent on the combination in which the other drugs are administered. Senn (2002) provides a mathematical explanation to justify that this simple carryover model is implausible in the pharmacokinetic study (see Chapter $10$ for details).

These problems can be addressed by including an interaction between treatment and carryover effects (Kershner & Federer, 1981). However, including an interaction term is feasible for a two-treatment design (d.f. for the interaction is 1), but for a study with several treatments, many interaction terms will be present in the model. To tackle this difficulty, Afsarinejad & Hedayat (2002) introduced a model that differentiates self-carryover and mixed carryover effects. Kempton et al. (2001) introduced another model where the carryover effect of a treatment in the current period is proportional to the treatment effect in the previous period. Another approach is to include a second-order carryover term in the model. That means the response at period $j$ will also be affected by the treatment at period $j-2,j\geq 3$. An excellent repository of models with different carryover effects can be found in Ozan & Stufken (2010). We assume that all the parameters are estimable through within-subject differences. A general model is

$\displaystyle Y_{ijk}=\mu+\pi_{j}+\alpha_{i}+\beta_{d(i,j)}+S_{ik}+{\bm{\mathbf{{\gamma}}}}^{\top}{\bm{\mathbf{{B}}}}_{ik}+\epsilon_{ijk}$

(14)

where $\alpha_{i}$ is the effect of ith group, $\beta_{d(i,j)}$ is the effect of the treatment and carryover associated with $j$th period for $i$th group and ${\bm{\mathbf{{B}}}}_{ik}$ be the vector of baseline covariates measured for subject $k$ in group $i$. Any baseline covariates could also be included in the basic model (3), defined for AB-BA design. Depending on the design, one can include treatment by period interaction term and other factors in the model. Different forms of $\beta_{d(i,j)}$ are summarized in chapter 4 of Jones & Kenward (2014).

In a crossover study with $t$ treatments, we aim to efficiently estimate any pair of treatment contrast $\tau_{i}-\tau_{j}$. To accomplish that, we want to design the study in such a way that the SE of the estimator, $\textrm{Var}(\hat{\tau}_{i}-\hat{\tau}_{j})$ is minimized compared to that for any competing design. Kiefer (1975) introduced different optimality criteria based on the total variance of the estimator of all possible treatment contrasts. This is essentially a function of the trace of information matrix $C_{d}$ under the model assumed. Similarly, optimality criteria for an estimation of the carryover effect can also be defined. Hedayat et al. (1978) pioneered finding optimal designs by considering a subclass of design called uniform designs. In a uniform design, each treatment is assigned to the same number of subjects in each period, and for each subject, every treatment is applied in the same number of periods. Uniformity is appealing because the treatment and carryover effect becomes orthogonal to each other, and it implies $p=\lambda_{1}t$ and $g=\lambda_{2}t$, where $\lambda_{1},\lambda_{2}$ is a positive integer. However, assuming $p=t$ and $g=\lambda_{2}t$, they proved that balanced uniform (BU) designs are universally optimal for estimating treatment and carryover effects over the class of uniform designs. In a balanced design, each treatment is preceded by every other treatment an equal number of times, and a treatment can’t be preceded by itself. Balance makes the standard error of all pairs of treatment contrasts the same.

Cheng et al. (1980) introduced a new class of design called strongly balanced designs and proved that any strongly balanced uniform design is universally optimal for estimating direct and carryover effects over the entire class of designs $d(t,p,g)$. A strongly balanced design is balanced with the relaxation that treatment can precede itself. If a design is strongly balanced, the number of groups (g) must be a multiple of $t^{2}$ and $p/t=m>1$. An example of a two-treatment strongly balanced design is ABBA-BAAB-AABB-BBAA (Kershner & Federer, 1981). This means when $p>t$, these designs repeat the same treatment in consecutive periods. Kunert et al. (1984) worked with $p=t$ and claims that if there exists a BU design, then the ratio of the trace of the information matrices of BU design and any competing design has a lower bound of $1-(t^{2}-t-1)^{-2}$. This means that BU designs are nearly optimal for $t>2$ but not strictly optimal. Kunert also proved that for $t=2=p$, the optimal designs for estimating direct treatment and carryover effects are different, and none is uniform.

Sen & Mukerjee (1987) proved that when an interaction between treatment and carryover is present in the model, strongly balanced designs are still optimal for estimating direct treatment effects because they are orthogonal. Kempton et al. (2001) developed optimality results when carryover effects are proportional to treatment effects. Kunert & Stufken (2002) developed optimal designs in the presence of self and mixed carryover effects by introducing a new class of designs called totally balanced designs. When the objective is to compare a bunch of test treatments with a control treatment, Hedayat et al. (2005) obtained optimality results with $p=t+1$, called extra-period designs. Many works have been done to find optimal designs when $p<t$, we refer the reader to chapter 3 of Bose & Dey (2009). We conclude our discussion on optimal designs by making two points: First, one should prioritize the effects or the contrasts that are of primary concern before choosing a design. Second, in the possibility of subject drop-out, will the design still be nearly optimal if the experimenter has to stop the study before the last period? R-package Crossover (Rohmeyer, 2017) offers an excellent interface to choose between several designs under different carryover models by comparing the average efficiencies of the competing design under different correlation structures of error such as independence, compound symmetry (CS), etc.

Depending on the design and the nature of experiments, the analyst might choose any model from (14) for the treatment and carryover effects with some baseline covariates recorded at the start of the study. For now, we assume that a single response is observed within a period for every subject. As stated before, crossover design with $p$ periods can be considered split-plot experiments where the whole-plots are the subjects, and measurements are observed within each subject at different periods. Standard ANOVA $F$-test techniques for split-plots using random subject effect (whole plot error) can be conducted to test for treatment or carryover effects under a Gaussian assumption using proc mixed or R package lmerTest (Kuznetsova et al., 2017)

Alternatively, crossover design can also be viewed as a repeated measurements study where measurements are observed for each subject over consecutive periods. Defining, ${\bm{\mathbf{{Y}}}}_{ik}=(Y_{i1k},\cdots,Y_{ipk})^{\top}$ as the vector of response for the kth subject in sequence $i$, the random subject effect model in (3) induces an exchangeable or CS covariance structure for the response vector i.e. $\textrm{Var}({\bm{\mathbf{{Y}}}}_{ik})={\bm{\mathbf{{\Sigma}}}}=\sigma^{2}{\bm{\mathbf{{I}}}}_{p}+\sigma_{s}^{2}{\bm{\mathbf{{J}}}}_{p}\forall\;i,k$. However, this CS structure might be too parsimonious to capture the true covariance. To generalize this we specify the mean model as $E({\bm{\mathbf{{Y}}}}_{ik}\;|\;{\bm{\mathbf{{X}}}})={\bm{\mathbf{{X}}}}_{ik}{\bm{\mathbf{{\theta}}}}$. The fixed parameters ${\bm{\mathbf{{\theta}}}}$ contain all the fixed parameters specified in the mean model assumed, and the design matrix ${\bm{\mathbf{{X}}}}_{ik}$ is specified accordingly. Furthermore, the covariance model is $\textrm{Var}({\bm{\mathbf{{Y}}}}_{ik}\;|\;{\bm{\mathbf{{X}}}})={\bm{\mathbf{{V}}}}_{ik}(\xi)$, which need not be compound symmetric. The vectors ${\bm{\mathbf{{\theta}}}}$ and ${\bm{\mathbf{{\xi}}}}$ are called mean and covariance parameters. This method directly captures the effect of treatment or carryover on the population mean response, which is called a population-averaged (PA) model. The influence of any random subject effect is captured in the complex covariance matrix. We assume that the responses for different subjects are mutually uncorrelated.

Under Gaussian assumption, the mean parameter ${\bm{\mathbf{{\theta}}}}$ is estimated through the generalized least square (GLS) approach and covariance parameters can be estimated using the restricted maximum likelihood (REML) method to adjust for the bias in the estimation of covariance. If normality is not assumed, then the mean component is estimated using linear estimating equations. SAS procedure proc mixed, proc glimmix and R-package nlme (Pinheiro et al., 2018) offers flexible PA models. The covariance structure can be varied according to any baseline covariate, e.g., the subject’s sex. Many different covariance matrix structures can be assumed, such as Autoregressive (AR), CS, Ante, and Unstructured (Un), the most general one. The final choice of the covariance model is made by AIC or BIC criterion. If the model is correctly specified, then GLS estimator of ${\bm{\mathbf{{\theta}}}}$ with its asymptotic distribution is given by

$\displaystyle\hat{{\bm{\mathbf{{\theta}}}}}=\left(\sum_{i,k}{\bm{\mathbf{{X}}}}_{ik}^{\top}{\bm{\mathbf{{V}}}}_{ik}^{-1}(\hat{\xi}){\bm{\mathbf{{X}}}}_{ik}\right)^{-1}\sum_{i,k}{\bm{\mathbf{{X}}}}_{ik}^{\top}{\bm{\mathbf{{V}}}}_{ik}^{-1}(\hat{\xi}){\bm{\mathbf{{Y}}}}_{ik}\sim\textrm{Asymptotically N}\left({\bm{\mathbf{{\theta}}}},\hat{{\bm{\mathbf{{\Sigma}}}}}_{M}\right)$

(15)

where $\hat{\xi}$ is the REML/ML estimator of $\xi$ and $\hat{{\bm{\mathbf{{\Sigma}}}}_{M}}=\left(\sum_{i,k}{\bm{\mathbf{{X}}}}_{ik}^{\top}{\bm{\mathbf{{V}}}}_{ik}^{-1}(\hat{\xi}){\bm{\mathbf{{X}}}}_{ik}\right)^{-1}$ is called the model-based SE. If the sample size in each group is not large, adjustment due to Kenward & Roger (1997) can be done. However, this adjustment still understates the SE of the estimator if the covariance structure is incorrectly specified. Diggle et al. (2002) pointed out in chapter 4 that the assumption of diagonal covariance does not affect the estimator of the mean parameter ${\bm{\mathbf{{\theta}}}}$. Still, it overstates or understates the SE of the estimator. The robust or empirical SE introduced in the context of non-normal response using a generalized estimating equation (GEE) gives protection against wrongly specified covariance (Zeger et al., 1988). The authors proved that under certain assumptions, as long as the mean model is correctly specified, the GLS estimator of the mean parameter still follows an asymptotic normal distribution with the covariance matrix replaced by $\hat{{\bm{\mathbf{{\Sigma}}}}}_{R}$ (Yuan & Jennrich, 1998), which is greater than the model-based SE (in the positive definite sense). The empirical option in SAS proc mixed allows inference on mean parameters using robust standard error. Adjustment to the empirical SE due to Mancl & DeRouen (2001) can be done for a small sample. Wald-type test for testing $H_{0}:{\bm{\mathbf{{C}}}}{\bm{\mathbf{{\theta}}}}=\phi$ and approximate CI can be constructed using the test-statistic

$\displaystyle T=\left({\bm{\mathbf{{C}}}}\hat{{\bm{\mathbf{{\theta}}}}}-\phi\right)^{\top}\left({\bm{\mathbf{{C}}}}\hat{{\bm{\mathbf{{\Sigma}}}}}{\bm{\mathbf{{C}}}}^{\top}\right)^{-1}\left({\bm{\mathbf{{C}}}}\hat{{\bm{\mathbf{{\theta}}}}}-\phi\right)\sim\chi^{2}_{r}\quad\textrm{asymptotically}\;\;;\quad\;\;r=\textrm{rank}({\bm{\mathbf{{C}}}})$

(16)

If the primary questions of interest are to investigate the effect of treatment on each subject, then a subject-specific (SS) linear mixed effect model is appropriate. If ${\bm{\mathbf{{b}}}}$ is the vector of random effects then a SS model can be written as $E({\bm{\mathbf{{Y}}}}_{ik}\;|\;{\bm{\mathbf{{X}}}},{\bm{\mathbf{{b}}}})={\bm{\mathbf{{X}}}}_{ik}{\bm{\mathbf{{\theta}}}}+{\bm{\mathbf{{Z}}}}_{ik}{\bm{\mathbf{{b}}}}$. Usually, ${\bm{\mathbf{{b}}}}$ is assumed to follow Gaussian with mean $0$ and unknown covariance ${\bm{\mathbf{{D}}}}$. The PA model that generates from the SS model by taking an expectation over ${\bm{\mathbf{{b}}}}$ is $E({\bm{\mathbf{{Y}}}}_{ik}\;|\;{\bm{\mathbf{{X}}}})={\bm{\mathbf{{X}}}}_{ik}{\bm{\mathbf{{\theta}}}}$, $\textrm{Var}({\bm{\mathbf{{Y}}}}_{ik}\;|\;{\bm{\mathbf{{X}}}})={\bm{\mathbf{{Z}}}}_{ik}{\bm{\mathbf{{D}}}}{\bm{\mathbf{{Z}}}}_{ik}^{\top}+{\bm{\mathbf{{R}}}}_{ik}$, where ${\bm{\mathbf{{R}}}}_{ik}$ is the covariance of the measurement error, possibly diagonal. This generates an equivalent PA model. However, in pharmacokinetic studies or when the carryover effect is assumed to be proportional to the treatment effect, the model becomes non-linear (Kempton et al., 2001). An SS model can be defined as $E({\bm{\mathbf{{Y}}}}_{ik}\;|\;{\bm{\mathbf{{X}}}},{\bm{\mathbf{{b}}}})={\bm{\mathbf{{f}}}}({\bm{\mathbf{{Z}}}}_{ik},{\bm{\mathbf{{\theta}}}},{\bm{\mathbf{{b}}}})$ and an appropriate SS covariance model with parameter ${\bm{\mathbf{{\psi}}}}$. A PA model can be derived from this by integrating the likelihood over the distribution of ${\bm{\mathbf{{b}}}}$. Numerical approximations such as the adaptive quadrature method, a first-order approximation of likelihood are typically used to approximate the marginal distribution, and the estimation of ${\bm{\mathbf{{\theta}}}}$ and ${\bm{\mathbf{{\xi}}}}$ is carried out through the GEE method or maximizing the approximated likelihood directly. See Davidian & Giltinan (2003) for a series of methodologies available in this context. SAS procedure proc nlmixed implements this method with great efficiency.