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

posted by mittyri – Russia, 2016-11-07 07:07 (2699 d 00:50 ago) – Posting: # 16782
Views: 19,323

Hi Astea,

❝ Here is my rough code:


Nice work! :ok:
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 :-D

Kind regards,
Mittyri

Complete thread:

UA Flag
Activity
 Admin contact
22,957 posts in 4,819 threads, 1,636 registered users;
111 visitors (0 registered, 111 guests [including 4 identified bots]).
Forum time: 07:58 CET (Europe/Vienna)

With four parameters I can fit an elephant,
and with five I can make him wiggle his trunk.    John von Neumann

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