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

posted by Helmut Homepage – Vienna, Austria, 2019-08-05 16:23  – Posting: # 20474
Views: 882

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.

Cheers,
Helmut Schütz
[image]

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

Complete thread:

Activity
 Mix view
Bioequivalence and Bioavailability Forum |  Admin contact
19,983 posts in 4,228 threads, 1,373 registered users;
online 12 (1 registered, 11 guests [including 6 identified bots]).
Forum time (Europe/Vienna): 17:27 CET

I believe that a scientist looking at nonscientific problems
is just as dumb as the next guy.    Richard Feynman

The BIOEQUIVALENCE / BIOAVAILABILITY FORUM is hosted by
BEBAC Ing. Helmut Schütz
HTML5