lmer: Method B is ready for scaling [🇷 for BE/BA]
❝ 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"
"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
- lmer: Method B is ready for scalingAstea 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