Linear Mixed Models — pamlmixed • pamlj

power analysis for linear mixed model

pamlmixed(
  aim = "n",
  find = "n",
  syntax = NULL,
  model = NULL,
  focus = NULL,
  clusterpars = list(),
  categorical = list(),
  model_type = "linear",
  sigma2 = 1,
  power = 0.9,
  sig.level = 0.05,
  algo = "mc",
  mcR = 500,
  parallel = TRUE,
  set_seed = FALSE,
  seed = 42,
  run = TRUE,
  verbose = getOption("pamlj.messages"),
  show_data = FALSE,
  ...
)

Arguments

aim: The aim of the analysis: n (default) sample size, power to estimate power
find: When aim='n', indicates whether to find number of clusters find='k' or number of cases within each cluster find='n' (default).
syntax: The model to be analysed with possible options. A random-effect term can attach a one-letter symbol to a coefficient, e.g. `(a*1*1+b*1*x|cluster)`, so that two such symbols can be correlated with a `cor(a,b)=<correlation>` or `cov(a,b)=<covariance>` command line; both symbols must label single-value (non-categorical) random terms of the same cluster.
model: A fitted lme4::lmer() or lme4::glmer() (binomial) model. When supplied, syntax, clusterpars, categorical, sigma2 and model_type are extracted from it (any of these the user also passes explicitly are kept instead), and aim defaults to "power" for the design observed in the fitted model (cluster sizes/levels read off the fit). Random effects, including any covariance between single-value (non-categorical) random terms of the same group, are read off `VarCorr(model)` and reproduced exactly (as `cor()`/`cov()` commands in the extracted syntax, see syntax above); covariance touching a categorical random term cannot be represented and is dropped, with a warning. Categorical predictors are simulated with the exact contrast coding read off the fitted model's own data.frame (see categorical's `coding="custom"` below), so the extracted fixed-effect coefficients remain correctly interpretable. Only Gaussian and binomial mixed models are supported. When the fit has two or more (crossed) clustering/grouping factors (e.g. `subject` and `item`), only the first carries the shared replication factor needed to reproduce the fit's total sample size, so the simulated dataset is not needlessly squared in size; a predictor found to be constant within every level of one of the clustering factors (e.g. a stimulus property fixed per `item`) is automatically marked `between: var|cluster`, so it isn't required to vary within every combination of clusters. A predictor that genuinely varies within every cluster (a true trial-level, counterbalanced factor) has no such shortcut and still requires enough real replication per combination of clusters to realize its levels; if the fit doesn't have that, model construction fails with "N per cluster too small for the planned design".
focus: Only used together with model: the name of the fixed-effect term (as it appears in attr(terms(model), "term.labels"), e.g. "Days", "grp", or an interaction label such as "x:z"; "1"/"(Intercept)" selects the intercept) to focus the search/power computation on, instead of the default of using the worst-case (minimum) power across all effects. Equivalent to attaching a symbol to a term and adding a test: command when writing syntax by hand.
clusterpars: A named list of the form `list(cluster1=c(n=n1,k=k1))`, where `cluster1` is the name of the clustering variable in the model, `n1` is the expcted number of cases within each cluster, and `k1` is the expcted number of clusters. if aim=n, `n1` is used as starting point for sample size. If aim=clusters, `k1` is used as starting point for number of clusters.
categorical: A named list specifying which variable is categorical. Each entry can be either a bare number of levels, `list(varname=x)`, which uses `deviation` (sum-to-zero) coding, or `list(varname=list(levels=x, coding=c))`, which additionally selects a contrast coding scheme `c`: `deviation` (default), `simple`, `dummy`, `difference`, `helmert`, `repeated`, or `polynomial`. `coding="custom"` together with `contrasts=<matrix>` (an x by x-1 numeric contrast matrix, e.g. as returned by `contr.sum()`/`contr.treatment()`/etc.) assigns that matrix verbatim; this is how `model=` extracts a fitted model's own coding, and can also be supplied directly. Any variable in the model not mentioned in `categorical` is assumed to be numeric.
model_type: The model type or family: `linear` (default) for linear mixed model, `logistic` for binomial logistic mixed model.
sigma2: Residual variance. Ignored for `model_type="logistic"`
power: Minimal desired power
sig.level: Type I error rate (significance cut-off or alpha)
algo: The algorithm to use: `mc` (default) for Monte Carlo simulation, `raw` for raw approximation based on Chi-squared (fast but not very accurate)
mcR: Number of repetitions for Monte Carlo method
parallel: Logical: should parallel computing be used for the Monte Carlo method
seed: the seed for Monte Carlo simulations, default=42.
run: TRUE (default) run the simulations, otherwise print out the model without results
verbose: (Boolean) `getOption("pamlj.messages")` (default). Print out updates of the simulation steps.
...: Used for internal purposes
Number: of repetitions for Monte Carlo method

Value

A results object containing:

`results$intro`					a html
`results$extrainfo`					a html
`results$issues`					a html
`results$initnotes`					a html
`results$infotab`					a table
`results$powertab`					a table
`results$plotnotes`					a html

Tables can be converted to data frames with asDF or as.data.frame. For example:

results$infotab$asDF

as.data.frame(results$infotab)