### Pairwise differences in cross-over design

November 8, 2016, 11:13 am A colleague asked me how to present the results of a study in gastric emptying. As you have seen multiple time here, the emptying time t50 is the target variable. But read on even if you work with blood pressure (also often misused as example), rehabilitation, or pharmacology. The methods are the same.

12 subjects received 4 different meals on separate occasions in random order. The gastric half-emptying time t50 measured with MRI techniques was the target variable. How do I compute the confidence intervals of t50 and the 6 possible pairwise differences between meals with their confidence intervals?

### Generate simulated study results

See the source code how to generate the data.

'data.frame':   48 obs. of  3 variables:
$subject: Factor w/ 12 levels "S01","S02","S03",..: 1 1 1 1 2 2 2 2 3 3 ...$ meal   : Factor w/ 4 levels "a","b","c","d": 1 2 3 4 1 2 3 4 1 2 ...
$t50 : num 80.8 84.5 85.3 83.1 82.3 ... ## Mixed model and contrasts with confidence intervals # Contrasts to get confidence intervals of emptying times absContrast = rbind("a" = c(1,0,0,0), "b" = c(1,1,0,0), "c" = c(1,0,1,0), "d" = c(1,0,1,0)) # Linear mixed model fit with meal as fixed and # subject as random # lme is a function in package nlme d_lme = lme(t50~meal, random = ~1|subject, data = d) summary(d_lme)  Linear mixed-effects model fit by REML Data: d AIC BIC logLik 212 223 -100 Random effects: Formula: ~1 | subject (Intercept) Residual StdDev: 0.78 2 Fixed effects: t50 ~ meal Value Std.Error DF t-value p-value (Intercept) 81 0.61 33 132 0.0000 mealb 2 0.81 33 2 0.0417 mealc 2 0.81 33 2 0.0334 meald 3 0.81 33 4 0.0009 Correlation: (Intr) mealb mealc mealb -0.66 mealc -0.66 0.50 meald -0.66 0.50 0.50 Standardized Within-Group Residuals: Min Q1 Med Q3 Max -2.01 -0.57 0.14 0.72 1.71 Number of Observations: 48 Number of Groups: 12  ### t50 with confidence intervals When the design is complete, i.e. no data missing, the estimate of t50 is equal to the treatment mean. The 95% confidence intervals give an approximate range where future means of studies will be expected. Note that I avoid the formulation "where with 95% percent probability the result will be expected"; watch my favorite youtube video to understand the subtle differences. # glht is a function in package multcomp glh_abs = glht(d_lme, linfct = absContrast) ci_abs = data.frame(round(confint(glh_abs)$confint,1))
ci_abs
  Estimate lwr upr
a       81  79  82
b       82  81  84
c       83  81  84
d       83  81  84

### Pairwise differences with confidence intervals

The pairwise differences are computed with the Tukey correction.

# mcp is a function in package multcomp
tukey_contrast = mcp(meal = "Tukey")
glh_tukey = glht(d_lme, linfct = tukey_contrast)
ci_tukey = data.frame(round(confint(glh_tukey)$confint,1)) ci_tukey  Estimate lwr upr b - a 1.7 -0.4 3.8 c - a 1.8 -0.3 3.9 d - a 3.0 0.9 5.0 c - b 0.1 -2.0 2.2 d - b 1.2 -0.8 3.3 d - c 1.2 -0.9 3.2 See the vignette coming with package multcomp for alternatives, such as the Dunnett correction; it is useful if you plan to compare the treatments against one reference treatment, here a, and gives slightly narrower confidence interval. dunnet_contrast = mcp(meal = "Dunnet") glh_dunnet = glht(d_lme, linfct = dunnet_contrast) ci_dunnet = data.frame(round(confint(glh_dunnet)$confint,1))
ci_dunnet
      Estimate  lwr upr
b - a      1.7 -0.2 3.6
c - a      1.8 -0.1 3.7
d - a      3.0  1.1 4.9

See the paper by Gelman, Hill and Yajima why you could also not worry about multiple corrections when using a function like lme for estimation of the contrasts. Some reviewers will object, though, but give it a try.

### How to put that into a publication

I leave it to you to neatly arrange the pairwise differences in a triangular table or use the results inline.

### And which pairing is significant?

Only one pair is "significant": both lower and upper limit have the same sign.

### But my boss wants real p-values!

No problems, these can easily be calculated (hint: summary(glh_tukey)$test$pvalues), but I do not show the details here so it gets harder for you to produce the horrible p-values table. Once there are p-values, nobody will look at the confidence intervals of the differences. Better journals with statistical reviewers will applaud when you only present confidence intervals, lesser journals with medical reviewers who have build their reputation on p-values will ask for you to p. In case you do not get the "reputation" part: you have not watched the youtube video video yet. 