Supplementary Material #3: simulation results for linear mixed-effects modeling
We used the following linear mixed-effect model in simulation to investigate the type I error rate and power performance of likelihood ratio test (LRT) and EBEs-based test in the context of linear mixed-effect modeling,
For null hypothesis , the likelihood ratio test is constructed as the -2log of likelihood of the models withand. EBEs-based test is constructed by two steps, that is, the empirical Bayes predictors of are firstly obtained from the base model without covariate and then the simple linear regression is applied using the EBEs as the dependent variable and centered lWTi as the independent variable. In the simulation studies, sample size N=20. Measurement times per individual nt is 2, 3 and 12 respectively. For nt=2, the measuring time points are 0.05 and 1; For nt=3, the time points are 0.05, 0.3, and 1 and for nt=12, time points are 0.01, 0.05, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, and 1. The covariate effect is 0, 0.15 and 0.2 and and For each scenario, 10,000 data sets are generated. We use the lme function in R to do estimation and testing for LRT and to generate the EBEs. The t-test statistics for H0 based on EBEs are reported. The percentage that the chi-square or t statistics is greater than the corresponding threshold is reported as the empirical type I error and power. The R code is attached bellow for reference.
# linear model for mean profile
flin <- function(base, t, SL) base + SL * t
#simuation function
foo.sim.est <- function(N, nt, WTP, sd.ETA, sd.Sigma) {
if(nt==2) time.v=c(0.05, 1)
if(nt==3) time.v=c(0.05, 0.3, 1)
if(nt==12) time.v=c(0.01, 0.05, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1)
ID <- rep(seq(1:N),each=nt)
Time <- rep(time.v,N)
YIJ <- rep(0,N*nt)
lWT.v <- rnorm(N, log(70), 0.3)
lWT <- rep(lWT.v, each=nt)
BSSim <- rep(rep(0,N)+sd.ETA*rnorm(N),each=nt)
SLSim <- rep(rep(0.5,N)+sd.ETA*rnorm(N)+WTP*(lWT.v-log(70)),each=nt)
myData <- as.data.frame(list(ID=ID,Time=Time,lWT=lWT, YIJ=YIJ))
myData$IPRE <- flin(t=myData$Time, base=BSSim, SL=SLSim)
myData$YIJ <- myData$IPRE + sd.Sigma*rnorm(dim(myData)[1])
myData <- groupedData( YIJ ~ Time | ID,
data = myData,
labels = list( x = "Time", y = "observation"))
fit1 <- try(lme(YIJ ~ Time, data = myData, random = pdDiag(~Time), method="ML"), silent=TRUE)
fit2 <- try(lme(YIJ ~ Time+ Time:I(lWT-log(70)), data = myData, random = pdDiag(~Time), method="ML"), silent=TRUE)
if(inherits(fit1,"try-error")|inherits(fit2,"try-error"))
{
output <- c(rep(NA, 6))
}
else {
ebe.parm <- ranef( fit1, augFrame = T, data=myData)
ebe.t <- summary(lm(Time ~ I(lWT-log(70)), data=ebe.parm))$coef[2, c(1, 4)]
#####################################
# shrinkage for CL
sh <- 1-var(ebe.parm$Time)/as.numeric(VarCorr(fit1)[2,1])#(sd.ETA^2)
# log likelihood ratio test
mod.lrt <- 1-pchisq(-2*fit1$logLik - -2*fit2$logLik, 1)
output <- c(ebe.t, mod.lrt, sh)
}
return(output)
1
Page 3 of 3