Bioequivalence and Bioavailability Forum

Dear ElMaestro, dear All,

❝ Bonus questions:

❝ 1. Is the estimated Swr unbiased?

❝ 2. If yes, is it even reflecting max likelihood?

During wondering about your questions and fiddling with some numerical experiments I came across with a peculiarity I didn't expect.

Here part of my R code tried to implement the FDA method.
It assumes that the data.frame PKdata has subject, sequence, code, repl and the PK-metrics AUC, Cmax. Key is here the variable repl which is the replicate number within subject and code (T or R).
avar contains the name of the PK metric (log-transformed).
# tmt2=T1, ... based on replicate number
PKdata$tmt2 <- paste(PKdata$code,PKdata$repl,sep="")
# reshape to subject, sequence, T1, ..., R1, ...
PKdata <- reshape(PKdata[c("subject","sequence","tmt2", avar)],
                  direction="wide", timevar="tmt2", v.names=avar,
                  idvar=c("subject","sequence"))
# change names as f.i. logAUC.T1 to T1
names(PKdata) <- sub(paste(avar,".",sep=""),"",names(PKdata))
# get the T1, T2, ... and R's correspondingly
Ts <- names(PKdata)[substring(names(PKdata),1,1)=="T"]
Rs <- names(PKdata)[substring(names(PKdata),1,1)=="R"]
# "basic estimator" aka within subject contrasts of T-R
# attention: coercion to a vector if only one T or R!
# rowMeans then does not work!
if (length(Ts)>1) Tv <- rowMeans(PKdata[Ts], na.rm=TRUE) else Tv <- PKdata[Ts]
if (length(Rs)>1) Rv <- rowMeans(PKdata[Rs], na.rm=TRUE) else Rv <- PKdata[Rs]
# attention: in case of only one T,R the column TR will be auto-named T1
# if assignment PKdata$TR <- is used!
PKdata["TR"] <- (Tv - Rv)
# now we can make the 'Proc GLM' call
# with standard contr.treatment we get not the wanted intercept!
# with that the intercept is intercept+seq1
oc <- options("contrasts")
options(contrasts=c("contr.sum","contr.poly"))
mlm <- lm(TR ~ sequence, data=PKdata)
# SAS gives us the same answer using Proc GLM with sequence as effect and ESTIMATE 'T-R' intercept 1;
CIs <- data.frame(point=coef(mlm)[1],
                  lower=confint(mlm,1,level=1-2*alpha)[1],
                  upper=confint(mlm,1,level=1-2*alpha)[2])
if (logtrans) CIs <- exp(CIs)
options(oc)  #restore
Note the attention: comments to see that R is also sometimes a hard to repressing beast.

Using the bear example dataset for replicate studies (4-period-2-sequence study with sequences RTRT/TRTR and 14 subjects balanced over sequence groups) I get:
      point      lower      upper
  1.01479     0.99336     1.03667    Intra-subject contrasts (MoM)
# compare to:
  1.01479     0.95223     1.08146    lme() like bear
  1.01479     0.95494     1.07839    SAS with FDA code Proc MIXED but CSH


Unbiased estimate of point estimator T-R but much tighter confidence interval of MoM!
Whats going on here? As I have commented in the above code "The power to know" gives us the same results. Thus the implementation seems error free. Someone out there to prove me wrong?

Then I remembered Stephen Senn's example 7.3 ("Cross-over Trials in Clinical Research" page 237 ff) who showed a similar effect due to a negative correlation of the two within-subject contrasts T1-R1 and T2-R2 which lead to a negative subject by formulation interaction if the constraints that variance-covariance components must >0 are relaxed.
In testing this for the bear data set I got a correlation coefficient r=-0.9366165, a very strong negative correlation, why ever.
Proc Mixed with type=UN as covariance matrix gives then a negative SxF interaction and the confidence interval:
      point      lower      upper
  1.01479     0.99336     1.03667    Proc Mixed, Type=UN
exactly the same as via the intra-subject contrasts.

BTW: Dear EM, just to cite you: "...you ask a question here in this forum to gain insight, but when you ask about X you get a story about Y ..."
Therefore here the results for the intra-individual variances (in the log domain):
         MoM     Proc Mixed FDA (REML)
s2WT    0.030130  0.02010 
s2WR    0.028125  0.01974