Title: | Power Analysis for Random Effects in Mixed Models |
---|---|
Description: | Simulation functions to assess or explore the power of a dataset to estimates significant random effects (intercept or slope) in a mixed model. The functions are based on the "lme4" and "lmerTest" packages. |
Authors: | Julien Martin [aut, cre] |
Maintainer: | Julien Martin <[email protected]> |
License: | GPL-2 | GPL-3 |
Version: | 1.122 |
Built: | 2024-11-20 03:34:22 UTC |
Source: | https://github.com/juliengamartin/pamm_r |
Simulation functions to assess or explore the power of a dataset to estimates significant random effects (intercept or slope) in a mixed model. The functions are based on the "lme4" and "lmerTest" packages.
Maintainer: Julien Martin [email protected]
Useful links:
Report bugs at https://github.com/JulienGAMartin/pamm_R/issues
Given a specific sample size, fixed number of group and replicates per group, the function simulate different variance-covariance structure and assess p-values and power of random intercept and random slope
EAMM( numsim, group, repl, fixed = c(0, 1, 0), VI = seq(0.05, 0.95, 0.05), VS = seq(0.05, 0.5, 0.05), CoIS = 0, relIS = "cor", n.X = NA, autocorr.X = 0, X.dist = "gaussian", intercept = 0, heteroscedasticity = c("null"), mer.sim = TRUE, mer.model = NULL )
EAMM( numsim, group, repl, fixed = c(0, 1, 0), VI = seq(0.05, 0.95, 0.05), VS = seq(0.05, 0.5, 0.05), CoIS = 0, relIS = "cor", n.X = NA, autocorr.X = 0, X.dist = "gaussian", intercept = 0, heteroscedasticity = c("null"), mer.sim = TRUE, mer.model = NULL )
numsim |
number of simulation for each step |
group |
number of group |
repl |
number of replicates per group |
fixed |
vector of lenght 3 with mean, variance and estimate of
fixed effect to simulate. Default: |
VI |
variance component of intercept. Could be specified as a
vector. Default: |
VS |
Variance component of the slope or IxE. Could be specified as a vector.
Default: |
CoIS |
value of correlation or covariance between random intercept and random slope. Default: 0 |
relIS |
"cor" or "cov" set the type of relation give in CoIS. By default the relation is set to correlation |
n.X |
number of different values to simulate for the fixed effect (covariate).
If |
autocorr.X |
correlation between two successive covariate value for a group. Default: |
X.dist |
specify the distribution of the fixed effect. Only "gaussian" (normal distribution) and
"unif" (uniform distribution) are accepted actually. Default: |
intercept |
a numeric value giving the expected intercept value. Default: 0 |
heteroscedasticity |
a vector specifying heterogeneity in residual
variance across X. If |
mer.sim |
Use the simluate.merMod function to simulate the data. Potentially faster for large dataset but more restricted in terms of options |
mer.model |
Simulate the data based on a existing data and model structure from a lmer object. Should be specified as a list of 3 components: a mer object fitted via lmer, an environmental covariate for which to test the random slope, a random effect (e.g. |
P-values for random effects are estimated using a log-likelihood ratio test between two models with and without the effect. Power represent the percentage of simulations providing a significant p-value for a given random structure. Residual variance, e, is calculted as 1-VI.
data frame reporting estimated P-values and power with CI for random intercept and random slope
[PAMM()], [SSF()] for other simulations [plot.EAMM()] for plotting output
## Not run: ours <- EAMM( numsim = 10, group = 10, repl = 4, fixed = c(0, 1, 1), VI = seq(0.1, 0.3, 0.05), VS = seq(0.05, 0.2, 0.05) ) plot(ours, "both") (fm1 <- lmer(Reaction ~ Days + (Days | Subject), sleepstudy)) ours2 <- EAMM( numsim = 10, mer.model = list(model = fm1, env = "Days", random = "Subject"), VI = seq(0.3, 0.5, 0.1), VS = seq(0.05, 0.2, 0.05) ) plot(ours2, "both") ## End(Not run)
## Not run: ours <- EAMM( numsim = 10, group = 10, repl = 4, fixed = c(0, 1, 1), VI = seq(0.1, 0.3, 0.05), VS = seq(0.05, 0.2, 0.05) ) plot(ours, "both") (fm1 <- lmer(Reaction ~ Days + (Days | Subject), sleepstudy)) ours2 <- EAMM( numsim = 10, mer.model = list(model = fm1, env = "Days", random = "Subject"), VI = seq(0.3, 0.5, 0.1), VS = seq(0.05, 0.2, 0.05) ) plot(ours2, "both") ## End(Not run)
Given a specific varaince-covariance structure for random effect, the function simulate different group size and assess p-values and power of random intercept and random slope
PAMM( numsim, group, repl, randompart, fixed = c(0, 1, 0), n.X = NA, autocorr.X = 0, X.dist = "gaussian", intercept = 0, heteroscedasticity = c("null"), ftype = "lmer", mer.sim = FALSE )
PAMM( numsim, group, repl, randompart, fixed = c(0, 1, 0), n.X = NA, autocorr.X = 0, X.dist = "gaussian", intercept = 0, heteroscedasticity = c("null"), ftype = "lmer", mer.sim = FALSE )
numsim |
number of simulation for each step |
group |
number of group. Could be specified as a vector |
repl |
number of replicates per group . Could be specified as a vector |
randompart |
vector of lenght 4 or 5, with 1: variance component
of intercept, VI; 2: variance component of slope, VS; 3: residual
variance, VR; 4: relation between random intercept and random
slope; 5: "cor" or "cov" determine if the relation 4 between I ans S is a correlation or a covariance. Default: |
fixed |
vector with mean, variance and estimate of fixed effect to simulate. Default: |
n.X |
number of different values to simulate for the fixed effect (covariate).
If |
autocorr.X |
correlation between two successive covariate value for a group. Default: |
X.dist |
specify the distribution of the fixed effect. Only "gaussian" (normal distribution) and
"unif" (uniform distribution) are accepted actually. Default: |
intercept |
a numeric value giving the expected intercept value. Default:0 |
heteroscedasticity |
a vector specifying heterogeneity in residual variance
across X. If |
ftype |
character value "lmer", "lme" or "MCMCglmm" specifying the function to use to fit the model. Actually "lmer" only is accepted |
mer.sim |
simulate the data using simulate.merMod from lme4. Faster for large sample size but not as flexible. |
P-values for random effects are estimated using a log-likelihood ratio test between two models with and without the effect. Power represent the percentage of simulations providing a significant p-value for a given random structure
data frame reporting estimated P-values and power with CI for random intercept and random slope
\@seealso [EAMM()], [SSF()], [plot.PAMM()]
## Not run: ours <- PAMM(numsim = 10, group = c(seq(10, 50, 10), 100), repl = c(3, 4, 6), randompart = c(0.4, 0.1, 0.5, 0.1), fixed = c(0, 1, 0.7)) plot(ours,"both") ## End(Not run)
## Not run: ours <- PAMM(numsim = 10, group = c(seq(10, 50, 10), 100), repl = c(3, 4, 6), randompart = c(0.4, 0.1, 0.5, 0.1), fixed = c(0, 1, 0.7)) plot(ours,"both") ## End(Not run)
Provide graphic interpretation of the simulation results
## S3 method for class 'EAMM' plot(x, graphtype = "both", vi, vs, fun3D = "wireframe", ...)
## S3 method for class 'EAMM' plot(x, graphtype = "both", vi, vs, fun3D = "wireframe", ...)
x |
an EAMM object |
graphtype |
"VI", "VS", or "both" "VI" give graphs with varying variance component of intercept and with a fixed variance component for slope specified in vs argument "VS" give graphs with varying variance component for slope and with a fixed variance component of intercept specified in vi argument "both" 3-D plot see also fun3D argument |
vi |
VI for which plots the output. Necessary for "VS" type of graph |
vs |
VS for which plots the output. Necessary for "VI" type of graph |
fun3D |
plot function used to plot the 3D graph. "wireframe" uses lattice, "persp" uses graphics, and "open3d" uses rgl. |
... |
potentially further arguments to pass to methods |
## Not run: ours <- EAMM(numsim=10, group=10, repl=4, VI=seq(0.1,0.95,0.05), VS=c(0.05,0.1) ) plot(ours, "both") plot(ours, "VI",vs=0.1) plot(ours,"VS",vi=0.2) ## End(Not run)
## Not run: ours <- EAMM(numsim=10, group=10, repl=4, VI=seq(0.1,0.95,0.05), VS=c(0.05,0.1) ) plot(ours, "both") plot(ours, "VI",vs=0.1) plot(ours,"VS",vi=0.2) ## End(Not run)
Provide graphic interpretation of the simulation results
## S3 method for class 'PAMM' plot(x, graphtype = "both", nbgroup, nbrepl, fun3D = "wireframe", ...)
## S3 method for class 'PAMM' plot(x, graphtype = "both", nbgroup, nbrepl, fun3D = "wireframe", ...)
x |
a PAMM object |
graphtype |
"group", "repl" or "both" "group" give graphs with varying number of ID and with a fixed number of replicates specified in nbrepl "repl" give graphs with varying number of replicates and with a fixed number of ID specified in nbgroup "both" 3-D plot. see also fun3D argument. Note: useful only with multiple group size and multiple number of replicates. |
nbgroup |
number of group for which plots the output. Necessary for "repl" type of graph |
nbrepl |
number of replicates for which plots the output. Necessary for "group" type of graph |
fun3D |
plot function used to plot the 3D graph. "wireframe" uses lattice, "persp" uses graphics and "open3d" uses rgl |
... |
potentially further arguments to pass to methods |
Parameters phi, theta, ticktype ("simple" or "detailed"), nticks, nbcol from [persp()] function could also be specified for 3D plots. In addition, color schemes ("grey", "cm.colors" and "rainbow") and coltype ("restricted" or "0-1") parameters could also be specified for 3D plots.
## Not run: ours <- PAMM(numsim=10,group=c(seq(10,50,10),100),repl=c(3,4,6), randompart=c(0.4,0.1,0.5,0.1),fixed=c(0,1,0.7)) plot(ours, "both") plot(ours, "group",nbrepl=4) plot(ours,"repl",nbgroup=20) ## End(Not run)
## Not run: ours <- PAMM(numsim=10,group=c(seq(10,50,10),100),repl=c(3,4,6), randompart=c(0.4,0.1,0.5,0.1),fixed=c(0,1,0.7)) plot(ours, "both") plot(ours, "group",nbrepl=4) plot(ours,"repl",nbgroup=20) ## End(Not run)
Provide graphic interpretation of the simulation results
## S3 method for class 'SSF' plot(x, ...)
## S3 method for class 'SSF' plot(x, ...)
x |
an SSF object |
... |
potentially further arguments to pass to methods |
## Not run: oursSSF <- SSF(50,100,10,c(0.4,0.1,0.6,0)) plot(oursSSF) ## End(Not run)
## Not run: oursSSF <- SSF(50,100,10,c(0.4,0.1,0.6,0)) plot(oursSSF) ## End(Not run)
Given a specific total number of observations and variance-covariance structure for random effect, the function simulates different association of number of group and replicates, giving the specified sample size, and assess p-values and power of random intercept and random slope
SSF( numsim, tss, nbstep = 10, randompart, fixed = c(0, 1, 0), n.X = NA, autocorr.X = 0, X.dist = "gaussian", intercept = 0, exgr = NA, exrepl = NA, heteroscedasticity = c("null") )
SSF( numsim, tss, nbstep = 10, randompart, fixed = c(0, 1, 0), n.X = NA, autocorr.X = 0, X.dist = "gaussian", intercept = 0, exgr = NA, exrepl = NA, heteroscedasticity = c("null") )
numsim |
number of simulation for each step |
tss |
total sample size, nb group * nb replicates |
nbstep |
number of group*replicates associations to simulate |
randompart |
vector of lenght 4 or 5 with 1: variance component
of intercept, VI; 2: variance component of slope, VS; 3: residual
variance, VR; 4: relation between random intercept and random
slope; 5: "cor" or "cov" determine id the relation between I ans S is
correlation or covariance, set to |
fixed |
vector of lenght 3 with mean, variance and estimate of fixed effect to simulate |
n.X |
number of different values to simulate for the fixed effect (covariate).
If |
autocorr.X |
correlation between two successive covariate value for a group. Default: |
X.dist |
specify the distribution of the fixed effect. Only "gaussian" (normal distribution) and
"unif" (uniform distribution) are accepted actually. Default: |
intercept |
a numeric value giving the expected intercept value. Default:0 |
exgr |
a vector specifying minimum and maximum value for number of group.
Default: |
exrepl |
a vector specifying minimum and maximum value for number
of replicates. Default: |
heteroscedasticity |
a vector specifying heterogeneity in residual variance
across X. If |
P-values for random effects are estimated using a log-likelihood ratio test between two models with and without the effect. Power represent the percentage of simulations providing a significant p-value for a given random structure
data frame reporting estimated P-values and power with CI for random intercept and random slope
PAMM(), EAMM() for other simulation functions plot.SSF() for plotting
## Not run: oursSSF <- SSF(10, 100, 10, c(0.4, 0.1, 0.6, 0)) plot(oursSSF) ## End(Not run)
## Not run: oursSSF <- SSF(10, 100, 10, c(0.4, 0.1, 0.6, 0)) plot(oursSSF) ## End(Not run)