Bioequivalence and Bioavailability Forum

Hi ElAlumno,

❝ I agree with the EMA Biostatistics Working Party […]. 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.

It is really existing though David Brown will be brexited (is this already a word?) unceremoniously.

❝ But I digress.

That’s fine. Adds some spice to that dry stuff.

❝ […] I tried analyzing Patterson & Jones' Example 4.5 with the Two‐at‐a‐Time Principle using R 1.1.456.

What? 1.1.456? 1.1.4 was the first one I used in 2002. Current is 3.5.3!

❝ 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,…

That’s not the point. The devil is in the details. I wrote above:

❝ ❝ ❝ ❝ Phoenix/WinNonlin 8, mixed effects

Why? Well, Patterson & Jones analyzed the data with SAS Proc MIXED and I wanted to reproduce it in PHX/WNL (restricted maximum likelihood). OK, they used in SAS ddfm=KENWARDROGER¹ and yours truly Satterthwaite’s df.
When I switched in PHX/WNL to a fixed effects model (i.e., like lm() in R and SAS Proc GLM) I got exactly your results. Hence, your code is fine though you can streamline it (specify NAs and classes):

dta <- read.table("exam45.dat", header=TRUE, na.strings="99999",
                  colClasses=c(rep("factor", 4), rep("numeric", 2)))

Checks:

sum(is.na(dta$AUC))
[1] 6
sum(is.na(dta$CMAX))
[1] 2
str(dta)
'data.frame':   186 obs. of  6 variables:
 $ subject : Factor w/ 62 levels "001","002","004",..: […]
 $ sequence: Factor w/ 6 levels "1","2","3","4",..: […]
 $ period  : Factor w/ 3 levels "1","2","3": […]
 $ formula : Factor w/ 3 levels "R","S","T": […]
 $ AUC     : num  4089 7411 5513 2077 3684 ...
 $ CMAX    : num  906 1711 1510 504 845 ...

I would not use option(..., digits=4) cause IMHO, this may affect also the precision of internal calculations. Better to round yourself at the end, f.i.

IBD[1, 2:5] <- round(c(foo, bar, baz), 2).

So what do we have (T/R of AUC)?
$$\small{\begin{array}{llccc}\hline
\text{Method} & \text{Setup} & \text{PE} & \text{90% CI} & \text{CV (%)} \\
\hline
\text{pooled} & \text{SAS Proc MIXED (KENWARDROGER)} & \text{116.2 } & \text{109.0 , 123.9 } & \text{21.3}^2 \\
& \text{SAS Proc MIXED (SATTERTHWAITE)} & 116.15 & \text{108.97, 123.81} & 21.20 \\
& \text{PHX/WNL REML (Satterthwaite)} & 116.15 & \text{108.97, 123.81} & 21.20 \\
& \text{R lmer() Satterthwaite} & 116.15 & \text{108.97, 123.81} & 21.20 \\
& \text{R lmer() Kenward-Roger} & 116.15 & \text{108.97, 123.81} & 21.20 \\\hline
& \text{PHX/WNL LM (all fixed)} & 116.24 & \text{109.04, 123.92} & 21.22 \\
& \text{R lm()} & 116.24 & \text{109.04, 123.92} & 21.22 \\
\hline
\text{IBD} & \text{SAS Proc MIXED (SATTERTHWAITE)} & 116.05 & \text{108.92, 123.65} & 20.84 \\
& \text{PHX/WNL REML (Satterthwaite)} & 116.05 & \text{108.92, 123.65} & 20.84 \\\hline
& \text{R lmer() Satterthwaite} & 116.05 & \text{108.97, 123.58} & 20.84 \\
& \text{R lmer() Kenward-Roger} & 116.05 & \text{108.97, 123.58} & 20.84 \\\hline
& \text{PHX/WNL LM (all fixed)} & 116.14 & \text{108.98, 123.78} & 20.88 \\
& \text{R lm()} & 116.14 & \text{108.98, 123.78} & 20.88 \\
\hline
\end{array}}$$If Satterthwaite’s df are used, the mixed effects model in SAS and PHX/WNL agree. R lm() agrees with PHX/WNL’s fixed effects model. A little bit surprising to me that Proc MIXED with ddfm=KENWARDROGER³ gives the same result as a fixed effects model. OK, but the FDA recommends Satterthwaite’s df for ages.⁴ A mixed effects model may give a lower residual variance – like in this case – because information of incomplete data is recovered. R lmer() is set up according to the EMA’s ‘Method B’ which treats subject as a random effect but is not a true mixed effects model. Given only for completeness (for the code see below).
A request to the SASians of the forum: Please run Proc MIXED with DDFM=SATTERTHWAITE of the pooled data and both DDFM-options on the data where S is excluded. THX to mittyri!

Mixed effects in R are a swamp. The closest I could get was this:

library(nlme)
muddle1 <- lme(log(AUC) ~ Trt + Per + Seq, data = dta[dta$Trt != "S", ],
                          random = ~Trt-1|Subj, weights = varIdent(form = ~1|Trt),
                          method = "REML", na.action = na.exclude,
                          control = lmeControl(opt = "optim", msMaxIter = 200))
anova(muddle1, type="marginal")
summary(muddle1)

Extracting the required stuff is beyond me. We had many discussions already in the forum and the verdict was “mixed effects models for BE are not possible in R”.
So if you are happy with a fixed effect model (many jurisdictions like the EMA), fine. If you have to deal with the FDA or Health Canada, maybe not.

---

1. Jones B, Kenward MG. Design and Analysis of Cross-Over Trials. Boca Raton: CRC Press; 3^rd ed. 2015.
   
2. Limited precision of results in the log-domain of the reference, CV not given. Estimated:
   library(PowerTOST)
round(100*CVfromCI(lower = 1.09, upper = 1.239, design = "3x6x3",
                   n = c(8, 11, 11, 10, 10, 10)), 1)
   
3. SAS tells me: This method changes output in the following tables: Contrast, CorrB, CovB, Diffs, Estimates, InvCovB, LSMeans, Slices, SolutionF, SolutionR, Tests1–Tests3. The OUTP= and OUTPM= data sets are also affected.
   
4. FDA/CDER. Guidance for Industry. Statistical Approaches to Establishing Bioequivalence. January 2001. download.