Bioequivalence and Bioavailability Forum

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

posted by mittyri  – Russia, 2016-11-07 07:07 (3147 d 12:37 ago) – Posting: # 16782
Views: 21,345

Hi Astea,

❝ Here is my rough code:

Nice work!
Here's mine. I tried to prepare the dataframe for validation (Results)

library(readxl)
library(lmerTest)
library(dplyr)
options(contrasts=c("contr.treatment", "contr.poly"))
methodB <- function(Dataset){
  muddle.lmer <- (lmer(log(Data)~Sequence+Period+Formulation + (1|Subject), data=Dataset))
  cat("\n", "lmer treats Subject, Sequence, Period and Formulation as factors, Random is /~1|Subject/ (Method B), then")
  FormulationT.lmer <- summary(muddle.lmer)$coefficients["FormulationT",]
  Results<-data.frame(exp(as.numeric(FormulationT.lmer[1])),
    exp(FormulationT.lmer["Estimate"] + FormulationT.lmer["Std. Error"]*qt(0.05, FormulationT.lmer["df"])),
    exp(FormulationT.lmer["Estimate"] - FormulationT.lmer["Std. Error"]*qt(0.05, FormulationT.lmer["df"])))
  colnames(Results) <- c("PE", "Lower_CI", "Upper_CI")
  cat("\n", "PE is", sprintf("%.2f", Results$PE*100), "%\n")
  cat("CI is", sprintf("%.2f", Results$Lower_CI*100), "-", sprintf("%.2f", Results$Upper_CI*100), "%\n")
  # get reference dataset ready
  Ref<-Dataset %>%
    filter(Formulation=="R") %>%
    group_by(Subject) %>%
    filter(max(n())==2)
  Ref.lm <- lm(log(Data)~Sequence+Period+Subject, data=Ref)
  MSerror <-sum(Ref.lm$residuals^2)/Ref.lm$df.residual
  # get Stud Residuals ready
  Ref$StudResiduals <- rstandard(Ref.lm)
  StudResidualsFirst <- Ref %>%
    filter(Period=="1"|Period=="2")
  # plotting
  bxpdat <- boxplot(StudResidualsFirst$StudResiduals, range = 3, main ="Studentized Residuals Boxplot, IQR = 3")
  text(bxpdat$group,                                              # the x locations
       bxpdat$out,                                                # the y values
       StudResidualsFirst[match(StudResidualsFirst$StudResiduals, bxpdat$out, nomatch=0)>0, "Subject"]$Subject,  # the labels
       pos = 4) 
  qqnorm(StudResiduals)
  # estimating CVr
  Results$CV_R <- (exp(MSerror)-1)^(1/2)*100
  cat("CV_R is ", sprintf("%.2f",Results$CV_R), "%\n")
  # decision tree for CV_Rs
  if(Results$CV_R>30)  {
    if(Results$CV_R<=50)    {
      cat("certainly we need scaling\n")
      Results$Ext_L_Bound<-exp(-MSerror^(1/2)*0.760)
      Results$Ext_U_Bound<-exp(+MSerror^(1/2)*0.760) 
      cat("the scaled bounds are", sprintf("%.2f", Results$Ext_L_Bound*100), "-", sprintf("%.2f", Results$Ext_U_Bound*100), "%\n" )
    } else  {
      cat("certainly we need scaling, but the bounds should be fixed\n")
      Results$Ext_L_Bound<-0.6984
      Results$Ext_U_Bound<-1.4319
      cat("the scaled bounds are", sprintf("%.2f", Results$Ext_L_Bound*100), "-", sprintf("%.2f", Results$Ext_U_Bound*100), "%\n" )
    }
  }  else  {
    cat("we do not need scaling (bounds are 80.00-125.00 %)\n")
    Results$Ext_L_Bound<-80.00
    Results$Ext_U_Bound<-125.00
  }
  # print results
  if(Results$Lower_CI>Results$Ext_L_Bound & Results$Upper_CI<Results$Ext_U_Bound) {
    cat("BE is met, congrats")
  } else {
    cat("BE is not met, sorry")
  }
  return(Results)
}

# downloading the reference file from BEBAC; that's the most convenient way to get the file
download.file('http://bebac.at/downloads/Validation Replicate Design EMA.xls', "Dataset.xls", cacheOK = FALSE, mode="wb")
# read it without Perl and Java; thanks, Hadley!
Dataset<-read_excel("Dataset.xls", sheet = 1)
# prepare dataset
Dataset$Formulation<-factor(Dataset$Formulation, levels = c("R", "T"))
Dataset$Sequence<-factor(Dataset$Sequence, levels = c("TRTR", "RTRT"))
Dataset$Period<-factor(Dataset$Period)
Dataset$Subject<-factor(Dataset$Subject)
Results<- methodB(Dataset)

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

I suspect deviations from SAS, because some values are not equal in full precision...
OK, that would be my next task

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 scaling Astea 2016-11-06 11:50
                  • lmer: Method B is ready for scalingmittyri 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