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

posted by Astea – Russia, 2016-11-06 12:50 (2699 d 06:55 ago) – Posting: # 16781
Views: 19,205

❝ 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:

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

Nothing shows a lack of mathematical education more
than an overly precise calculation.    Carl Friedrich Gauß

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