Compare treatments within the same cow

Class 10: a complete crossover analysis

Gavin Simpson

Aarhus University

Wednesday, 11 November 2026

Start with the cow study, not a formula

13:55–14:05. In the synthetic study, 24 cows were assigned to AB or BA. Every cow was measured under both rations.

The question is: how much does mean intake differ under B compared with A, in kg/day?

We assume the earlier ration no longer affects intake when the next outcome is recorded. The outcome file cannot establish this assumption.

Read the cow-period records

library("lme4")
library("performance")
library("marginaleffects")
x <- readr::read_csv("data/cow-crossover.csv")
head(x, 4)

Each row is a cow’s average intake over the final three days of one treatment period.

Read the cow-period records

# A tibble: 4 × 5
  cow   sequence period treatment intake_kg_day
  <chr> <chr>     <dbl> <chr>             <dbl>
1 W01   AB            1 A                  21.6
2 W01   AB            2 B                  22.7
3 W02   BA            1 B                  22.5
4 W02   BA            2 A                  20.0

Preserve the assigned order and period

library("dplyr")

x <- x |>
  mutate(
    cow = factor(cow),
    sequence = factor(sequence, levels = c("AB", "BA")),
    period = factor(period, levels = c(1, 2)),
    treatment = factor(treatment, levels = c("A", "B"))
  )
x |>
  distinct(cow, sequence) |>
  count(sequence)
# A tibble: 2 × 2
  sequence     n
  <fct>    <int>
1 AB          12
2 BA          12

Twelve cows were randomly assigned to each order. Each order must remain in the data even when we rearrange the rows.

Check that both treatments occur in both periods

library("dplyr")

x |> count(period, treatment)
# A tibble: 4 × 3
  period treatment     n
  <fct>  <fct>     <int>
1 1      A            12
2 1      B            12
3 2      A            12
4 2      B            12
c(cows = nlevels(x$cow), cow_period_rows = nrow(x))
           cows cow_period_rows 
             24              48

Researchers measured the same 24 cows under A and under B. Each cow’s two outcomes form a pair; the two treatment columns do not describe different sets of cows.

Put each cow’s two observations side by side

tidyr::pivot_wider() rearranges existing observations. It does not fit a model.

library("dplyr")

paired <- x |>
  select(cow, treatment, intake_kg_day) |>
  tidyr::pivot_wider(
    names_from = treatment,
    values_from = intake_kg_day,
    names_prefix = "intake_kg_day."
  )
head(paired, 3)
# A tibble: 3 × 3
  cow   intake_kg_day.A intake_kg_day.B
  <fct>           <dbl>           <dbl>
1 W01              21.6            22.7
2 W02              20.0            22.5
3 W03              19.3            21.5

We can now compare A with B for the same cow. The full model will also allow for period and assigned order.

Look at the paired observations

plot(
  paired$intake_kg_day.A,
  paired$intake_kg_day.B,
  xlab = "Same cow under A (kg/day)",
  ylab = "Same cow under B (kg/day)",
  pch = 16
)
abline(a = 0, b = 1, lty = 2)

Above the dashed line means the cow ate more under B. Why might pairing help when usual intake differs between cows?

Build the model from the design

crossover_fit <- lmer(
  intake_kg_day ~ treatment + period + sequence + (1 | cow),
  data = x,
  REML = TRUE
)
summary(crossover_fit)

(1 | cow) allows cows’ usual differences in appetite.

period allows a common first-to-second-period change. sequence allows an average difference between the two assigned-order groups. That group difference may reflect chance differences among cows; it does not identify a lingering treatment effect.

Build the model from the design

Linear mixed model fit by REML ['lmerMod']
Formula: intake_kg_day ~ treatment + period + sequence + (1 | cow)
   Data: x

REML criterion at convergence: 180.8

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-1.70501 -0.37577  0.03845  0.41809  1.80745 

Random effects:
 Groups   Name        Variance Std.Dev.
 cow      (Intercept) 5.0542   2.2482  
 Residual             0.7033   0.8386  
Number of obs: 48, groups:  cow, 24

Fixed effects:
            Estimate Std. Error t value
(Intercept)  22.4265     0.6927  32.377
treatmentB    1.0401     0.2421   4.296
period2      -0.7238     0.2421  -2.990
sequenceBA   -2.2258     0.9492  -2.345

Correlation of Fixed Effects:
           (Intr) trtmnB perid2
treatmentB -0.175              
period2    -0.175  0.000       
sequenceBA -0.685  0.000  0.000

Ask the model for B minus A

crossover_gap <- avg_comparisons(
  crossover_fit,
  variables = "treatment",
  re.form = NA
)
crossover_gap

The estimate is in kg/day. The interval describes uncertainty in the mean comparison, after allowing for period and assigned order.

re.form = NA sets each cow’s fitted intake difference to zero. Both periods and assigned orders receive equal weight in these complete study records.

Ask the model for B minus A


 Estimate Std. Error   z Pr(>|z|)    S 2.5 % 97.5 %
     1.04      0.242 4.3   <0.001 15.8 0.566   1.51

Term: treatment
Type: response
Comparison: B - A

Predict mean intake under each ration

crossover_means <- avg_predictions(
  crossover_fit,
  by = "treatment",
  re.form = NA
)
crossover_means

These are mean intakes, averaging over the periods and assigned orders in this balanced dataset. They are not ranges for an individual cow’s future intake.

Predict mean intake under each ration


 treatment Estimate Std. Error    z Pr(>|z|)   S 2.5 % 97.5 %
         A       21       0.49 42.8   <0.001 Inf    20   21.9
         B       22       0.49 44.9   <0.001 Inf    21   23.0

Type: response

First check that the calculation worked

check_convergence(crossover_fit)
[1] TRUE
attr(,"gradient")
[1] 6.234418e-08
check_singularity(crossover_fit)
[1] FALSE

TRUE on the first line: fitting completed successfully.

FALSE on the second: the estimated variation in cows’ usual intakes did not collapse to zero.

Neither result establishes that the washout was adequate.

Check unexplained differences visually

print(
  check_model(
    crossover_fit,
    check = c("linearity", "qq"),
    residual_type = "normal"
  )
)

Look for a pattern or changing spread against fitted values. In the other panel, compare the points with the reference line.

A crucial limit: the model cannot create washout

Including sequence records and allows for the assigned order. It does not remove an earlier ration’s lingering influence.

With just two periods, we cannot reliably separate a treatment’s current effect from all possible lingering effects using these outcomes alone.

Failing to detect an order difference does not establish that carryover is absent.

Write a recommendation for this study

The researcher wants to know whether ration B increases intake under conditions where earlier treatment effects have worn off.

In your group, write three sentences:

  1. Give the B–A estimate and interval in kg/day.
  2. Explain why period and cow appear in the model.
  3. State the scientific assumption that the output cannot check.

14:40–14:55; then break 14:55–15:10.

Project checkpoint — 25 minutes

15:10–15:35. Use your own group’s project case.

Save the analysis plan before fitting:

  1. Name the outcome and one precise comparison.
  2. Map treatment assignments, shared conditions and repeated measurements.
  3. Explain the model and two checks in words.
  4. Give one alternative and your next action.

Compare designs before choosing the analysis

A calf stays on one diet and is weighed five times. A cow receives both rations in two periods.

For each study, report: what was randomised, what a row means, and what links observations from the same animal.

What additional scientific assumption does the cow study require?

Individual exit — finish at 16:00

A colleague says: “There are 48 cow-period records, so the crossover has 48 separately chosen treatment orders.” Correct the claim.

Another says: “Adding sequence means washout no longer matters.” Explain the error.

After class: 60 minutes, including the calf change comparison and two growth-curve sketches for Class 11.