## Why not package bear? [🇷 for BE/BA]

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, "\u2013", ci),            ", 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)[]["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)[]["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.

Dif-tor heh smusma 🖖 Helmut Schütz The quality of responses received is directly proportional to the quality of the question asked. 🚮
Science Quotes 22,110 posts in 4,630 threads, 1,567 registered users;
online 5 (0 registered, 5 guests [including 4 identified bots]).
Forum time: Friday 20:05 CEST (Europe/Vienna)

We absolutely must leave room for doubt
or there is no progress and no learning.
There is no learning without having to pose a question.
And a question requires doubt.
People search for certainty.
But there is no certainty.    Richard Feynman

The Bioequivalence and Bioavailability Forum is hosted by Ing. Helmut Schütz 