Bioequivalence and Bioavailability Forum

Dear ElMaestro!

❝ ❝ For Example 2 I got:

❝ Where is example 2??

At the end of Madame Potvin’s paper. For convenience:
sub  <- as.factor(sort(rep(c(1:20),2)))
stg  <- as.factor(c(rep(1,12),rep(2,8)))
seq  <- as.factor(c(rep('TR',6),c(rep('RT',6)),
                    rep('TR',4),c(rep('RT',4))))
per  <- as.factor(rep(c(1,2),20))
trt  <- as.factor(c(rep(c('T','R'),6),rep(c('R','T'),6),
                    rep(c('T','R'),4),rep(c('R','T'),4)))
data <- data.frame(stg, sub, trt, per, seq)
data$resp <- c(4.8,4.6,1.4,1.3,3.9,4.4,1.5,2.4,3.8,2.9,3.3,2.5,
               1.5,1.7,1.9,1.3,2.4,2.1,13.9,11.4,2.3,1.5,2.4,2.6,
               4.0,4.2,0.7,0.5,10.3,11.9,4.7,4.8,5.5,5.9,5.4,3.8,2.1,4.3,2.4,3.6)

❝ This should absolutely work, I think. Either it is the syntax that was wrong or PHX/WNL is flawed?

Well, the flaw sits behind the LSMeans.

❝ You don't need a mixed model when you have only one sigma in play (this isn't a replicated study I presume?).

I know the Silly-O-Meter! In PHX/WNL everything goes the REML-way. The only way to mimick the behaviour of SAS Proc GLM is to specify all effects fixed, and filter incomplete data beforehand. Here this is not an issue, of course.

❝ In the simple non-replicated case fitting a mixed model with subject as random should give the same result as the normal linear model based on all fixed.

Yes.

❝ If there are subtle differences at some decimal then this sounds most likely just like either a rounding phenomenon, or a convergence issue where somehow the optimiser in this case is capable of finding smaller sigma when the mixed model is applied.

Yep – though the setup (convergence criteria, step size, ) was the same. Doesn’t bother me too much (yet).

❝ The normal linear model can be solved exactly with e.g. R at least to the extent we can rely on its numerical stability for inversion of X^tX; I would take a look at sigma that way and see what the truth might look like.

I’m try to dig out some code from the grave.