afex (version 0.6-82)

mixed: Obtain p-values for a mixed-model from lmer().

Description

Fits and calculates p-values for all effects in a mixed model fitted with lmer. The default behavior calculates type 3 like p-values using the Kenward-Rogers approximation for degrees-of-freedom implemented in KRmodcomp (for LMMs only), but also allows for parametric bootstrap (method = "PB") (for LMMs and GLMMs). print, summary, and anova methods for the returned object of class "mixed" are available (all return the same data.frame).

Usage

mixed(formula, data, type = 3, method = c("KR", "PB",
    "LRT"), per.parameter = NULL, args.test = list(),
    check.contrasts = TRUE, progress = TRUE, cl = NULL,
    ...)

Arguments

formula

a formula describing the full mixed-model to be fitted. As this formula is passed to lmer, it needs at least one random term.

data

data.frame containing the data. Should have all the variables present in fixed, random, and dv as columns.

type

type of test on which effects are based. Only type 3 tests (3 or "III") are correctly implemented (see Details).

method

character vector indicating which methods for obtaining p-values should be used. "KR" (the default) corresponds to the Kenward-Rogers approximation for degrees of freedom (only working with linear mixed models). "PB"

per.parameter

character vector specifying for which variable tests should be run for each parameter (instead for the overall effect). Can be useful e.g., for testing ordered factors. Relatively untested so results should be compared with a seco

args.test

list of arguments passed to the function calculating the p-values. See details.

check.contrasts

logical. Should contrasts be checked and (if necessary) changed to be "contr.sum". See details.

progress

if TRUE, shows progress with a text progress bar

cl

A vector identifying a cluster; used for distributing the estimation of the different models using several cores. See examples.

...

further arguments (such as weights) passed to lmer.

Value

• An object of class "mixed" (i.e., a list) with the following elements:
  1. anova.tablea data.frame containing the statistics returned fromKRmodcomp.
  2. full.modelthe"lmerMod"object returned from fitting the full mixed model.
  3. restricted.modelsa list of"lmerMod"objects from fitting the restricted models (i.e., each model lacks the corresponding effect)
  4. testsa list of objects returned by the function for obtaining the p-values.
  5. typeThetypeargument used when calling this function.
  6. methodThemethodargument used when calling this function.
  The following methods exist for objects of class "mixed": print (which uses rounding and only returns the results with F.scaling = 1), summary, and anova (the latter two return the same data.frame).

Details

Type 3 tests are obtained by comparing a model in which only the corresponding effect is missing with the full model (containing all effects). This corresponds to the (type 3) Wald tests given by car::Anova for "lmerMod" models. Type 2 tests are obtained by comparing a model in which the corresponding effect and all higher oder effect (e.g., all three-way interactions for a two-way interaction) are missing with a model in which all effects of the relevant order are present and all higher order effects absent. Consequently, the results for lower order effects are identical of wether or not higher order effects are part of the model or not, which is rather dubious (but https://stat.ethz.ch/pipermail/r-sig-mixed-models/2012q3/018992.html{I didn't find a better way} of implementing the Type 2 tests). This does not correspond to the (type 2) Wald Test reported by car::Anova. If you want type 2 Wald tests, use car::Anova with test = "F" (from version 2.0-13) instead of this function. For an introduction to mixed-modeling for experimental designs using p-values see Judd, Westfall, and Kenny (2012). Further introductions to mixed-modeling for experimental designs are given by Baayen and colleagues (Baayen, 2008; Baayen, Davidson & Bates, 2008; Baayen & Milin, 2010). Recommendations on how to specify the random effects structure for experimental designs can be found in Barr and colleagues (2013). p-values are per default calculated via methods from pbkrtest. When method = "KR", the Kenward-Rogers approximation for degrees-of-freedom is calculated using KRmodcomp, which is only applicable to linear-mixed models. The possible argument to pbkrtest via args.test is details. method = "PB" calculates p-values using parametric bootstrap using PBmodcomp. This can be used for linear and also generalized linear mixed models (GLMM) by specifiying a family argument to mixed. Note that you should specify further arguments to PBmodcomp via args.test, especially nsim (the number of simulations to form the reference distribution. For other arguments see PBmodcomp. Note that REML (argument to lmer) will be set to FALSE if method is PB. method = "LRT" calculates p-values via likelihood ratio tests implemented in the anova method for "merMod" objects. This is recommended by Barr et al. (2013) which did not test the other methods implemented here. Furthermore, this is not recommended by the http://glmm.wikidot.com/faq{lme4 faq} which instead recommends the other methods implemented here. If check.contrasts = TRUE, contrasts will be set to "contr.sum" for all factors in the formula if default contrasts are not equal to "contr.sum" or attrib(factor, "contrasts") != "contr.sum".

References

Baayen, R. H. (2008). Analyzing linguistic data: a practical introduction to statistics using R. Cambridge, UK; New York: Cambridge University Press. Baayen, R. H., Davidson, D. J., & Bates, D. M. (2008). Mixed-effects modeling with crossed random effects for subjects and items. Journal of Memory and Language, 59(4), 390-412. doi:10.1016/j.jml.2007.12.005 Baayen, R. H., & Milin, P. (2010). Analyzing Reaction Times. International Journal of Psychological Research, 3(2), 12-28. Barr, D. J., Levy, R., Scheepers, C., & Tily, H. J. (2013). Random effects structure for confirmatory hypothesis testing: Keep it maximal. Journal of Memory and Language, 68(3), 255–278. doi:10.1016/j.jml.2012.11.001 Judd, C. M., Westfall, J., & Kenny, D. A. (2012). Treating stimuli as a random factor in social psychology: A new and comprehensive solution to a pervasive but largely ignored problem. Journal of Personality and Social Psychology, 103(1), 54-69. doi:10.1037/a0028347

See Also

ez.glm and aov.car for convenience functions to analyze experimental deisgns with classical ANOVA or ANCOVA wrapping Anova.

Examples

# use the obk.long data (mildly reasonable)
data(obk.long)
mixed(value ~ treatment * phase + (hour|id), obk.long)

# Examples for using the per.parammeter argument:
data(obk.long, package = "afex")
obk.long$hour <- ordered(obk.long$hour)

# tests only the main effect parameters of hour individually per parameter.
mixed(value ~ treatment*phase*hour +(1|id), per.parameter = "^hour$", data = obk.long)

# tests all parameters including hour individually
mixed(value ~ treatment*phase*hour +(1|id), per.parameter = "hour", data = obk.long)

# tests all parameters individually
mixed(value ~ treatment*phase*hour +(1|id), per.parameter = ".", data = obk.long)

# example data from package languageR:
# Lexical decision latencies elicited from 21 subjects for 79 English concrete nouns,
# with variables linked to subject or word.
data(lexdec, package = "languageR")

# using the simplest model
m1 <- mixed(RT ~ Correct + Trial + PrevType * meanWeight +
    Frequency + NativeLanguage * Length + (1|Subject) + (1|Word), data = lexdec)

m1
# gives:
##                   Effect df1       df2      Fstat p.value
## 1            (Intercept)   1   96.6379 13573.1410  0.0000
## 2                Correct   1 1627.7303     8.1452  0.0044
## 3                  Trial   1 1592.4301     7.5738  0.0060
## 4               PrevType   1 1605.3939     0.1700  0.6802
## 5             meanWeight   1   75.3919    14.8545  0.0002
## 6              Frequency   1   76.0821    56.5348  0.0000
## 7         NativeLanguage   1   27.1213     0.6953  0.4117
## 8                 Length   1   75.8259     8.6959  0.0042
## 9    PrevType:meanWeight   1 1601.1850     6.1823  0.0130
## 10 NativeLanguage:Length   1 1555.4858    14.2445  0.0002

# Fitting a GLMM using parametric bootstrap:
require("mlmRev") # for the data, see ?Contraception

gm1 <- mixed(use ~ age + I(age^2) + urban + livch + (1 | district),
 family = binomial, data = Contraception, args.test = list(nsim = 10), method = "PB")

## using multicore
require(parallel)
(nc <- detectCores()) # number of cores
cl <- makeCluster(rep("localhost", nc)) # make cluster
# to keep track of what the function is doind, redirect output to outfile:
# cl <- makeCluster(rep("localhost", nc), outfile = "cl.log.txt")

# obtain fits with multicore:
mixed(value ~ treatment*phase*hour +(1|id), data = obk.long, method = "LRT", cl = cl)