sury
☆    

India,
2019-08-05 12:30

Posting: # 20472
Views: 1,006
 

 R-Code for 2×2 Bio­equi­valence study [R for BE/BA]

Hii all,

I am very beginner for the R-Coding and have searched for the code for 2X2 Bio-equivalence code in the forum

Can anyone help me in this regards as i got to know that lm () do the same thing as that of the PROC GLM in the SAS but i am nt aware of how to use the function for my bio equivalence study. As i cant afford the licenced SAS Software, i want to run my analysis in r.

Any help in this regards is highly appreciated.
yjlee168
★★  
avatar
Homepage
Kaohsiung, Taiwan,
2019-08-05 14:23

@ sury
Posting: # 20473
Views: 903
 

 R-Code for 2×2 Bioequivalence study

Hi sury,

If you do "Search...", you can find this gem. The codes are included in that .rtf file.

» ... As i cant afford the licenced SAS Software, i want to run my analysis in r.

All the best,
-- Yung-jin Lee
bear v2.8.6:- created by Hsin-ya Lee & Yung-jin Lee
Kaohsiung, Taiwan http://pkpd.kmu.edu.tw/bear
Download link (updated) -> here
sury
☆    

India,
2019-08-06 11:35

@ yjlee168
Posting: # 20477
Views: 820
 

 R-Code for 2×2 Bioequivalence study

» If you do "Search...", you can find this gem. The codes are included in that .rtf file.

Ya got it, thanks for the help.:-)

» » ... As i cant afford the licenced SAS Software, i want to run my analysis in r.
Helmut
★★★
avatar
Homepage
Vienna, Austria,
2019-08-05 16:23

@ sury
Posting: # 20474
Views: 884
 

 Why not package bear?

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
sury
☆    

India,
2019-08-06 11:33

@ Helmut
Posting: # 20476
Views: 822
 

 Why not package bear?

Thanks helmut, tried the above code and it is matching with the SAS Results.

» Always a good idea to learn the basics of R.

Ya, You are right and working on it to improve my r basics, but however i have searched for many websites to learn the Bio equivalence packages, but not able to find how can i implement it.

» However, consider Yung-jin’s nice package bear.

Yes, i heard about the BEAR Package, and also installed the same package. when i am trying to use the package, it showing the warning message "package ‘bear’ was built under R version 3.5.2". i dont know how to clear this error and how can i call the package for using in the Bio equivalence study. Any help in this regards is highly appreciated.

Also want to know about the PKNCA package on R. How can i use the package. Does the package covers all the parameters as that of the WinNonlin. Does it gives the same result as that in the Phoenix WinNonlin Software

Thanks in Advance..


Edit: Full quote removed. Please delete everything from the text of the original poster which is not necessary in understanding your answer; see also this post #5[Helmut]
ElMaestro
★★★

Belgium?,
2019-08-06 12:10

@ sury
Posting: # 20478
Views: 804
 

 Why not package bear?

Hi sury,

» Yes, i heard about the BEAR Package, and also installed the same package. when i am trying to use the package, it showing the warning message "package ‘bear’ was built under R version 3.5.2 ". i dont know how to clear this error and how can i call the package for using in the Bio equivalence study. Any help in this regards is highly appreciated.

It sounds like you need to study some basics of R: How your GUI or commandline works, what packages are, how to use a package (notably also learning to decipher the help files that are a standard feature of packages), what the difference is between a warning and an error, and so forth. It is not too hard if you have a little idea about coding, but it may indeed take a lot of time.
If you are serious, that time is well invested.

I could be wrong, but...
Best regards,
ElMaestro
yjlee168
★★  
avatar
Homepage
Kaohsiung, Taiwan,
2019-08-06 12:31

@ sury
Posting: # 20479
Views: 801
 

 how to build bear yourself from source tarball

Hi sury,

It still works like a charm, though the current bear release was built under R v3.5.2 (I guess your platform should be macOS). If you like to build it yourself for current R v3.6.1 from source, you can download the source tarball from Sourceforge, and place the tarball file under Home directory. Type R CMD INSTALL bear_2.8.4.tar.gz under a terminal. That'll be all. Good luck.

» ... it showing the warning message "package ‘bear’ was built under R version 3.5.2 ". i dont know how to clear this error and how can i call the package for using in the Bio equivalence study. Any help in this regards is highly appreciated.

All the best,
-- Yung-jin Lee
bear v2.8.6:- created by Hsin-ya Lee & Yung-jin Lee
Kaohsiung, Taiwan http://pkpd.kmu.edu.tw/bear
Download link (updated) -> here
Helmut
★★★
avatar
Homepage
Vienna, Austria,
2019-08-06 12:35

@ sury
Posting: # 20480
Views: 794
 

 R: warning ≠ error

Hi sury,

» Yes, i heard about the BEAR Package, and also installed the same package. when i am trying to use the package, it showing the warning message "package ‘bear’ was built under R version 3.5.2". i dont know how to clear this error

This is a warning, not an error. There is no problem using any package in a higher version of R. If a package tries to use a function which was deprecated* in your R-version, you will get a specific warning pointing to the function and generally an information which function should be used instead or how the call should be modified. Then you have two options:
  1. If you are an expert, modify the source code and build the package locally.
  2. Ask the package’s maintainer for an update. Probably the best cause sooner or later this warning will turn into an error.
Only if a package tries to use a function which was removed in a higher version of R, you will get an error. If the two options above are not possible, install an earlier version of R from CRAN (Windows binaries, sources). You can have any number of versions of R simultaneously installed on your machine. Use this as a last resort because generally binary data files (.rda) are associated with the last version you installed.

» … and how can i call the package for using in the Bio equivalence study. Any help in this regards is highly appreciated.

It’s straightforward. Though Yung-jin can explain it better than I could. ;-)

» Also want to know about the PKNCA package on R. How can i use the package.

PKNCA is a nice one as well. Why don’t you consult the man-pages as usual?

library(PKNCA)
help(package=PKNCA)


» Does the package covers all the parameters as that of the WinNonlin. Does it gives the same result as that in the Phoenix WinNonlin Software

AFAIK, yes. Once you started the main man-pages (see above) click
  • User guides, package vignettes and other documentation.


  • Deprecated functions are still working in the current version of R but will be defunct (removed) in the next version. All changes back to 3.0.0 here.

Cheers,
Helmut Schütz
[image]

The quality of responses received is directly proportional to the quality of the question asked. ☼
Science Quotes
Activity
 Thread view
Bioequivalence and Bioavailability Forum |  Admin contact
19,983 posts in 4,228 threads, 1,373 registered users;
online 14 (1 registered, 13 guests [including 10 identified bots]).
Forum time (Europe/Vienna): 17:31 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