Analysing a 2x2 Crossover: Separating Treatment From Period and Carry-Over

Clinical Epidemiology ResearchMethodology and Research DesignUniqcret doctor knowledges
Analysing a 2x2 Crossover: Separating Treatment From Period and Carry-Over
On this page

อ่านฉบับภาษาไทย (Thai version)

Abstract

A 2x2 crossover trial gives every participant both treatments, one per period, in order AB or BA. Readings move with treatment, with a period effect shared by everyone in period 2, and possibly with carry-over: the first drug still acting after the switch. Half the difference between the two sequences' mean period differences (each participant's period-1 reading minus period-2 reading) estimates the treatment effect free of any period effect. Carry-over is inseparable from the sequence comparison. In a simulated trial of 24 adults with hypertension, the estimated effect of drug A versus drug B on systolic blood pressure is -8.02 mmHg (95% CI -11.74 to -4.30). With the simulated carry-over of -2 mmHg the estimator targets -7 mmHg; this draw lies near the true -8 by chance. The carry-over test has a power (chance of detecting a real effect) of 0.054. A washout, a treatment-free gap long enough for the first drug to fade, is the defence against carry-over, not a preliminary test.


Visual summary. Simulated data.

Two orders, two periods, one drug comparison

Twenty-four adults with hypertension join a trial of two fictional blood pressure drugs, A and B. Each takes both drugs in turn, with systolic blood pressure measured at the end of each treatment period. Twelve are randomised to the order AB and twelve to the order BA.

At the steering meeting, the interim summary shows the two orders moving in opposite directions after the switch (simulated data). Mean systolic blood pressure rises by 3.88 mmHg from period 1 to period 2 in the group that took A first (sequence AB) and falls by 12.15 mmHg in the group that took B first (sequence BA). Averaged over both orders, period 2 sits lower.

Part of each change belongs to the drug, part may belong to period 2 itself, and part to whatever the first drug left behind. The committee wants one number for the effect of A versus B. The analysis must separate that drug effect from the period and from any effect of the first drug that lingers after the switch.

What a crossover trial is, and what moves its readings

A crossover trial gives each participant every treatment under study, one at a time, in a randomised order. Each treatment interval is a period, and the order a participant follows is that participant's sequence. The 2x2 design has two periods and two sequences, AB and BA.

The design rests on a within-participant comparison: each person is measured on both drugs, so stable differences between people, such as their usual blood pressure, cancel. When those differences are large, a crossover can match the precision of a much larger parallel-group trial, in which each participant receives only one treatment [1, 2]. The price is that time and treatment order now enter the data.

Four effects can be named. The treatment effect $\Delta$ (delta) is the difference between drug A and drug B within a participant. The period effect $\pi$ (pi) is a shift shared by everyone in period 2, whichever drug they take.

The carry-over effect $\lambda$ (lambda) is an effect of the period-1 drug that persists into period 2. A washout period, a treatment-free interval between periods long enough for the first drug's effect to fade, is the design defence against carry-over.

A sequence effect is any systematic difference between the AB and BA groups as wholes. The table shows what each effect follows, and the next sections show why the last two cannot be separated.

What each effect follows in a 2x2 crossover.
EffectIt followsSymbolCan it be estimated?
Treatmentthe drug taken in the current period$\Delta$Yes, within participants (shifted by minus half the carry-over, if any)
Periodwhether the reading is in period 1 or period 2$\pi$Yes, within participants (shifted by half the carry-over, if any)
Carry-overthe drug taken in the previous period$\lambda$Only between participants, and only imprecisely
Sequencethe order, AB or BAnone of its ownNo: the design gives it the same pattern in the data as carry-over (the two are aliased), so no analysis can separate them

A model for the four cell means

Write $Y_{ij}$ for the outcome of participant $i$ in period $j$, where $j$ is 1 or 2. The design has four cells, one for each sequence in each period, and a cell mean is the average reading in one cell. The standard model for a 2x2 crossover is

$$Y_{ij} = \mu + s_i + \pi\,[\text{period 2}] + \Delta\,[\text{on A}] + \lambda\,[\text{period 2 of AB}] + e_{ij}$$

Each bracket equals 1 when its condition holds and 0 otherwise. The constant $\mu$ (mu) is the mean on drug B in period 1, $s_i$ is participant $i$'s own level relative to it, and $e_{ij}$ is the noise of a single reading. Treating $s_i$ as a draw from a normal distribution makes it a random intercept, the term a mixed model (a regression with a random level for each participant) uses for each participant's own level.

Strictly, $\lambda$ is the carry-over of A minus the carry-over of B. Carry-over common to both drugs looks exactly like a period effect, so the model carries only the difference, placed on period 2 of sequence AB.

Expected cell means under the model. Drug B is the reference, and the participant levels $s_i$ average out.
SequencePeriod 1Period 2
AB$\mu + \Delta$ (on A)$\mu + \pi + \lambda$ (on B)
BA$\mu$ (on B)$\mu + \pi + \Delta$ (on A)

Period differences: the estimator

For each participant, subtract the period-2 reading from the period-1 reading: $d_i = Y_{i1} - Y_{i2}$. The participant's own level $s_i$ appears in both readings and cancels. Write $\bar d_{AB}$ and $\bar d_{BA}$ for the mean of $d_i$ in each sequence.

From the cell means, $\bar d_{AB}$ estimates $\Delta - \pi - \lambda$ and $\bar d_{BA}$ estimates $-\Delta - \pi$. Their difference removes $\pi$, and their sum removes $\Delta$:

$$\hat\Delta = \tfrac{1}{2}\left(\bar d_{AB} - \bar d_{BA}\right)$$

The treatment estimate $\hat\Delta$ (delta-hat; a hat marks an estimate) is half the difference between the two sequences' mean period differences.

$$\hat\pi = -\tfrac{1}{2}\left(\bar d_{AB} + \bar d_{BA}\right)$$

The period estimate $\hat\pi$ is minus half their sum, the shift from period 1 to period 2. With no carry-over ($\lambda = 0$), $\hat\Delta$ estimates $\Delta$ and $\hat\pi$ estimates $\pi$.

The standard error (SE) of $\hat\Delta$ is half the standard error of a two-sample t-test comparing $d_i$ between the sequences (Hills and Armitage [3]). The SE is the spread of the estimate over many repetitions of the same trial. That t-test estimates $\bar d_{AB} - \bar d_{BA}$, which is twice $\hat\Delta$, so its t statistic and P value apply to $\hat\Delta$ unchanged.

Hand example: a period effect with no drug difference

This hand example uses a symptom score from 0 to 10, where higher is worse. In sequence AB the mean score is 8 on A in period 1 and 5 on B in period 2. In sequence BA it is 8 on B in period 1 and 5 on A in period 2.

  1. Mean period difference in AB

    \[ \bar d_{AB} = 8 - 5 = 3 \]

    Period 1 minus period 2 in sequence AB.

  2. Mean period difference in BA

    \[ \bar d_{BA} = 8 - 5 = 3 \]

    The same subtraction in sequence BA.

  3. Treatment effect

    \[ \hat\Delta = \tfrac{1}{2}(3 - 3) = 0 \]

    Drug A and drug B do not differ.

  4. Period effect

    \[ \hat\pi = -\tfrac{1}{2}(3 + 3) = -3 \]

    Scores are 3 points lower in period 2, whichever drug is taken.

Result: The whole fall from 8 to 5 belongs to period 2, and none of it to either drug.

Inside sequence AB alone, the fall coincides with the switch to B and would be credited to B. Inside sequence BA alone, the same fall would be credited to A. Only the two sequences together separate the drug from the period.

Carry-over is aliased with sequence

Now let the first drug linger. Put the cell means into the two estimators and take expectations, written $E(\cdot)$ for the average over many repetitions of the same trial:

$$E(\hat\Delta) = \Delta - \tfrac{\lambda}{2}$$

So the treatment estimate is shifted by minus half the carry-over. The period estimate is shifted the other way, $E(\hat\pi) = \pi + \lambda/2$.

The two mean period differences are the only within-participant information, and they give two equations, $\Delta - \pi - \lambda$ and $-\Delta - \pi$, in three unknowns. No combination of them isolates $\Delta$, or $\pi$, free of $\lambda$. The one estimate of $\Delta$ that is free of $\lambda$ compares period 1 alone (AB period 1 minus BA period 1), a between-participant comparison; the model with a carry-over term, later in this article, does exactly that.

Two effects are aliased when the design gives them the same pattern in the data, so no analysis can separate them.

In the 2x2 crossover, carry-over, a sequence effect and a treatment-by-period interaction (a drug difference that changes between periods) all show up as one quantity: the difference between the sequences in their participant totals (each participant's period-1 reading plus period-2 reading) [1, 2]. The data cannot say which of the three it measures.

The simulated trial below was generated with $\Delta = -8$, $\pi = -2$ and $\lambda = -2$ mmHg. The treatment estimate therefore targets $-8 - (-2)/2 = -7$ mmHg, a bias of 1 mmHg towards no effect. The period estimate targets $-2 + (-2)/2 = -3$ mmHg.

Move $\Delta$, $\pi$ and $\lambda$ and watch the four cell means, the two mean period differences and $\hat\Delta$. The bias of $\hat\Delta$ is always $-\lambda/2$. The defaults are the values used to generate the simulated data.

Why the carry-over test has little power

Carry-over can only be estimated from the participant totals, $t_i = Y_{i1} + Y_{i2}$. Their means differ between the sequences by $\lambda$, and a two-sample t-test on the totals is the carry-over test. Each total, however, contains $2s_i$, so the whole between-participant variation stays in it:

$$\operatorname{SE}(\hat\lambda) = \sqrt{\left(4\sigma_s^2 + 2\sigma_e^2\right)\left(\tfrac{1}{n_{AB}} + \tfrac{1}{n_{BA}}\right)}$$

Here $\sigma_s$ is the between-participant standard deviation (SD), the spread of the levels $s_i$, and $\sigma_e$ is the within-participant SD of a single reading. The counts $n_{AB}$ and $n_{BA}$ are the participants in each sequence. The treatment estimate works on differences, where $\sigma_s$ cancels, so its standard error depends on $\sigma_e$ alone.

With the design values of the simulated trial, $\sigma_s = 12$ mmHg, $\sigma_e = 6$ mmHg and 12 participants per sequence, the standard error of the carry-over estimate is 10.39 mmHg.

Power is the probability that a test declares an effect when an effect of a stated size truly exists. The power of a two-sided test at the 5% level against the true $\lambda = -2$ mmHg is 0.054, barely above the 5% it would have with no carry-over at all. The carry-over that this design would detect in four trials out of five, a usual planning target for power, is about 30.5 mmHg, far larger than the treatment effect.

The carry-over test compares participant totals between sequences, a between-participant comparison with little power. A non-significant result does not show that carry-over is absent; carry-over is prevented by an adequate washout, not removed by a preliminary test.

Testing first, then using period 1 only: why the two-stage procedure fails

A once-standard approach, proposed by Grizzle [4], runs in two stages. First, carry-over is tested on the participant totals at a lenient significance level. If that test is significant, period 2 is discarded and the drugs are compared in period 1 as a parallel-group trial; otherwise, the period differences are used.

Freeman [5] showed that this procedure inflates the type I error, the probability of declaring a drug difference when none exists, and biases the final estimate. The carry-over test and the period-1 comparison share the period-1 readings and are strongly correlated, so the first stage selects which result the second stage reports.

Methodological work now advises against the two-stage procedure [1, 5]. Carry-over is handled at the design stage instead: a condition that is stable over the trial, a washout period between treatments long enough for the first drug's effect to fade, and outcomes measured at the end of each period. The crossover design and analysis guide covers the design side, including the washout, in more depth.

The analysis model: a mixed model with treatment and period

The usual analysis is a linear mixed model: a regression with fixed terms for treatment and period and a random intercept for each participant. It is fitted by REML (restricted maximum likelihood, the usual way to estimate the spread between and within participants in small samples).

The degrees of freedom of a t-test set the exact shape of the t distribution behind its P value and CI; Kenward-Roger degrees of freedom, a small-sample adjustment, are used for this model's t-tests.

With complete data, the treatment coefficient of this model equals $\hat\Delta$ exactly, and its standard error is half the standard error of the pooled-variance (Student) two-sample t-test on the period differences, the version that assumes one variance shared by both sequences. The t statistic, the P value and the 22 degrees of freedom are identical. A Welch t-test, which does not pool the two variances, would give slightly different degrees of freedom.

The model earns its place with less tidy data: it takes baseline covariates and keeps a participant with one missing period. The mixed-model series covers the model in depth.

Adding a carry-over term does not rescue the analysis. With $\lambda$ in the model, the treatment coefficient is estimated from period 1 alone, a between-participant comparison that gives up the crossover's main advantage.

The simulated trial in Stata and R

The dataset (simulated data) has 48 rows, one per participant per period: 24 participants, 12 in each sequence, with systolic blood pressure in mmHg. It was generated with $\Delta = -8$, $\pi = -2$ and $\lambda = -2$ mmHg, a between-participant SD of 12 mmHg and a within-participant SD of 6 mmHg.

The mean period differences are -3.88 mmHg in sequence AB and 12.15 mmHg in sequence BA. Half their difference gives $\hat\Delta = -8.02$ mmHg (standard error 1.79; 95% CI -11.74 to -4.30; P = 0.0002). The mixed model returns the same coefficient, standard error and interval. The period estimate is -4.13 mmHg (95% CI -7.85 to -0.41).

In this one draw, $\hat\Delta$ happens to lie near the true $\Delta$ of -8, although its target is -7. A single trial of 24 participants cannot reveal a bias of 1 mmHg against a standard error of 1.79.

The carry-over estimate is 7.45 mmHg (standard error 5.45; 95% CI -3.85 to 18.75; P = 0.19), against a true value of -2 mmHg. It is far from the truth and not significant, exactly what low power predicts. Its standard error is about half the design value of 10.39 because the fitted between-participant SD in this draw, 5.19 mmHg, is well below the 12 mmHg used to generate the data.

The model with a carry-over term estimates the treatment effect from period 1 only: -4.29 mmHg (95% CI -10.89 to 2.31).

Estimates from the simulated trial (simulated data), in mmHg, beside what each estimator targets under the values used to generate the data.
QuantityEstimate (mmHg)95% CI (mmHg)Target (mmHg)
Mean period difference, sequence AB-3.88not shown-4
Mean period difference, sequence BA12.15not shown10
Treatment effect, period differences-8.02-11.74 to -4.30-7
Treatment effect, mixed model-8.02-11.74 to -4.30-7
Period effect-4.13-7.85 to -0.41-3
Carry-over, totals AB minus BA7.45-3.85 to 18.75-2
Treatment effect, model with carry-over term-4.29-10.89 to 2.31-8

Stata: the whole analysis, from the cell means to the power of the carry-over test

Stata code crossover_stata.do
* Analysis of a 2x2 AB/BA crossover (simulated data, not evidence about any real drug)
* Outcome sbp = systolic blood pressure in mmHg; drug A versus drug B; 12 + 12 participants, two periods.
version 18
clear all
set more off
set linesize 120

* read the simulated dataset (not published with this article)
import delimited using crossover.csv, clear varnames(1)
describe, short
list in 1/6, noobs sepby(id)

* participants per sequence (one row per participant)
egen byte first = tag(id)
tabulate seq if first

* cell means by sequence and period
table seq period, statistic(mean sbp) nformat(%9.2f)

* mixed model: treatment and period as fixed effects, a random intercept per participant
* (REML, Kenward-Roger degrees of freedom). Fitted quietly; the lines below print its treatment
* and period rows and the two standard deviations. R prints the same model in full.
quietly mixed sbp i.trt i.period || id:, reml dfmethod(kroger)
matrix T = r(table)
scalar mm_d = T[rownumb(T, "b"), colnumb(T, "sbp:1.trt")]
scalar mm_se = T[rownumb(T, "se"), colnumb(T, "sbp:1.trt")]
scalar mm_lo = T[rownumb(T, "ll"), colnumb(T, "sbp:1.trt")]
scalar mm_hi = T[rownumb(T, "ul"), colnumb(T, "sbp:1.trt")]
scalar mm_df = T[rownumb(T, "df"), colnumb(T, "sbp:1.trt")]
scalar mm_p = T[rownumb(T, "pvalue"), colnumb(T, "sbp:1.trt")]
scalar mm_pi = T[rownumb(T, "b"), colnumb(T, "sbp:2.period")]
scalar mm_pilo = T[rownumb(T, "ll"), colnumb(T, "sbp:2.period")]
scalar mm_pihi = T[rownumb(T, "ul"), colnumb(T, "sbp:2.period")]
scalar mm_sdid = exp(_b[lns1_1_1:_cons])
scalar mm_sdres = exp(_b[lnsig_e:_cons])
display "mixed model, treatment (A minus B): " %8.4f mm_d "  SE " %7.4f mm_se "  95% CI " %8.4f mm_lo " to " %8.4f mm_hi
display "mixed model, treatment: df " %7.4f mm_df "  P " %6.4f mm_p
display "mixed model, period (2 minus 1): " %8.4f mm_pi "  95% CI " %8.4f mm_pilo " to " %8.4f mm_pihi
display "mixed model, SD between participants " %7.4f mm_sdid "  SD within participants " %7.4f mm_sdres

* mixed model with an explicit carry-over term: treatment is then estimated from period 1 only
quietly mixed sbp i.trt i.period i.carry || id:, reml dfmethod(kroger)
matrix C = r(table)
scalar mmc_d = C[rownumb(C, "b"), colnumb(C, "sbp:1.trt")]
scalar mmc_lo = C[rownumb(C, "ll"), colnumb(C, "sbp:1.trt")]
scalar mmc_hi = C[rownumb(C, "ul"), colnumb(C, "sbp:1.trt")]
scalar mmc_l = C[rownumb(C, "b"), colnumb(C, "sbp:1.carry")]
scalar mmc_lse = C[rownumb(C, "se"), colnumb(C, "sbp:1.carry")]
scalar mmc_lp = C[rownumb(C, "pvalue"), colnumb(C, "sbp:1.carry")]
display "carry-over model, treatment (A minus B): " %8.4f mmc_d "  95% CI " %8.4f mmc_lo " to " %8.4f mmc_hi
display "carry-over model, carry-over term: " %8.4f mmc_l "  SE " %7.4f mmc_lse "  P " %6.4f mmc_lp

* one row per participant: period difference d = period 1 minus period 2, and the participant total
keep id seq period sbp
reshape wide sbp, i(id) j(period)
generate double d = sbp1 - sbp2
generate double tot = sbp1 + sbp2
generate byte grp = cond(seq == "AB", 1, 2)
label define grp 1 "AB" 2 "BA"
label values grp grp

* treatment effect from period differences: Delta-hat = (mean d in AB - mean d in BA)/2, two-sample t-test
ttest d, by(grp)
scalar dbar_ab = r(mu_1)
scalar dbar_ba = r(mu_2)
scalar df_d = r(df_t)
scalar delta = (r(mu_1) - r(mu_2)) / 2
scalar delta_se = r(se) / 2
scalar delta_p = r(p)
scalar tc = invttail(df_d, 0.025)
display "treatment (A minus B): " %8.4f delta "  SE " %7.4f delta_se "  95% CI " %8.4f delta - tc * delta_se " to " %8.4f delta + tc * delta_se "  P " %6.4f delta_p

* period effect (period 2 minus period 1) from the same differences: -(mean d in AB + mean d in BA)/2
scalar per = -(dbar_ab + dbar_ba) / 2
display "period (2 minus 1): " %8.4f per "  95% CI " %8.4f per - tc * delta_se " to " %8.4f per + tc * delta_se

* carry-over test: compare participant totals between sequences (between-participant, so low power)
ttest tot, by(grp)

* design values used to simulate: sigma_s = 12, sigma_e = 6, lambda = -2
* power of the carry-over test at lambda = -2, and the carry-over it detects with 80 percent power
* (two-sided 5 percent level); SE of lambda-hat = sqrt((4 sigma_s^2 + 2 sigma_e^2) (1/n_AB + 1/n_BA))
scalar sig_s = 12
scalar sig_e = 6
scalar lam = -2
count if grp == 1
scalar n_ab = r(N)
count if grp == 2
scalar n_ba = r(N)
scalar des_se = sqrt((4 * sig_s^2 + 2 * sig_e^2) * (1 / n_ab + 1 / n_ba))
scalar des_df = n_ab + n_ba - 2
scalar des_t = invttail(des_df, 0.025)
scalar des_power = nt(des_df, lam / des_se, -des_t) + nttail(des_df, lam / des_se, des_t)
scalar des_mde80 = (des_t + invttail(des_df, 0.20)) * des_se
display "design SE of the carry-over estimate: " %7.4f des_se
display "power of the carry-over test at lambda = -2: " %6.4f des_power
display "carry-over detected with 80% power: " %7.4f des_mde80
Output of the run crossover_stata.log
. * Analysis of a 2x2 AB/BA crossover (simulated data, not evidence about any r
> eal drug)
. * Outcome sbp = systolic blood pressure in mmHg; drug A versus drug B; 12 + 1
> 2 participants, two periods.
. version 18

. clear all

. set more off

. set linesize 120

.
. * read the simulated dataset (not published with this article)
. import delimited using crossover.csv, clear varnames(1)
(encoding automatically selected: ISO-8859-1)
(7 vars, 48 obs)

. describe, short

Contains data
 Observations:            48
    Variables:             7
Sorted by:
     Note: Dataset has changed since last saved.

. list in 1/6, noobs sepby(id)

  +----------------------------------------------------+
  | id   seq   period   treatm~t   trt   carry     sbp |
  |----------------------------------------------------|
  |  1    BA        1          B     0       0   146.2 |
  |  1    BA        2          A     1       0   134.5 |
  |----------------------------------------------------|
  |  2    AB        1          A     1       0   152.4 |
  |  2    AB        2          B     0       1   149.6 |
  |----------------------------------------------------|
  |  3    AB        1          A     1       0   160.7 |
  |  3    AB        2          B     0       1   156.8 |
  +----------------------------------------------------+

.
. * participants per sequence (one row per participant)
. egen byte first = tag(id)

. tabulate seq if first

        seq |      Freq.     Percent        Cum.
------------+-----------------------------------
         AB |         12       50.00       50.00
         BA |         12       50.00      100.00
------------+-----------------------------------
      Total |         24      100.00

.
. * cell means by sequence and period
. table seq period, statistic(mean sbp) nformat(%9.2f)

-----------------------------------
        |           period
        |       1        2    Total
--------+--------------------------
seq     |
  AB    |  145.72   149.60   147.66
  BA    |  150.01   137.86   143.93
  Total |  147.86   143.73   145.80
-----------------------------------

.
. * mixed model: treatment and period as fixed effects, a random intercept per participant
. * (REML, Kenward-Roger degrees of freedom). Fitted quietly; the lines below print its treatment
. * and period rows and the two standard deviations. R prints the same model in full.
. quietly mixed sbp i.trt i.period || id:, reml dfmethod(kroger)

. matrix T = r(table)

. scalar mm_d = T[rownumb(T, "b"), colnumb(T, "sbp:1.trt")]

. scalar mm_se = T[rownumb(T, "se"), colnumb(T, "sbp:1.trt")]

. scalar mm_lo = T[rownumb(T, "ll"), colnumb(T, "sbp:1.trt")]

. scalar mm_hi = T[rownumb(T, "ul"), colnumb(T, "sbp:1.trt")]

. scalar mm_df = T[rownumb(T, "df"), colnumb(T, "sbp:1.trt")]

. scalar mm_p = T[rownumb(T, "pvalue"), colnumb(T, "sbp:1.trt")]

. scalar mm_pi = T[rownumb(T, "b"), colnumb(T, "sbp:2.period")]

. scalar mm_pilo = T[rownumb(T, "ll"), colnumb(T, "sbp:2.period")]

. scalar mm_pihi = T[rownumb(T, "ul"), colnumb(T, "sbp:2.period")]

. scalar mm_sdid = exp(_b[lns1_1_1:_cons])

. scalar mm_sdres = exp(_b[lnsig_e:_cons])

. display "mixed model, treatment (A minus B): " %8.4f mm_d "  SE " %7.4f mm_se "  95% CI " %8.4f mm_lo " to " %8.4f mm_
> hi
mixed model, treatment (A minus B):  -8.0167  SE  1.7938  95% CI -11.7367 to  -4.2966

. display "mixed model, treatment: df " %7.4f mm_df "  P " %6.4f mm_p
mixed model, treatment: df 22.0000  P 0.0002

. display "mixed model, period (2 minus 1): " %8.4f mm_pi "  95% CI " %8.4f mm_pilo " to " %8.4f mm_pihi
mixed model, period (2 minus 1):  -4.1333  95% CI  -7.8534 to  -0.4133

. display "mixed model, SD between participants " %7.4f mm_sdid "  SD within participants " %7.4f mm_sdres
mixed model, SD between participants  5.1869  SD within participants  6.2139

.
. * mixed model with an explicit carry-over term: treatment is then estimated from period 1 only
. quietly mixed sbp i.trt i.period i.carry || id:, reml dfmethod(kroger)

. matrix C = r(table)

. scalar mmc_d = C[rownumb(C, "b"), colnumb(C, "sbp:1.trt")]

. scalar mmc_lo = C[rownumb(C, "ll"), colnumb(C, "sbp:1.trt")]

. scalar mmc_hi = C[rownumb(C, "ul"), colnumb(C, "sbp:1.trt")]

. scalar mmc_l = C[rownumb(C, "b"), colnumb(C, "sbp:1.carry")]

. scalar mmc_lse = C[rownumb(C, "se"), colnumb(C, "sbp:1.carry")]

. scalar mmc_lp = C[rownumb(C, "pvalue"), colnumb(C, "sbp:1.carry")]

. display "carry-over model, treatment (A minus B): " %8.4f mmc_d "  95% CI " %8.4f mmc_lo " to " %8.4f mmc_hi
carry-over model, treatment (A minus B):  -4.2917  95% CI -10.8943 to   2.3110

. display "carry-over model, carry-over term: " %8.4f mmc_l "  SE " %7.4f mmc_lse "  P " %6.4f mmc_lp
carry-over model, carry-over term:   7.4500  SE  5.4483  P 0.1853

.
. * one row per participant: period difference d = period 1 minus period 2, and the participant total
. keep id seq period sbp

. reshape wide sbp, i(id) j(period)
(j = 1 2)

Data                               Long   ->   Wide
-----------------------------------------------------------------------------
Number of observations               48   ->   24
Number of variables                   4   ->   4
j variable (2 values)            period   ->   (dropped)
xij variables:
                                    sbp   ->   sbp1 sbp2
-----------------------------------------------------------------------------

. generate double d = sbp1 - sbp2

. generate double tot = sbp1 + sbp2

. generate byte grp = cond(seq == "AB", 1, 2)

. label define grp 1 "AB" 2 "BA"

. label values grp grp

.
. * treatment effect from period differences: Delta-hat = (mean d in AB - mean d in BA)/2, two-sample t-test
. ttest d, by(grp)

Two-sample t test with equal variances
------------------------------------------------------------------------------
   Group |     Obs        Mean    Std. err.   Std. dev.   [95% conf. interval]
---------+--------------------------------------------------------------------
      AB |      12   -3.883336    2.973998    10.30223   -10.42906     2.66239
      BA |      12       12.15    2.006485    6.950669    7.733753    16.56624
---------+--------------------------------------------------------------------
Combined |      24    4.133331    2.423217    11.87129   -.8794746    9.146137
---------+--------------------------------------------------------------------
    diff |           -16.03333    3.587569                -23.4735   -8.593171
------------------------------------------------------------------------------
    diff = mean(AB) - mean(BA)                                    t =  -4.4691
H0: diff = 0                                     Degrees of freedom =       22

    Ha: diff < 0                 Ha: diff != 0                 Ha: diff > 0
 Pr(T < t) = 0.0001         Pr(|T| > |t|) = 0.0002          Pr(T > t) = 0.9999

. scalar dbar_ab = r(mu_1)

. scalar dbar_ba = r(mu_2)

. scalar df_d = r(df_t)

. scalar delta = (r(mu_1) - r(mu_2)) / 2

. scalar delta_se = r(se) / 2

. scalar delta_p = r(p)

. scalar tc = invttail(df_d, 0.025)

. display "treatment (A minus B): " %8.4f delta "  SE " %7.4f delta_se "  95% CI " %8.4f delta - tc * delta_se " to " %8
> .4f delta + tc * delta_se "  P " %6.4f delta_p
treatment (A minus B):  -8.0167  SE  1.7938  95% CI -11.7367 to  -4.2966  P 0.0002

.
. * period effect (period 2 minus period 1) from the same differences: -(mean d in AB + mean d in BA)/2
. scalar per = -(dbar_ab + dbar_ba) / 2

. display "period (2 minus 1): " %8.4f per "  95% CI " %8.4f per - tc * delta_se " to " %8.4f per + tc * delta_se
period (2 minus 1):  -4.1333  95% CI  -7.8534 to  -0.4132

.
. * carry-over test: compare participant totals between sequences (between-participant, so low power)
. ttest tot, by(grp)

Two-sample t test with equal variances
------------------------------------------------------------------------------
   Group |     Obs        Mean    Std. err.   Std. dev.   [95% conf. interval]
---------+--------------------------------------------------------------------
      AB |      12    295.3167    3.709261    12.84926    287.1526    303.4807
      BA |      12    287.8667    3.990715    13.82424    279.0832    296.6502
---------+--------------------------------------------------------------------
Combined |      24    291.5917      2.7752    13.59565    285.8507    297.3326
---------+--------------------------------------------------------------------
    diff |            7.449999    5.448341               -3.849168    18.74917
------------------------------------------------------------------------------
    diff = mean(AB) - mean(BA)                                    t =   1.3674
H0: diff = 0                                     Degrees of freedom =       22

    Ha: diff < 0                 Ha: diff != 0                 Ha: diff > 0
 Pr(T < t) = 0.9073         Pr(|T| > |t|) = 0.1853          Pr(T > t) = 0.0927

.
. * design values used to simulate: sigma_s = 12, sigma_e = 6, lambda = -2
. * power of the carry-over test at lambda = -2, and the carry-over it detects with 80 percent power
. * (two-sided 5 percent level); SE of lambda-hat = sqrt((4 sigma_s^2 + 2 sigma_e^2) (1/n_AB + 1/n_BA))
. scalar sig_s = 12

. scalar sig_e = 6

. scalar lam = -2

. count if grp == 1
  12

. scalar n_ab = r(N)

. count if grp == 2
  12

. scalar n_ba = r(N)

. scalar des_se = sqrt((4 * sig_s^2 + 2 * sig_e^2) * (1 / n_ab + 1 / n_ba))

. scalar des_df = n_ab + n_ba - 2

. scalar des_t = invttail(des_df, 0.025)

. scalar des_power = nt(des_df, lam / des_se, -des_t) + nttail(des_df, lam / des_se, des_t)

. scalar des_mde80 = (des_t + invttail(des_df, 0.20)) * des_se

. display "design SE of the carry-over estimate: " %7.4f des_se
design SE of the carry-over estimate: 10.3923

. display "power of the carry-over test at lambda = -2: " %6.4f des_power
power of the carry-over test at lambda = -2: 0.0539

. display "carry-over detected with 80% power: " %7.4f des_mde80
carry-over detected with 80% power: 30.4717
Stata code and its full log (simulated data). Stata fits the two mixed models without printing its full table, then prints their treatment, period and carry-over rows; R below prints both models in full. The t-test on d prints the difference in mean period differences, which is twice $\hat\Delta$, with twice its standard error. The t-test on totals is the carry-over test. The last lines compute the power of that test, and the carry-over it would detect in four trials out of five, from the design values used to simulate the data: a between-participant SD of 12 mmHg, a within-participant SD of 6 mmHg and a carry-over of -2 mmHg. The 48-row dataset is simulated and not published; the code shows every step, so run the same commands on your own crossover data and compare the structure of the output, not these numbers.

R: the same analysis with lme4 and lmerTest

R code crossover_r.R
# Analysis of a 2x2 AB/BA crossover (simulated data, not evidence about any real drug)
# Outcome sbp = systolic blood pressure in mmHg; drug A versus drug B; 12 + 12 participants, two periods.
# Packages: lme4 and lmerTest, plus pbkrtest for the Kenward-Roger degrees of freedom.
suppressPackageStartupMessages(library(lmerTest))

# read the simulated dataset (not published with this article)
dat <- read.csv("crossover.csv")
str(dat)
print(head(dat, 6))

# participants per sequence (one row per participant)
ids <- unique(dat[, c("id", "seq")])
print(table(ids$seq))

# cell means by sequence and period
print(round(tapply(dat$sbp, list(seq = dat$seq, period = dat$period), mean), 2))

# mixed model: treatment and period as fixed effects, a random intercept per participant (REML, Kenward-Roger df)
dat$trt <- factor(dat$trt)
dat$period <- factor(dat$period)
dat$carry <- factor(dat$carry)
mm <- lmer(sbp ~ trt + period + (1 | id), data = dat, REML = TRUE)
smm <- summary(mm, ddf = "Kenward-Roger")
print(smm)
cm <- coef(smm)
tq <- qt(0.975, cm["trt1", "df"])
# 95% CI of the treatment coefficient (A minus B)
print(round(cm["trt1", "Estimate"] + c(lo = -1, hi = 1) * tq * cm["trt1", "Std. Error"], 4))

# mixed model with an explicit carry-over term: treatment is then estimated from period 1 only
mmc <- lmer(sbp ~ trt + period + carry + (1 | id), data = dat, REML = TRUE)
smmc <- summary(mmc, ddf = "Kenward-Roger")
print(smmc)
cc <- coef(smmc)
tqc <- qt(0.975, cc["trt1", "df"])
# 95% CI of the treatment coefficient in this model
print(round(cc["trt1", "Estimate"] + c(lo = -1, hi = 1) * tqc * cc["trt1", "Std. Error"], 4))

# one row per participant: period difference d = period 1 minus period 2, and the participant total
w <- reshape(dat[, c("id", "seq", "period", "sbp")], idvar = c("id", "seq"),
             timevar = "period", direction = "wide")
w$d <- w$sbp.1 - w$sbp.2
w$tot <- w$sbp.1 + w$sbp.2
w$seq <- factor(w$seq, levels = c("AB", "BA"))

# treatment effect from period differences: Delta-hat = (mean d in AB - mean d in BA)/2, two-sample t-test
td <- t.test(d ~ seq, data = w, var.equal = TRUE)
print(td)
dbar_ab <- unname(td$estimate[1])
dbar_ba <- unname(td$estimate[2])
delta <- (dbar_ab - dbar_ba) / 2
delta_se <- td$stderr / 2
tc <- qt(0.975, unname(td$parameter))
cat(sprintf("treatment (A minus B): %.4f, SE %.4f, 95%% CI %.4f to %.4f\n",
            delta, delta_se, delta - tc * delta_se, delta + tc * delta_se))

# period effect (period 2 minus period 1) from the same differences: -(mean d in AB + mean d in BA)/2
per <- -(dbar_ab + dbar_ba) / 2
cat(sprintf("period (2 minus 1): %.4f, 95%% CI %.4f to %.4f\n", per, per - tc * delta_se, per + tc * delta_se))

# carry-over test: compare participant totals between sequences (between-participant, so low power)
tt <- t.test(tot ~ seq, data = w, var.equal = TRUE)
print(tt)

# design values used to simulate: sigma_s = 12, sigma_e = 6, lambda = -2
# power of the carry-over test at lambda = -2, and the carry-over it detects with 80 percent power
# (two-sided 5 percent level); SE of lambda-hat = sqrt((4 sigma_s^2 + 2 sigma_e^2) (1/n_AB + 1/n_BA))
sigma_s <- 12; sigma_e <- 6; lambda <- -2
n_ab <- sum(ids$seq == "AB"); n_ba <- sum(ids$seq == "BA")
des_se <- sqrt((4 * sigma_s^2 + 2 * sigma_e^2) * (1 / n_ab + 1 / n_ba))
des_df <- n_ab + n_ba - 2
des_t <- qt(0.975, des_df)
des_power <- pt(-des_t, des_df, lambda / des_se) + 1 - pt(des_t, des_df, lambda / des_se)
des_mde80 <- (des_t + qt(0.80, des_df)) * des_se
cat(sprintf("design SE of the carry-over estimate: %.4f\n", des_se))
cat(sprintf("power of the carry-over test at lambda = -2: %.4f\n", des_power))
cat(sprintf("carry-over detected with 80%% power: %.4f\n", des_mde80))
Output of the run crossover_r.log
> suppressPackageStartupMessages(library(lmerTest))

> dat <- read.csv("crossover.csv")

> str(dat)
'data.frame':	48 obs. of  7 variables:
 $ id       : int  1 1 2 2 3 3 4 4 5 5 ...
 $ seq      : chr  "BA" "BA" "AB" "AB" ...
 $ period   : int  1 2 1 2 1 2 1 2 1 2 ...
 $ treatment: chr  "B" "A" "A" "B" ...
 $ trt      : int  0 1 1 0 1 0 1 0 1 0 ...
 $ carry    : int  0 0 0 1 0 1 0 1 0 1 ...
 $ sbp      : num  146 134 152 150 161 ...

> print(head(dat, 6))
  id seq period treatment trt carry   sbp
1  1  BA      1         B   0     0 146.2
2  1  BA      2         A   1     0 134.5
3  2  AB      1         A   1     0 152.4
4  2  AB      2         B   0     1 149.6
5  3  AB      1         A   1     0 160.7
6  3  AB      2         B   0     1 156.8

> ids <- unique(dat[, c("id", "seq")])

> print(table(ids$seq))

AB BA
12 12

> print(round(tapply(dat$sbp, list(seq = dat$seq, period = dat$period),
+     mean), 2))
    period
seq       1      2
  AB 145.72 149.60
  BA 150.01 137.86

> dat$trt <- factor(dat$trt)

> dat$period <- factor(dat$period)

> dat$carry <- factor(dat$carry)

> mm <- lmer(sbp ~ trt + period + (1 | id), data = dat,
+     REML = TRUE)

> smm <- summary(mm, ddf = "Kenward-Roger")

> print(smm)
Linear mixed model fit by REML. t-tests use Kenward-Roger's method [
lmerModLmerTest]
Formula: sbp ~ trt + period + (1 | id)
   Data: dat

REML criterion at convergence: 321

Scaled residuals:
    Min      1Q  Median      3Q     Max
-2.0487 -0.4621 -0.1869  0.6003  1.4972

Random effects:
 Groups   Name        Variance Std.Dev.
 id       (Intercept) 26.90    5.187
 Residual             38.61    6.214
Number of obs: 48, groups:  id, 24

Fixed effects:
            Estimate Std. Error      df t value Pr(>|t|)
(Intercept)  151.871      1.880  44.797  80.784  < 2e-16 ***
trt1          -8.017      1.794  22.000  -4.469 0.000192 ***
period2       -4.133      1.794  22.000  -2.304 0.031028 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
        (Intr) trt1
trt1    -0.477
period2 -0.477  0.000

> cm <- coef(smm)

> tq <- qt(0.975, cm["trt1", "df"])

> print(round(cm["trt1", "Estimate"] + c(lo = -1, hi = 1) *
+     tq * cm["trt1", "Std. Error"], 4))
      lo       hi
-11.7367  -4.2966

> mmc <- lmer(sbp ~ trt + period + carry + (1 | id),
+     data = dat, REML = TRUE)

> smmc <- summary(mmc, ddf = "Kenward-Roger")

> print(smmc)
Linear mixed model fit by REML. t-tests use Kenward-Roger's method [
lmerModLmerTest]
Formula: sbp ~ trt + period + carry + (1 | id)
   Data: dat

REML criterion at convergence: 313.9

Scaled residuals:
    Min      1Q  Median      3Q     Max
-2.1949 -0.5026 -0.0767  0.4785  1.5487

Random effects:
 Groups   Name        Variance Std.Dev.
 id       (Intercept) 25.22    5.022
 Residual             38.61    6.214
Number of obs: 48, groups:  id, 24

Fixed effects:
            Estimate Std. Error      df t value Pr(>|t|)
(Intercept)  150.008      2.306  38.059  65.041   <2e-16 ***
trt1          -4.292      3.262  38.059  -1.316   0.1961
period2       -7.858      3.262  38.059  -2.409   0.0209 *
carry1         7.450      5.448  22.000   1.367   0.1853
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
        (Intr) trt1   perid2
trt1    -0.707
period2  0.279 -0.698
carry1  -0.591  0.835 -0.835

> cc <- coef(smmc)

> tqc <- qt(0.975, cc["trt1", "df"])

> print(round(cc["trt1", "Estimate"] + c(lo = -1, hi = 1) *
+     tqc * cc["trt1", "Std. Error"], 4))
      lo       hi
-10.8943   2.3110

> w <- reshape(dat[, c("id", "seq", "period", "sbp")],
+     idvar = c("id", "seq"), timevar = "period", direction = "wide")

> w$d <- w$sbp.1 - w$sbp.2

> w$tot <- w$sbp.1 + w$sbp.2

> w$seq <- factor(w$seq, levels = c("AB", "BA"))

> td <- t.test(d ~ seq, data = w, var.equal = TRUE)

> print(td)

	Two Sample t-test

data:  d by seq
t = -4.4691, df = 22, p-value = 0.0001918
alternative hypothesis: true difference in means between group AB and group BA is not equal to 0
95 percent confidence interval:
 -23.473498  -8.593169
sample estimates:
mean in group AB mean in group BA
       -3.883333        12.150000


> dbar_ab <- unname(td$estimate[1])

> dbar_ba <- unname(td$estimate[2])

> delta <- (dbar_ab - dbar_ba)/2

> delta_se <- td$stderr/2

> tc <- qt(0.975, unname(td$parameter))

> cat(sprintf("treatment (A minus B): %.4f, SE %.4f, 95%% CI %.4f to %.4f\n",
+     delta, delta_se, delta - tc * delta_se, delta + tc * delta_se))
treatment (A minus B): -8.0167, SE 1.7938, 95% CI -11.7367 to -4.2966

> per <- -(dbar_ab + dbar_ba)/2

> cat(sprintf("period (2 minus 1): %.4f, 95%% CI %.4f to %.4f\n",
+     per, per - tc * delta_se, per + tc * delta_se))
period (2 minus 1): -4.1333, 95% CI -7.8534 to -0.4133

> tt <- t.test(tot ~ seq, data = w, var.equal = TRUE)

> print(tt)

	Two Sample t-test

data:  tot by seq
t = 1.3674, df = 22, p-value = 0.1853
alternative hypothesis: true difference in means between group AB and group BA is not equal to 0
95 percent confidence interval:
 -3.849168 18.749168
sample estimates:
mean in group AB mean in group BA
        295.3167         287.8667


> sigma_s <- 12

> sigma_e <- 6

> lambda <- -2

> n_ab <- sum(ids$seq == "AB")

> n_ba <- sum(ids$seq == "BA")

> des_se <- sqrt((4 * sigma_s^2 + 2 * sigma_e^2) * (1/n_ab +
+     1/n_ba))

> des_df <- n_ab + n_ba - 2

> des_t <- qt(0.975, des_df)

> des_power <- pt(-des_t, des_df, lambda/des_se) + 1 -
+     pt(des_t, des_df, lambda/des_se)

> des_mde80 <- (des_t + qt(0.8, des_df)) * des_se

> cat(sprintf("design SE of the carry-over estimate: %.4f\n",
+     des_se))
design SE of the carry-over estimate: 10.3923

> cat(sprintf("power of the carry-over test at lambda = -2: %.4f\n",
+     des_power))
power of the carry-over test at lambda = -2: 0.0539

> cat(sprintf("carry-over detected with 80%% power: %.4f\n",
+     des_mde80))
carry-over detected with 80% power: 30.4717
R code and its full log (simulated data). With var.equal = TRUE, t.test gives the pooled-variance test whose t statistic, P value and degrees of freedom match the mixed model. Its confidence interval is for the difference in mean period differences, which is twice $\hat\Delta$. The last lines compute the power of the carry-over test, and the carry-over it would detect in four trials out of five, from the design values used to simulate the data. The 48-row dataset is simulated and not published; the code shows every step, so run the same commands on your own crossover data and compare the structure of the output, not these numbers.

Pitfalls

  • Crediting a period-2 change to the period-2 drug

    The change from period 1 to period 2 inside one sequence is read as the effect of the drug taken in period 2, although part of it may belong to period 2 itself.

    Fix: Compare the sequences, not the periods within one sequence. Half the difference between the two mean period differences removes any shift shared by period 2.

  • Treating the sequence effect as a separate effect that can be estimated

    In a 2x2 crossover, a difference between the AB and BA groups cannot be split into sequence, carry-over and treatment-by-period interaction.

    Fix: Report the comparison of totals as one aliased quantity, and prevent carry-over by design.

  • "A non-significant carry-over test shows there is no carry-over."

    In the simulated trial the test gives P = 0.19, although a carry-over of -2 mmHg was built into the data.

    Fix: The carry-over test compares participant totals between sequences, a between-participant comparison with little power. A non-significant result does not show that carry-over is absent; carry-over is prevented by an adequate washout, not removed by a preliminary test.

  • Testing carry-over first, then falling back to period 1

    The two-stage procedure inflates the type I error and biases the estimate of the treatment effect [5].

    Fix: Prespecify the period-difference analysis, or the equivalent mixed model, and plan an adequate washout.

What to do in your own analysis

Glossary

crossover trial (การศึกษาแบบไขว้)
A trial in which each participant receives every treatment under study, one per period, in a randomised order.
sequence (ลำดับ)
The order in which a participant receives the treatments; in a 2x2 crossover, AB or BA.
period effect (ผลของช่วงเวลา)
A shift in the outcome shared by all participants in one period compared with another, whichever treatment they take.
carry-over effect (ผลตกค้าง)
An effect of an earlier period's treatment that persists into a later period.
washout period (ช่วงล้างยา)
A treatment-free interval between periods that lets the first treatment's effect fade.
aliasing (การปนกันของผล)
Two effects are aliased when the design gives them the same pattern in the data, so no analysis can separate them.
within-participant comparison (การเปรียบเทียบภายในผู้เข้าร่วม)
A comparison of readings taken on the same person, so stable differences between people cancel.
power (อำนาจการทดสอบ)
The probability that a test declares an effect when an effect of a stated size truly exists.
CONSORT extension (ส่วนขยายของ CONSORT)
A version of the CONSORT reporting guideline for randomised trials adapted to one design; the crossover version appeared in 2019.

References

  1. Senn S. Cross-over trials in clinical research. 2nd ed. Wiley; 2002. https://doi.org/10.1002/0470854596
  2. Jones B, Kenward MG. Design and analysis of cross-over trials. 3rd ed. Chapman and Hall/CRC; 2014. https://doi.org/10.1201/b17537
  3. Hills M, Armitage P. The two-period cross-over clinical trial. Br J Clin Pharmacol. 1979;8(1):7-20. https://doi.org/10.1111/j.1365-2125.1979.tb05903.x
  4. Grizzle JE. The two-period change-over design and its use in clinical trials. Biometrics. 1965;21(2):467-80. https://doi.org/10.2307/2528104
  5. Freeman PR. The performance of the two-stage analysis of two-treatment, two-period crossover trials. Stat Med. 1989;8(12):1421-32. https://doi.org/10.1002/sim.4780081202
  6. Dwan K, Li T, Altman DG, Elbourne D. CONSORT 2010 statement: extension to randomised crossover trials. BMJ. 2019;366:l4378. https://doi.org/10.1136/bmj.l4378

Key takeaways

  • In a 2x2 crossover, half the difference between the two sequences' mean period differences estimates the treatment effect, and any period effect cancels.
  • With complete data, a random-intercept mixed model with treatment and period estimates half the difference in mean period differences, with the P value and degrees of freedom of the pooled-variance t-test.
  • Carry-over is aliased with sequence: it shifts the treatment estimate by minus half its size and can be estimated only from comparisons between participants.
  • The carry-over test has little power, so a non-significant result does not show that carry-over is absent.
  • Prevent carry-over with an adequate washout, and do not let a preliminary carry-over test choose the analysis.

Related in the wiki: [[treatment-by-time-interaction-longitudinal-trials]] [[crossover-trials-design-analysis-guide]] [[mixed-model-series-guide]]

0
Message for International and Thai ReadersUnderstanding My Medical Context in ThailandRead more →Message for International and Thai ReadersUnderstanding My Broader Content Beyond MedicineRead more →

Comments

No comments yet. Be the first to share your thoughts.

Sign in to comment