## Two‐at‐a‐Time analysis in R [Design Issues]

Thank you Helmut for your quick and detailed response!

I'm glad to hear that pseudo-periods have gone the way of the dodo. I agree with the EMA Biostatistics Working Party that it makes more sense (statistically) to maintain the original codings when leaving out data from irrelevant treatments (not to mention being less prone to human error than recoding). Also, I love their name; I bet the EMA Biostatistics Working Party could form a coalition government with the Slightly Silly Party, the Surprise Party, and the Rent Is Too Damn High Party. But I digress.

I don't have WinNonlin, but I tried analyzing Patterson & Jones' Example 4.5 with the Two‐at‐a‐Time Principle using R 1.1.456. My results were similar to the results you posted above, but not as close as I would have expected. I will post the code & table below. Do you see anything I am doing wrong or any obvious explanation of the differences? I don't know how WinNonlin handles missing values, so I tried a few variations (subjects with complete cases in the full dataset, and subjects with complete cases in each of the 3 incomplete-block subsets), but the differences persisted.

 # 3-treatment design from Patterson & Jones (2017) dta = read.table("exam45.dat", header=T) dta[dta==99999] = NA dta$Subj = factor(dta$subject) dta$Per = factor(dta$period) dta$Seq = factor(dta$sequence) dta$Trt = factor(dta$formula) options(contrasts=c("contr.treatment","contr.poly"), digits=4) # Data frame to store results IBD = data.frame(Test = paste(rep(c("T vs. R", "S vs. R", "T vs. S"), each=2),                               rep(c(" AUC", "Cmax"), times=3)),                  PE = NA, LCI = NA, UCI = NA, CV = NA) # AUC T vs R muddle = lm(log(AUC)~Trt+Per+Seq+Subj, data=dta[dta$Trt!="S",]) IBD[1, 2:5] = c(100*exp(coef(muddle)["TrtT"]), 100*exp(confint(muddle,c("TrtT"), level=.9)), 100*sqrt(exp(summary(muddle)$sigma^2)-1)) # Cmax T vs R muddle = lm(log(CMAX)~Trt+Per+Seq+Subj, data=dta[dta$Trt!="S",]) IBD[2, 2:5] = c(100*exp(coef(muddle)["TrtT"]), 100*exp(confint(muddle,c("TrtT"), level=.9)), 100*sqrt(exp(summary(muddle)$sigma^2)-1)) # AUC S vs R muddle = lm(log(AUC)~Trt+Per+Seq+Subj, data=dta[dta$Trt!="T",]) IBD[3, 2:5] = c(100*exp(coef(muddle)["TrtS"]), 100*exp(confint(muddle,c("TrtS"), level=.9)), 100*sqrt(exp(summary(muddle)$sigma^2)-1)) # Cmax S vs R muddle = lm(log(CMAX)~Trt+Per+Seq+Subj, data=dta[dta$Trt!="T",]) IBD[4, 2:5] = c(100*exp(coef(muddle)["TrtS"]), 100*exp(confint(muddle,c("TrtS"), level=.9)), 100*sqrt(exp(summary(muddle)$sigma^2)-1)) # AUC T vs S muddle = lm(log(AUC)~Trt+Per+Seq+Subj, data=dta[dta$Trt!="R",]) IBD[5, 2:5] = c(100*exp(coef(muddle)["TrtT"]), 100*exp(confint(muddle,c("TrtT"), level=.9)), 100*sqrt(exp(summary(muddle)$sigma^2)-1)) # Cmax T vs S muddle = lm(log(CMAX)~Trt+Per+Seq+Subj, data=dta[dta$Trt!="R",]) IBD[6, 2:5] = c(100*exp(coef(muddle)["TrtT"]), 100*exp(confint(muddle,c("TrtT"), level=.9)), 100*sqrt(exp(summary(muddle)$sigma^2)-1)) # Print table print(IBD) 

This produced the following table:

           Test     PE    LCI    UCI    CV 1 T vs. R  AUC 116.14 108.98 123.78 20.88 2 T vs. R Cmax 129.82 119.50 141.04 27.87 3 S vs. R  AUC 140.73 131.02 151.15 23.50 4 S vs. R Cmax 160.04 144.76 176.94 34.06 5 T vs. S  AUC  83.51  78.73  88.58 19.07 6 T vs. S Cmax  82.30  75.96  89.18 26.67