fixed & mixed (dammit!) and a request to SASians [Design Issues]

posted by Helmut Homepage – Vienna, Austria, 2019-03-23 01:03  – Posting: # 20070
Views: 3,359

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.

:lol3: 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=KENWARDROGER1 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)?

                                            PE        90% CI      CV%
pooled SAS Proc MIXED (ddfm=KENWARDROGER)  116.2   109.0  123.9   21.3 
2
       SAS Proc MIXED (ddfm=SATTERTHWAITE) 116.15  108.97 123.81  21.20 ╮
       PHX/WNL REML (DF Satterthwaite)     116.15  108.97 123.81  21.20 |
       R lmer() ddf="Satterthwaite"        116.15  108.97 123.81  21.20 |
       R lmer() ddf="Kenward-Roger"        116.15  108.97 123.81  21.20 ╯
       PHX/WNL fixed                       116.24  109.04 123.92  21.22 ╮
       R lm()                              116.24  109.04 123.92  21.22 ╯
IBD    SAS Proc MIXED (ddfm=SATTERTHWAITE) 116.05  108.92 123.65  20.84 ╮
       PHX/WNL REML (DF Satterthwaite)     116.05  108.92 123.65  20.84 ╯
       R lmer() ddf="Satterthwaite"        116.05  108.97 123.58  20.84 ╮
       R lmer() ddf="Kenward-Roger"        116.05  108.97 123.58  20.84 ╯
       PHX/WNL fixed                       116.14  108.98 123.78  20.88 ╮
       R lm()                              116.14  108.98 123.78  20.88 ╯


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=KENWARDROGER3 gives the same result as a fixed effects model. OK, but the FDA recommends Satterthwaite’s df for ages.4 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; 3rd 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.090, 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. online.

Cheers,
Helmut Schütz
[image]

The quality of responses received is directly proportional to the quality of the question asked. ☼
Science Quotes

Complete thread:

Activity
 Mix view
Bioequivalence and Bioavailability Forum |  Admin contact
19,780 posts in 4,196 threads, 1,357 registered users;
online 9 (0 registered, 9 guests [including 7 identified bots]).
Forum time (Europe/Vienna): 18:54 CEST

If you obey all the rules,
you will miss all the fun.    Katharine Hepburn

The BIOEQUIVALENCE / BIOAVAILABILITY FORUM is hosted by
BEBAC Ing. Helmut Schütz
HTML5