Bioequivalence and Bioavailability Forum • lmer: Method B is ready for scaling

lmer: Method B is ready for scaling [🇷 for BE/BA]

posted by Astea – Russia, 2016-11-06 12:50 (3092 d 02:39 ago) – Posting: # 16781
Views: 20,802

❝ If nobody does it during the next week I'll try to code it by myself.

Challenge accepted!

Here is my rough code:

library(readxl)

library(lmerTest)

options(contrasts=c("contr.treatment","contr.poly"))

Dataset<-read_excel("Dataset.xls", sheet = 1)



methodB.scABEL <- function(Dataset){

Dataset$Formulation<-factor(Dataset$Formulation, levels = c("R", "T")) 

Dataset$Sequence<-factor(Dataset$Sequence, levels = c("TRTR", "RTRT")) 

Dataset$Period<-factor(Dataset$Period)

mod.lmer <- (lmer(log(Data)~Sequence+Period+Formulation + (1|Subject/Sequence), data=Dataset))

FormulationEffect.lmer <- summary(mod.lmer)$coefficients["FormulationT","Estimate"]
  cat("\n", "PE is", round(exp(as.numeric(FormulationEffect.lmer))*100, digits = 2), "\n")

FormulationT<-summary(mod.lmer)$coefficients["FormulationT",]
 

CI_Lr<-round(exp(FormulationT[1]+FormulationT[2]*qt(0.05, FormulationT[3]))*100, 2)

CI_Ur<-round(exp(FormulationT[1]-FormulationT[2]*qt(0.05, FormulationT[3]))*100, 2)

cat("\n", "CI is", CI_Lr, "-", CI_Ur,"\n")



Dataset2<-Dataset[order(Dataset$Formulation), ]


n<-length(which(Dataset2$Formulation=="R"))



Dataset3<-Dataset2[-c((n+1):nrow(Dataset2)),]


for (i in 1:nrow(Dataset3)) {if(Dataset3$Period[i]==3) Dataset3$Formulation[i]<-"T" else if(Dataset3$Period[i]==2) Dataset3$Formulation[i]<-"T"}

for (i in 1:nrow(Dataset3)) {if(Dataset3$Period[i]==3) Dataset3$Period[i]<-2 else if(Dataset3$Period[i]==4) Dataset3$Period[i]<-2 else if(Dataset3$Period[i]==2) Dataset3$Period[i]<-1}



Dataset3$Formulation<-factor(Dataset3$Formulation, levels = c("R", "T")) 

Dataset3$Sequence<-factor(Dataset3$Sequence, levels = c("RTRT", "TRTR")) 

Dataset3$Period<-factor(Dataset3$Period)

Dataset3$Subject<-factor(Dataset3$Subject)

muddle   <- lm(log(Data)~Formulation+Period+Sequence+Subject, data=Dataset3)

CV_R<-(exp(anova(muddle)[5,3])-1)^(1/2)*100

cat("\n", "CV_R is", round(CV_R, 2),"\n")

{if (CV_R>30)

 CI_L<-exp(-anova(muddle)[5,3]^(1/2)*0.760)*100

 CI_U<-exp(anova(muddle)[5,3]^(1/2)*0.760)*100

cat("Wide CI:", round(CI_L,2),"-", round(CI_U,2))}

{if ((CI_Lr<CI_L)|(CI_Ur>CI_U)) cat("BE failed")

else cat("=> BE success")}

}



methodB.scABEL(Dataset)

gives

 PE is 115.73 



 CI is 107.17 - 124.97 



 CV_R is 46.96 

Wide CI: 71.23 - 140.4=> BE success>

as in Phoenix (see Helmut's RSABE tutorial). The code needs to be tested via another data...

—
"Being in minority, even a minority of one, did not make you mad"

Complete thread:

Bear vs. Phoenix & SAS Helmut 2015-04-20 17:34
- R vs. Phoenix & SAS? yjlee168 2015-04-20 19:36
  - R vs. Phoenix & SAS? Helmut 2015-04-21 01:02
    - lme() does not work with all fixed effects yjlee168 2015-04-21 23:41
      - lme() does not work with all fixed effects Astea 2016-11-04 00:13
        
        lmer: Method B (PE catched for imbalanced dataset!!!) and Method C mittyri 2016-11-05 17:38
        
        lmer: Method B (PE catched for imbalanced dataset!!!) and Method C Astea 2016-11-05 19:27
        
        lmer: Method B is ready for scaling mittyri 2016-11-05 20:01
        
        lmer: Method B is ready for scalingAstea 2016-11-06 11:50
        
        lmer: Method B is ready for scaling mittyri 2016-11-07 06:07
- info for lsmeans yjlee168 2015-04-20 21:34
  - info for lsmeans Helmut 2015-04-21 01:15
    - once more about R and replicate designes Astea 2016-11-02 23:43
      - once more about R and replicate designes VStus 2016-11-06 11:34
        
        Getting variance components from the lmer output StatR 2017-02-03 13:53
        
        Getting variance components from the lmer output VStus 2017-02-03 15:47
        
        Getting variance components from the lmer output StatR 2017-02-03 17:12
        
        Getting variance components d_labes 2017-02-07 11:16
        
        Getting variance components StatR 2017-02-07 11:36
        Getting variance components StatR 2017-02-08 08:41
        
        Getting variance components d_labes 2017-02-08 10:13
        
        Getting variance components StatR 2017-02-08 10:19
        
        Data structure Helmut 2017-02-08 10:33
        
        Data structure StatR 2017-02-08 10:49