Bioequivalence and Bioavailability Forum

Hi sury,

an example in plain R. Once you have the data, three lines of code (one for the model and one each for the PE and CI):

alpha     <- 0.05
sequence  <- c("TR", "TR", "RT", "RT", "TR", "RT", "TR", "RT",
               "RT", "RT", "TR", "TR", "RT", "RT", "TR", "TR",
               "TR", "TR", "RT", "RT", "TR", "TR", "TR", "TR",
               "TR", "TR", "TR", "TR", "RT", "RT", "RT", "RT",
               "RT", "RT", "TR", "TR", "RT", "RT", "TR", "TR",
               "RT", "RT", "RT", "RT", "TR", "TR", "RT", "RT")
subject    <- rep(1:24, each = 2)
period     <- c(1, 2, 1, 2, 1, 1, 2, 2,
                1, 2, 1, 2, 1, 2, 1, 2,
                1, 2, 1, 2, 1, 2, 1, 2,
                1, 2, 1, 2, 1, 2, 1, 2,
                1, 2, 1, 2, 1, 2, 1, 2,
                1, 2, 1, 2, 1, 2, 1, 2)
treatment  <- c("T", "R", "R", "T", "T", "R", "R", "T",
                "R", "T", "T", "R", "R", "T", "T", "R",
                "T", "R", "R", "T", "T", "R", "T", "R",
                "T", "R", "T", "R", "R", "T", "R", "T",
                "R", "T", "T", "R", "R", "T", "T", "R",
                "R", "T", "R", "T", "T", "R", "R", "T")
PK         <- c( 53.45,  54.08, 131.39, 207.19, 185.92, 164.26,
                103.83,  79.81,  91.19, 103.40, 100.20,  46.17,
                117.75, 117.15,  51.68,  42.12, 149.37, 147.14,
                163.32, 157.48,  59.66,  59.28,  87.49, 101.02,
                 88.31,  80.70, 118.84, 123.49, 238.65, 124.01,
                 69.12,  68.02, 189.58, 127.24, 131.29, 101.05,
                208.16, 175.65, 104.44, 106.58, 126.55, 151.19,
                109.26,  91.66,  75.98,  99.08, 176.13, 119.58)
data       <- data.frame(sequence = sequence, subject = subject,
                         period = period, treatment = treatment, PK = PK)
cols       <- c("subject", "period", "sequence", "treatment")
data[cols] <- lapply(data[cols], factor) # factorize all
mod        <- lm(log(PK) ~ sequence + subject%in%sequence + period + treatment,
                           data = data)
pe         <- 100*exp(coef(mod)[["treatmentT"]])
ci         <- 100*exp(confint(mod, "treatmentT", level=1-2*alpha))
print(anova(mod), digits=7)
cat(paste0("\n",100*(1-2*alpha), "% CI:",
           sprintf("%6.2f%s%6.2f%%", ci[1], "\u2013", ci[2]),
           ", PE:", sprintf("%6.2f%%", pe), "\n")) # rounded acc. to GLs


Which gives

Analysis of Variance Table

Response: log(PK)
                 Df   Sum Sq   Mean Sq  F value     Pr(>F)   
sequence          1 1.874652 1.8746525 44.37735 1.7428e-06 ***
period            1 0.221688 0.2216881  5.24787 0.03296812 * 
treatment         1 0.002563 0.0025631  0.06067 0.80794404   
sequence:subject 24 5.343510 0.2226462  5.27055 0.00018651 ***
Residuals        20 0.844869 0.0422435                       
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

90% CI: 89.00–110.22%, PE: 99.04%

Note that this is the fixed-effects model recommended by the EMA. However, the nested structure subject(sequence) is superfluous and the model over-specified. Try print(mod) and you will see many lines which are NA (not specific to R, you get the same in any software). If we change the model to the simple one

mod <- lm(log(PK) ~ sequence + subject + period + treatment,
                    data = data)

we get

Analysis of Variance Table

Response: log(PK)
          Df   Sum Sq   Mean Sq  F value     Pr(>F)   
sequence   1 1.874652 1.8746525 46.41489 9.7698e-07 ***
subject   23 5.474793 0.2380345  5.89354 6.3723e-05 ***
period     1 0.087104 0.0871039  2.15662    0.15678   
treatment  1 0.002563 0.0025631  0.06346    0.80356   
Residuals 21 0.848170 0.0403890                       
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

90% CI: 89.00–110.22%, PE: 99.04%

Identical PE and CI. No more NAs if you try print(mod).

If you want subjects as a random effect and degrees of freedom like SAS Proc Mixed DDFM=CONTAIN:

library(nlme)
options(contrasts = c("contr.treatment", "contr.poly"))
mod <- lme(log(PK) ~ sequence + period + treatment, random = ~1|subject,
                     na.action = na.omit, data = data)
pe  <- 100*exp(summary(mod)$tTable["treatmentT", "Value"])
ci  <- 100*exp(intervals(mod, which = "fixed",
                         level = 1-2*alpha)[[1]]["treatmentT", c(1, 3)])

            numDF denDF  F-value p-value
(Intercept)     1    23 4306.009  <.0001
sequence        1    21    4.000  0.0586
period          1    21    2.280  0.1459
treatment       1    21    0.053  0.8196

90% CI: 88.38–109.88%, PE: 98.55%

Similar but Satterthwaite’s degrees of freedom like Proc Mixed DDFM=SATTERTHWAITE:

library(lmerTest)
mod <- lmer(log(PK) ~ sequence + period + treatment + (1|subject),
                      data = data)
pe  <- summary(mod, ddf = "Satterthwaite")$coefficients["treatmentT", "Estimate"]
ci  <- 100*exp(intervals(mod, which = "fixed",
                         level = 1-2*alpha)[[1]]["treatmentT", c(1, 3)])

Type III Analysis of Variance Table with Satterthwaite's method
            Sum Sq  Mean Sq NumDF  DenDF F value Pr(>F) 
sequence  0.192161 0.192161     1 37.903  4.0001 0.0527 .
period    0.109546 0.109546     1 19.844  2.2803 0.1468 
treatment 0.002563 0.002563     1 19.076  0.0534 0.8198 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

90% CI: 88.34–109.94%, PE: 98.55%

Always a good idea to learn the basics of R. However, consider Yung-jin’s nice package bear.