An optimal multi-period crossover design for an application in paediatric nephrology

Children undergoing regular haemodialysis for renal failure need to attend the dialysis unit (DU) two or three times per week. On each visit the dialyser is connected to the patient’s circulation. To achieve this the patient usually has a surgically inserted venous central line which remains in place permanently and is accessed on each visit to the DU. During the interdialytic period the blood in the lumen of the line may clot. Unclotting a line is a time-consuming problem and if it happens regularly then the central line will have to be re-sited. Re-siting is both unpleasant for the patient and undesirable, as it is a procedure that cannot be repeated indefinitely because of the limited number of sites which can be used for this purpose. For these patients failure to achieve reliable venous access has dire consequences. Consequently, a small amount of anticoagulant (a “lock”) is injected into the lumen of the central line at the end of each dialysis session in order to prevent clotting. The usual anticoagulant is heparin (H). This article is concerned with design issues that arose in a trial to compare heparin with alteplase (A), an alternative anticoagulant.

The dialysis treatment of choice for children with renal failure is peritoneal dialysis, but there are a few patients for whom this is unsuitable, or for whom peritoneal dialysis is no longer viable. Consequently, in a paediatric population there are only a few children attending for haemodialysis at any centre: fewer than ten patients in a UK regional centre at any given time is typical. Moreover, differences in treatment protocols between centres mean that a multicentre study would be difficult to run. However, each patient who does attend the DU is obliged to do so three (or, for occasional patients, two) times a week. Therefore information on the comparison between A and H can best be obtained using a crossover design, as this allows each patient to provide numerous observations, and at the same time improve the precision of the estimated comparison between A and H by eliminating inter-patient differences. The trial will measure the weight of clot in the blood withdrawn from the central line at the start of each dialysis session. As this withdrawal of blood is part of the usual clinical procedure to flush the line prior to dialysis, it imposes no burden on the patient and as such can repeated many times.

Numerous crossover designs are available in the literature (Jones and Kenward, 2003, chapters 4 & 5) but designs with more than ten periods are unusual in most areas of application and almost unheard of in clinical settings. Moreover, what optimality properties these designs possess are usually based on one of a handful of rather general models which may not be the most appropriate choice for the current trial. The next section considers the technical issues which arise in developing a design for this study, in particular the choice of a realistic statistical model for this application and its implications for design selection. In Section 3 an optimal design is derived for this model and in Section 4 some practical issues are discussed.

To avoid causing practical difficulties on the DU it was decided to run the trial for $w$ weeks, so for most patients a $3w$-period design is required but for those patients who need to attend only twice per week, a $2w$-period design is needed. This is the first non-standard feature of the design problem: almost all “off-the-shelf” designs assume that each experimental subject is observed equally often.

The second issue is that most available designs assume that the response observed on subject $i$ in period $j$, $y_{ij}$, follows the model:

$$y_{ij}=\mu+\xi_{i}+\pi_{j}+\tau_{d}(i,j)+\epsilon_{ij},$$

(1)

where $1\leq i\leq n$ and $1\leq j\leq p$. Here $\mu$ is a general mean; $\xi_{i}$ is some form of subject effect, either fixed or random; $\pi_{j}$ is the effect of period $j$; $\tau_{t}$ is the effect of treatment $t$ and $d(i,j)$ is the treatment allocated to subject $i$ in period $j$, and $\epsilon_{ij}$ are independent residuals with a common variance. More specialised designs specific to the crossover setting have been derived for a model which explicitly acknowledges the possibility of the effect of a treatment carrying over to the following period, namely

$$y_{ij}=\mu+\xi_{i}+\pi_{j}+\tau_{d}(i,j)+\gamma_{d(i,j-1)}+\epsilon_{ij}.$$

(2)

where $\gamma_{t}$ is the carryover effect of $t$ and we set $\gamma_{d}(i,0)=0$: see (Jones and Kenward, 2003, chapters 4 and 5) for details.

For our purposes the former model (1) is more relevant because the effects of the anticoagulants will not carry over to subsequent periods. This is because only a small quantity of anticoagulant is used at the end of each dialysis session. The volume is just sufficient to fill the lumen of the central line and it largely stays in situ until the following session. Moreover, at the start of the next session the line is aspirated, thereby removing most of the anticoagulant; any small amount left in the line will be flushed out in the course of the dialysis. Systemic effects due to small quantities of anticoagulant diffusing from the lumen into the circulation will be negligible because i) the volume of anticoagulant involved will be very small compared with the volume of the circulation and ii) the agents will be metabolised long before the subsequent dialysis session.

Optimal designs for model (1) are row-column designs such as Latin square and Youden designs, in which each treatment occurs equally often on each subject and in each period (Shah and Sinha, 1989, Chapter 4). These conditions are rather restrictive in that they impose conditions on the divisibility of $p$ and $n$ by the number of treatments. Moreover, they arise from the very general form of period effect allowed by (1), where a separate parameter is specified for each period.

In the present application the condition that each treatment, i.e. A and H, should be applied equally often to each patient is unlikely to be troublesome. However, having A and H occur equally often in each period is potentially more troublesome to arrange. The model in (1) goes back at least to animal feeding trials (Cochran et al., 1941) where the appropriateness of the period effect would have been much clearer. Its widespread adoption in the theory of crossover designs has often had a precautionary flavour - analysts feeling that designs needed to be protected again the possible untoward influences of non-specific period effects (Matthews, 1994a, b), rather than that there was any clear-cut reason for the form adopted in (1).

In fact, the study of aspects of dialysis is an area where the inclusion of period effects in the statistical model can be justified. As pointed out by Matthews and Hoenich (1994), many variables which can be measured during or following dialysis will be affected by the point in the dialysis cycle when the observation is made. Most DUs operate a cycle of dialysis sessions whereby a patient attends on Mondays, Wednesdays and Fridays (or on Tuesdays, Thursdays and Saturdays but for definiteness we will assume the former, as the alternatives are equivalent). On Mondays, variables such as the size of clots, or the amount of urea available for removal from the blood will have accumulated over the three days since attending the previous Friday, so can be expected to be different from the corresponding quantity on a Wednesday or Friday, as this will have accumulated over only two days.

We therefore decided to adopt the following model for the expected weight of clot aspirated from a patient $i$ who attends the DU $\ell$ times per week. In period $j=1,\ldots,\ell w$.

$$\mathsf{E}(y_{ij})=\tau d(i,j)+\pi_{\ell}(i,j)+\xi_{i},$$

(3)

where $d(i,j)=1$ if H is allocated to patient $i$ in period $j$ and -1 if A is allocated. For patients attending three times per week

	$\displaystyle\pi_{3}(i,j)$	$\displaystyle=$	$\displaystyle\pi_{1},\text{if }j\text{ is a Monday}$
		$\displaystyle=$	$\displaystyle\pi_{2},\text{if }j\text{ is a Wednesday}$
		$\displaystyle=$	$\displaystyle\pi_{3},\text{if }j\text{ is a Friday.}$

Note that the observation $y_{ij}$ is the weight of clot aspirated following allocation of treatment $d(i,j)$ but this is only obtained when the lumen is aspirated at the following dialysis session. Therefore, for example, the observation $y_{i,\text{Friday}}$ which relates to the treatment allocated on the Friday, is not obtained until the following Monday. Given that the interval between Wednesday and Friday is the same as that between Monday and Wednesday, a possible simplification of (2) is to set $\pi_{1}=\pi_{2}$. We do not adopt this slightly more prescriptive form as it transpires that this extra restriction provides no important practical advantage in the specification of the design.

For those patients required to attend only on Mondays and Fridays we replace $\pi_{3}(i,j)$ by $\pi_{2}(i,j)$ to accommodate the different attendance schedule, namely

	$\displaystyle\pi_{2}(i,j)$	$\displaystyle=$	$\displaystyle\pi_{4},\text{if }j\text{ is a Monday}$
		$\displaystyle=$	$\displaystyle\pi_{3},\text{if }j\text{ is a Friday.}$

The parameter $\pi_{3}$ is common between $\pi_{2}(i,j)$ and $\pi_{3}(i,j)$ as both groups of patients attend on Monday having attended on the preceding Friday.

The final aspect of the model which needs consideration is the residual term $\epsilon_{ij}=y_{ij}-\mathsf{E}(y_{ij})$. This represents the within-patient variation and these terms will be assumed to be independent with common variance $\sigma^{2}$. The assumption of independence is questionable when each individual is observed serially over time: we will briefly consider the possible effect of dependence on the design calculations in the next section. It is also possible that the assumption of homoscedasticity will prove unfounded, although we would aim to address such issues by suitable transformation of the $y_{ij}$.

We suppose that $N_{3}$ patients dialysed three times per week and that $N_{2}$ patients dialysed twice weekly will be recruited, providing a total of $m=3wN_{3}+2wN_{2}$ observations. The aim is to determine the best treatment sequences to use to define the $d(i,j)$. A good review of this area can be found in Stufken (1996). There are a variety of optimality criteria which can be used but because we are comparing only two treatments these all coincide and optimal designs are chosen to minimise the variance of the estimator of $\tau$.

Although planning for the study assumed four patients treated three times per week and two twice-weekly, by the time the trial started, seven patients treated thrice-weekly, i.e. $N_{3}=7$ and two treated twice weekly, $N_{2}=2$, were available. Full details of the study are reported in Gittins et al. (2007).

There were some difficulties in the conduct of the trial. The most important of these was that although the allocation to the treatment for the final period was done correctly, confusion on the DU meant that the corresponding observation of clot weight was missed for all thrice weekly patients. That this observation was to have been made at the start of August, when the junior medical staff change over, may have been relevant. In addition there are some isolated discrepancies: patient 8, one of the twice-weekly patients, missed their final four observations (three H, one A) due the availability of a kidney transplant. Patients 1, 2 and 7 fell short of the revised total of 29 observations: patient 2 missed three observations due to holidays, while patients 1 and 7 missed one each due to errors on the DU. Although any deviation from the planned trial is undesirable, given the complexity of the design this level is probably no more than should be expected. While the omission of a few observations will reduce the power of the study, the reasons for the missing observations should not induce any noticeable bias. Moreover, the study was planned on the basis of $N_{3}=4,N_{2}=2$, but in fact $N_{3}=7,N_{2}=2$ were recruited, so the power would in fact have been higher than planned.

Fitting the model (3) gave an estimate of the mean semi-treatment difference, $\hat{\tau}$=2.79g, P=0.13, 95% confidence interval (-0.79,6.37), indicating a higher mean weight of clot on heparin compared with alteplase. A normal probability plot of residuals shown in Figure 1(a) indicates that the model assumptions were questionable. A randomization analysis based on 2000 re-randomizations of the treatment allocation using the scheme outlined in Section 3.3 gave P=0.123, confirming the model-based P-value.

The probability plot in Figure 1(a) indicates that the data are skewed. This is shown, separately for each treatment, in Figure 2 where it is apparent that in a high proportion of occasions (30% H, 50% A) no clot was retrieved from the central line. A natural response is to skewed data is to take logs of the observations but with so many zeroes we need to use a transformation of the form $\log(y_{ij}+k)$. Inspection of the data showed that the range of non-zero clot weights retrieved was 0.6g to 146g (H) and 0.7g to 114g (A). Rather than use a formal approach to estimating $k$, we chose $k=0.1g$. More formal approaches using a shifted Box-Cox transformation suggested a log transformation with $k\approx 0.001g$. The estimate of the treatment effect ($2\tau$, on log scale) becomes 0.830, $P=0.013$, 95% (0.175,1.485) ($k=0.1$) and 1.440, $P=0.018$, 95% (0.250,2.631) ($k=0.001$). This sensitivity to $k$, while typical of data with many zeroes, is clearly undesirable.

An alternative approach is to take the analysis in two stages. First the probability of a clot-free sample is assessed, and then the mean weights when clot is present are analysed. In both cases a generalized estimating equation (GEE) was used, with patient used to define the clustering variable. The first analysis gave an odds ratio of clotting as 2.34 times larger on H than A, $P<0.001$, 95% confidence interval (1.24,4.42). The second analysis addressed the considerable remaining skewness in the non-zero clot weights by fitting a GEE with gamma variances and a log link. This indicated that mean clots weights were 1.86 times heavier using H than A ($P<0.001$), 95% confidence interval (1.45,2.39).

This application required a crossover design with many more periods than is usually the case. While multi-period designs are available from families of designs already in the literature, these require balance across patients: for example by ensuring equal replication of treatments at each period. They also required that each patient be treated an equal number of times. In the present application the second of these requirements would have needed the trial to be extended for the twice-weekly patients which, while possible, would have been very inconvenient for the DU. Practical difficulties could arise with the requirement to balance over periods. Moreover, this would arise from the use of a model for period effects that is unnecessarily general. The use of the model (3) tailored to the application means that an optimal design can be constructed by ensuring that each patient follows a particular form of sequence, with no requirement to balance across patients. This is a helpful feature when implementing the design.

The optimality calculations made no distributional assumptions although it did make second-order assumptions, namely that the residual terms were uncorrelated and had constant variance. Even if these assumptions are untrue, randomization tests can provide a valid test of the hypothesis of no treatment effect. In the event, the data proved to be rather different from that anticipated, with many instances where no clot was formed, especially when alteplase was the used as the anticoagulant. A model-based analysis assuming Normal errors gave a $P$-value very similar to that obtained from a randomization test. However, the conclusion of the analysis was substantially different if $\log(y_{ij}+k)$, rather that the clot weights $y_{ij}$ were analysed. This suggests that there was also substantial heteroscedasticity in the data. Moreover, the analyses were sensitive to the value of $k$, a parameter poorly estimated by the data. Analyses of data of this form will seldom be adequately summarised by a single paramter such as a mean or median, so something similar to the two-stage method used will usually be necessary.

A more difficult question to address is the extent to which the properties of the design depart from optimal now that it is clear that the original assumptions are untenable. Optimal design theory for crossover trials started with the model (2) and has developed by considering departures from this model in terms of the form of the carryover treatment term and the dependence structure of the residuals. There is no work on design theory for crossover trials with other forms of outcome, such as a binary variable. Designs that are optimal for the mixed binary and continuous outcome that we have pursued here would be extremely challenging to determine.

Moreover, the real issue is not how far the optimal design for the known outcomes differs from that used at the outset but could we have chosen a design that would have had good, but not necessarily optimal properties for a range of assumptions about the outcome. There are many papers going back over many years, which derive designs which have good properties over a range of assumptions (Cook and Nachtsheim, 1982; DuMouchel and Jones, 1994), including recent contributions which address the problems posed by generalized linear models (Woods et al., 2006; Dror and Steinberg, 2006). However, most of these contributions address robustness to mis-specification of the linear predictor or of the link function, rather than the distribution of the the outcome variable.

Although derived using optimal design theory, the present design might have been arrived at by more general considerations of orthogonality. The design used in this study is optimal because it ensures the treatment effect is orthogonal to the period and patient effects. Such a property can be shown to have formal advantages when the outcomes are uncorrelated with constance variance. However, the notion that each treatment occurs equally often on each patient and on each treatment day, has more general appeal. However, rigorous confirmation of the value of this design for the type of data that actually arose in this trial is elusive.

I am grateful to Drs Malcolm Coulthard and Nicola Gittins for allowing me to use their data in this paper.

The matrices $B_{2}$ and $B_{1}$ are

$$B_{2}=\begin{pmatrix}I_{N_{3}}\otimes 1_{3w}&0(2wN_{3},N_{2})\\ 0(2wN_{2},N_{3})&I_{N_{2}}\otimes 1_{2w}\end{pmatrix}\qquad=\begin{pmatrix}I_{N_{3}}\otimes 1_{3w}&0\\ 0&I_{N_{2}}\otimes 1_{2w}\end{pmatrix}$$

where $I_{n}$ is the $n\times n$ identity matrix and $1_{n}$ is an $n\times 1$ vector of ones, and $0(a,b)$ denotes an $a\times b$ matrix of zeroes. Also

$$B_{1}=\begin{pmatrix}P\\ Q\end{pmatrix}=\begin{pmatrix}1_{wN_{3}}\otimes I_{3}&0(3wN_{3},1)\\ 0(2wN_{2},3)&I_{wN_{2}}\otimes K\end{pmatrix}$$

where

$$K=\begin{pmatrix}0&1\\ 1&0\end{pmatrix}$$

Therefore $B^{T}_{1}B_{1}=P^{T}P+Q^{T}Q$ can be written

$$\begin{pmatrix}wN_{3}\otimes I_{3}&0\\ 0&0\end{pmatrix}+\begin{pmatrix}0&0\\ 0&wN_{2}\otimes K^{T}K\end{pmatrix}=w\begin{pmatrix}N_{3}&&&\\ &N_{3}&&\\ &&N_{3}+N_{2}&\\ &&&N_{2}\end{pmatrix}$$

(10)

noting that $K^{T}K=I_{2}$.

The product $B_{1}^{T}B_{2}$ is

$$\begin{pmatrix}P^{T}&Q^{T}\end{pmatrix}\begin{pmatrix}I_{N_{3}}\otimes 1_{3w}&0\\ 0&I_{N_{2}}\otimes 1_{2w}\end{pmatrix}=\begin{pmatrix}P^{T}(I_{N_{3}}\otimes 1_{3w})&Q^{T}(I_{N_{2}}\otimes 1_{2w})\end{pmatrix}$$

Evaluating this matrix requires calculation of products such as $(1^{T}_{wN_{3}}\otimes I_{3})(I_{N_{3}}\otimes 1_{3w})$ and this is made easier if we note that, e.g., $1_{3w}=1_{w}\otimes 1_{3}$. We then obtain

$$B_{1}^{T}B_{2}=w\begin{pmatrix}1^{T}_{N_{3}}\otimes\begin{pmatrix}1_{3}\\ 0\end{pmatrix}&1^{T}_{N_{2}}\otimes\begin{pmatrix}0\\ 0\\ 1_{2}\end{pmatrix}\end{pmatrix}$$

To evaluate $A^{T}B_{1}$ let $A^{T}=(A^{T}_{3}\;\;A^{T}_{2})$ where $A_{\ell}$ is the part of $A$ relating to patients undergoing dialysis $\ell$ times per week. Then $A^{T}B_{1}=A^{T}_{3}P+A^{T}_{2}Q$ and we then obtain

$$A^{T}_{3}P=\begin{pmatrix}q^{3}_{M}&q^{3}_{W}&q^{3}_{F}&0\end{pmatrix}$$

and

$$A^{T}_{2}Q=\begin{pmatrix}0&0&q^{2}_{F}&q^{2}_{M}\end{pmatrix}.$$

Here

$$q^{3}_{M}=q^{3}_{M}(H)-q^{3}_{M}(A)$$

where $q^{3}_{M}(H)$ is the number of times H is allocated on a Monday to a thrice-weekly patient, and $q^{3}_{M}(A)$ is the corresponding count for alteplase: $q^{3}_{W}(\cdot)$ and $q^{3}_{F}(\cdot)$ refer to Wednesdays and Fridays, while $q^{2}_{X}(Y)$ are the corresponding quantities for patients who are dialysed twice-weekly. It follows that

$$A^{T}B_{1}=\begin{pmatrix}q^{3}_{M}&q^{3}_{W}&q^{3}_{F}+q^{2}_{F}&q^{2}_{M}\end{pmatrix}=q^{T}\text{, say}.$$

(11)

If we write $R=w^{-1}\text{diag}\begin{pmatrix}N_{3}^{-1}&N_{3}^{-1}&(N_{3}+N_{2})^{-1}&N_{2}^{-1}\end{pmatrix}$, then it follows that $A^{T}\mathscr{P}(B_{1})B_{2}$ can be expressed as $q^{T}R(B^{T}_{1}B_{2})$. Each element of $A^{T}\mathscr{P}(B_{1})B_{2}$ corresponds to a patient and each element is either $q^{T}RP_{3}$ for thrice-weekly patients or $q^{T}RP_{2}$ for twice weekly patients, where

$$P_{3}=\begin{pmatrix}1\\ 1\\ 1\\ 0\end{pmatrix}\qquad P_{2}=\begin{pmatrix}0\\ 0\\ 1\\ 1\end{pmatrix}$$