lmer: Method B is ready for scaling [🇷 for BE/BA]
Hi Astea,
Nice work!
Here's mine. I tried to prepare the dataframe for validation (Results)
I suspect deviations from SAS, because some values are not equal in full precision...
OK, that would be my next task
❝ 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

—
Kind regards,
Mittyri
Kind regards,
Mittyri
Complete thread:
- Bear vs. Phoenix & SAS Helmut 2015-04-20 17:34 [🇷 for BE/BA]
- 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
- lmer: Method B is ready for scaling Astea 2016-11-06 11:50
- lmer: Method B is ready for scaling mittyri 2016-11-05 20:01
- lmer: Method B (PE catched for imbalanced dataset!!!) and Method C Astea 2016-11-05 19:27
- lmer: Method B (PE catched for imbalanced dataset!!!) and Method C mittyri 2016-11-05 17:38
- lme() does not work with all fixed effects Astea 2016-11-04 00:13
- lme() does not work with all fixed effects yjlee168 2015-04-21 23:41
- R vs. Phoenix & SAS? Helmut 2015-04-21 01:02
- 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
- Getting variance components d_labes 2017-02-08 10:13
- Getting variance components d_labes 2017-02-07 11:16
- Getting variance components from the lmer output StatR 2017-02-03 17:12
- Getting variance components from the lmer output VStus 2017-02-03 15:47
- Getting variance components from the lmer output StatR 2017-02-03 13:53
- once more about R and replicate designes VStus 2016-11-06 11:34
- once more about R and replicate designes Astea 2016-11-02 23:43
- info for lsmeans Helmut 2015-04-21 01:15
- R vs. Phoenix & SAS? yjlee168 2015-04-20 19:36