permlme | R Documentation |
Permutation test for a single fixed effects of a linear mixed model.
This is essentially a re-implementation of the function predictmeans::perlmer
limited to fixed effects and random intercept models.
Should be faster than the original function and output the test statistic relative to each permutation.
permlme(
lme0,
lme1,
data = NULL,
seed = NULL,
statistic = "Wald",
perm.X = FALSE,
return.index.perm = FALSE,
nperm = 1000,
cpus = 1,
trace = TRUE
)
lme0 |
[lme] model under the null. |
lme1 |
[lme] model under the alternative. |
data |
[data.frame] dataset used to fit the models. Normally not need as the data is extracted from the models. |
seed |
[integer] specify a random number generator seed, for reproducible results. |
statistic |
[character] type of test statistic: Wald test and/or likelihood ratio test. |
perm.X |
[logical] should the covariates be also permuted? |
nperm |
[integer] number of permutations. |
cpus |
[cpus] number of cpus used to perform in parallel the permutations. |
trace |
[logical] should a progress bar displayed to follow the progress of the permutations? |
perm.index.perm |
[logical] should the index used to permute of the normalized residuals be output? |
a list containing
LRT: Likelihood ratio test for the coefficient. The p-value corresponding to the permutation test is denoted p-value.perm
.
Wald: Wald test for the coefficient. The p-value corresponding to the permutation test is denoted p-value.perm
.
n.perm: number of successful permutations
index.perm: index used to permute of the normalized residuals
Oliver E. Lee and Thomas M. Braun (2012), Permutation Tests for Random Effects in Linear Mixed Models. Biometrics, Journal 68(2).
n.perm <- 1000
## ** simulate data
library(lava)
library(data.table)
m <- lvm(Y1~1*eta+effect*conc,
Y2~1*eta+effect*conc,
Y3~1*eta+effect*conc,
Y4~1*eta+effect*conc,
Y5~1*eta+effect*conc)
latent(m) <- ~eta
transform(m, Id~eta) <- function(x,...){1:NROW(x)}
set.seed(1)
dtW <- data.table::data.table(lava::sim(m, p = c("effect" = 0.1), n = 50, latent = FALSE))
dtL <- data.table::melt(dtW, id.vars = c("Id","conc"))
setkeyv(dtL,"Id")
## ** "homemade" permutation test
library(nlme)
e.lmeH0 <- nlme::lme(value ~ variable, random =~1|Id, data = dtL)
e.lmeH1 <- nlme::lme(value ~ variable + conc, random=~1|Id, data = dtL)
e.permlme <- permlme(e.lmeH0, e.lmeH1, perm.X = TRUE, nperm = n.perm, seed = 10, statistic = "LRT")
e.permlme ## p=0.2197802
## ** permutation test via the function permlmer
library(lme4)
library(predictmeans)
e.lmerH0 <- lme4::lmer(value ~ variable + (1|Id), data = dtL, REML = FALSE)
e.lmerH1 <- lme4::lmer(value ~ variable + conc + (1|Id), data = dtL, REML = FALSE)
e.permlmer <- predictmeans::permlmer(e.lmerH0, e.lmerH1, nperm = n.perm, ncore = 1, seed = 10)
e.permlmer
## ** compare values
e.permlme$LRT[2,"p.value.perm"]-e.permlmer[["Perm-p"]][2] ## same!!!
## ** compare timing
library(microbenchmark)
microbenchmark(manual0 = permlme(e.lmeH0, e.lmeH1, nperm = 100, seed = 10, trace = FALSE),
manual = permlme(e.lmeH0, e.lmeH1, nperm = 100, seed = 10, trace = FALSE, statistic = c("LRT","Wald")),
package = predictmeans::permlmer(e.lmerH0, e.lmerH1, nperm = 100, ncore = 1, seed = 10),
times = 10)
## Unit: milliseconds (n=5)
## expr min lq mean median uq max neval cld
## manual0 689.9769 695.0086 752.4042 749.6891 770.365 918.2436 10 a
## manual 1963.7686 2053.5608 2111.0911 2135.7464 2147.678 2215.7608 10 b
## package 3433.9876 3456.7219 3591.3811 3641.0649 3662.663 3688.3372 10 c
## Unit: milliseconds
## expr min lq mean median uq max neval cld
## manual0 912.5592 939.0943 1117.806 1110.615 1278.572 1378.842 10 a
## manual 2646.7307 3291.1398 3368.645 3352.098 3466.784 3829.754 10 b
## package 4022.4287 4479.0607 4795.183 4530.408 4903.027 6810.617 10 c
## ** multiple testing
## *** unadjusted p-value
p.value <- e.permlme$LRT[["p.value.perm"]][2]
## *** perfect correlation between tests
e.permlme2 <- permlme(e.lmeH0, e.lmeH1, perm.X = TRUE, nperm = n.perm, seed = 10, statistic = "LRT")
e.permMax <- pmax(e.permlme$stat.perm[,"LRT"],e.permlme2$stat.perm[,"LRT"])
e.obsMax <- pmax(e.permlme$LRT[["L.Ratio"]][2],e.permlme2$LRT[["L.Ratio"]][2])
(sum(e.permMax > e.obsMax)+1)/(n.perm+1) - p.value ## should be 0
## *** perfect independence between tests
set.seed(11)
dtW2 <- data.table::data.table(lava::sim(m, p = c("effect" = 0.1), n = 50, latent = FALSE))
dtL2 <- data.table::melt(dtW2, id.vars = c("Id","conc"))
setkeyv(dtL2,"Id")
e.lmeH0.bis <- nlme::lme(value ~ variable, random =~1|Id, data = dtL2)
e.lmeH1.bis <- nlme::lme(value ~ variable + conc, random=~1|Id, data = dtL2)
e.permlme2 <- permlme(e.lmeH0.bis, e.lmeH1.bis, perm.X = TRUE, nperm = n.perm, seed = 10, statistic = "LRT")
p.value2 <- e.permlme2$LRT[["p.value.perm"]][2]
e.permMax <- pmax(e.permlme$stat.perm[,"LRT"],e.permlme2$stat.perm[,"LRT"])
e.obsMax <- pmax(e.permlme$LRT[["L.Ratio"]][2],e.permlme2$LRT[["L.Ratio"]][2])
(sum(e.permMax > e.obsMax)+1)/(n.perm+1) - 2*min(p.value ,p.value2) ## should be approximately true
## *** partially correlated tests
e.lmeH0.bis <- nlme::lme(value ~ variable, random =~1|Id, data = rbind(dtL,dtL2))
e.lmeH1.bis <- nlme::lme(value ~ variable + conc, random=~1|Id, data = rbind(dtL,dtL2))
e.permlme2 <- permlme(e.lmeH0.bis, e.lmeH1.bis, perm.X = TRUE, nperm = n.perm, seed = 10, statistic = "LRT")
p.value2 <- e.permlme2$LRT[["p.value.perm"]][2]
e.permMax <- pmax(e.permlme$stat.perm[,"LRT"],e.permlme2$stat.perm[,"LRT"])
e.obsMax <- pmax(e.permlme$LRT[["L.Ratio"]][2],e.permlme2$LRT[["L.Ratio"]][2])
(sum(e.permMax > e.obsMax)+1)/(n.perm+1)/min(p.value ,p.value2) ## should be between 1 and 2
## ** "homemade" permutation test in presence of missing values
set.seed(10)
dtL$valueMiss <- dtL$value
dtL$valueMiss[rbinom(NROW(dtL), size = 1, prob = 0.1)==1] <- NA
e.lmeH0 <- nlme::lme(valueMiss ~ variable, random =~1|Id, data = dtL, na.action = na.omit)
e.lmeH1 <- nlme::lme(valueMiss ~ variable + conc, random=~1|Id, data = dtL, na.action = na.omit)
e.permlme <- permlme(e.lmeH0, e.lmeH1, perm.X = TRUE, nperm = n.perm, seed = 10, statistic = "LRT")
e.permlme
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.