permlme: Permutation test for fixed effects of a linear mixed model.
In bozenne/butils: Various useful functions
permlme R Documentation
Permutation test for fixed effects of a linear mixed model.
Description
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.
Usage
permlme(
lme0,
lme1,
data = NULL,
seed = NULL,
statistic = "Wald",
perm.X = FALSE,
return.index.perm = FALSE,
nperm = 1000,
cpus = 1,
trace = TRUE
)
Arguments
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?
Value
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
References
Oliver E. Lee and Thomas M. Braun (2012), Permutation Tests for Random Effects in Linear Mixed Models. Biometrics, Journal 68(2).
Examples
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
bozenne/butils documentation built on Oct. 14, 2023, 6:19 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
| permlme | R Documentation |
Permutation test for fixed effects of a linear mixed model.
Description
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.
Usage
permlme(
lme0,
lme1,
data = NULL,
seed = NULL,
statistic = "Wald",
perm.X = FALSE,
return.index.perm = FALSE,
nperm = 1000,
cpus = 1,
trace = TRUE
)
Arguments
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? |
Value
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
References
Oliver E. Lee and Thomas M. Braun (2012), Permutation Tests for Random Effects in Linear Mixed Models. Biometrics, Journal 68(2).
Examples
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
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.