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

posted by Helmut Homepage – Vienna, Austria, 2019-08-05 18:23 (1919 d 23:06 ago) – Posting: # 20474
Views: 7,558

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.

Dif-tor heh smusma 🖖🏼 Довге життя Україна! [image]
Helmut Schütz
[image]

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

Complete thread:

UA Flag
Activity
 Admin contact
23,288 posts in 4,890 threads, 1,658 registered users;
21 visitors (0 registered, 21 guests [including 9 identified bots]).
Forum time: 16:29 CET (Europe/Vienna)

‘Normal science’ means research firmly based upon one or more
past scientific achievements, achievements that some particular
scientific community acknowledges for a time as supplying
the foundation for its further practice.    Thomas S. Kuhn

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