## Why not package bear? [R 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 Ing. Helmut Schütz 