Getting Started

This guide explains what each function in EDF9850_Functions.R does, what its options mean, and when to use it. Each entry explains not just how to call the function, but what the statistics it prints actually mean, since these functions assume no prior R experience.

Loading the functions

Before using any function below, run this once at the start of your R session (adjust the path to wherever you saved the file, or load it directly from the web):

source("EDF9850_Functions.R")
source("https://www.lesahoffman.com/EDF9850/EDF9850_Functions.R") # From the web

This also loads the R packages these functions depend on—nlme, lme4, lmerTest, performance, emmeans, and ordinal—installing them automatically the first time if they are not already on your computer, without the usual "masked" package-loading notices some of them would otherwise print (they are harmless naming overlaps between these packages, not something to act on).

Helpful error messages

If you call one of these functions on an object that does not exist yet (for example, a typo in a model's name, or a model you have not fit yet), or on the wrong kind of model (for example, handing an lmer model to glsSummary()), you will get a plain-English error message explaining the problem, instead of a confusing low-level R error. Just fix the problem the message describes and run the line again—nothing else is required.

Significance stars (optional)

By default these functions print p-values only. To add significance stars to every p-value column in every function at once, run this once in your R session:

options(EDF9850.signif.stars = TRUE)    # turn stars on for all of these functions
options(EDF9850.signif.stars = FALSE)   # turn them back off (the default)

Stars use R's usual cutoffs: *** for p < .001, ** for p < .01, * for p < .05, and a period (.) for p < .1. This setting is separate from R's own show.signif.stars option, so it is unaffected by how that is set. In each table with p-values, the stars go in their own unlabeled column at the end of each row, after all the other columns (in corrSummary()'s correlation matrix, which has no single row to end, each star follows its own p-value instead). The likelihood-ratio (LRT) tables never get stars, in any function, and the saved/returned objects are never changed—stars affect only what is printed.

A Few Terms Used Throughout This Guide

These functions report several quantities that are specific to multilevel/mixed-effects models. Brief plain-language definitions are collected here so each function's description below can stay short.

• Fixed effects — Model terms (the intercept and predictors) whose single estimated value is assumed to apply to every person/observation in the data.
• Random effects — The extent to which the intercept (and, if modeled, one or more slopes) varies from cluster to cluster (e.g., from school to school).
• ICC1 — Intraclass correlation: the proportion of the total outcome variance that lies between clusters rather than within them.
• ICC2 — The reliability of the cluster means—how trustworthy a given cluster's average is as an estimate of that cluster's true mean.
• LRT — Likelihood ratio test: a significance test for whether a more complex model fits significantly better than a simpler, nested model.
• Mixture p-value (p.mix) — A corrected p-value for an LRT that removes a variance component. Because a variance cannot go below 0, its true null distribution is a 50:50 mixture of chi-square(df) and chi-square(df−1), not plain chi-square(df) (Self & Liang, 1987)—for a 1-df variance test this works out to exactly half the ordinary p-value. It appears alongside the regular p-value in every LRT table these functions produce.
• Design effect / Effective N — The design effect is how much a cluster sample's precision is degraded by clustering, computed as 1 + (mean cluster size − 1) × ICC1. The effective N is the total sample size divided by the design effect—the "equivalent" sample size if the data were not clustered at all.
• Slope reliability (SR) — The reliability of a random slope: how trustworthy a given cluster's own slope is as an estimate of that cluster's true slope. Computed as slope variance / (slope variance + residual variance / (mean L1n × variance of the level-1 predictor)), where L1n is the number of level-1 observations per cluster; the predictor should be centered within cluster. Caveat: the residual-variance term approximates each cluster's own slope sampling variance as if the predictor's within-cluster sum of squares were a fixed, known quantity; in practice it varies from cluster to cluster, which understates that sampling variance somewhat (more so in small clusters) and inflates SR a little as a result. It also assumes the predictor is uncorrelated with any other predictor sharing the same random-effects block, and does not apply at all, even approximately, to a random slope on an interaction/product term (e.g. a random X1:X2 slope)—treat SR for those with real caution.
• Satterthwaite df — An approximate denominator degrees-of-freedom used for each fixed effect's significance test. Unlike ordinary regression, multilevel models do not have one simple, exact df, so this approximation is used instead.
• Cohen's d / partial r — Two standardized effect-size measures for a fixed effect, letting you judge its size on a common scale regardless of the outcome's original units.
• Odds ratio / rate ratio — The GLMM (glmerSummary()) equivalent of Cohen's d for a Wald z-test: exp(Estimate). For a logit-link (e.g. logistic) model this is an odds ratio; for a log-link (e.g. Poisson) model, a rate ratio. A value of 1 means no effect; above 1 means the outcome becomes more likely/frequent as the predictor increases, below 1 means less.
• CI — Confidence interval: a range of plausible values for a parameter, at a chosen confidence level (95% by default in these functions).
• Pseudo-R2 — The percent reduction in one variance component (e.g., the random intercept variance) when predictors are added, comparable in spirit to R2 in ordinary regression but computed separately for each variance component.
• Total R2 — The percent of the total outcome variance explained by a model's fixed effects only.

univSummary()

Reports the seven most basic descriptive statistics for one or more numeric variables in a data frame: N (the non-missing sample size), Mean, SD (standard deviation), Variance, Min, Max, and Median. This is a simpler, more focused alternative to functions like psych::describe() when this basic summary is all you need. Optionally split the results by one or more grouping variables (similar to psych::describeBy()), which produces one row per variable-by-group combination in a single table rather than a separate table for each group. The Median is the middle observation when N is odd; when N is even, medianType (see below) controls how the two middle observations are combined.

Basic usage

univSummary(MyData)
univSummary(MyData, vars=c("mpg","hp"), by="cyl")

Arguments

corrSummary()

Reports a Pearson correlation matrix for two or more numeric variables in a data frame, together with the p-value testing whether each correlation differs from zero. By default (p.separate=FALSE), the table is stacked: each variable's row of correlations is followed immediately by an indented row giving that row's p-values in parentheses, so every p-value sits directly beneath its correlation. Set p.separate=TRUE to instead print the correlations and p-values as two separate plain tables, with no stacking. Each pair of variables is correlated using only the rows where both variables are observed, so the sample size can differ from pair to pair when there is missing data.

Basic usage

corrSummary(MyData, vars=c("mpg","hp","wt"))
corrSummary(MyData, vars=c("mpg","hp","wt"), p.separate=TRUE)

Arguments

freqSummary()

Reports the frequency (N) and proportion of each observed level for one or more categorical variables in a data frame—analogous to univSummary(), but for categorical rather than numeric variables. Produces one row per variable-by-level combination in a single table (one row per variable-by-level-by-group combination when by is used), rather than a separate table for each variable. table()'s own level order is kept—a factor's own level order (including any level with zero occurrences, shown as N=0), or sorted order for a character/logical/numeric variable—and Proportion always sums to 1 within each variable (and group, when by is used).

Basic usage

freqSummary(MyData)
freqSummary(MyData, vars="lunch", by="schoolID")

Arguments

addUnitMeans()

Use this before fitting a model when you need cluster-level information added to your raw (level-1) data—most commonly, the number of level-1 units in each cluster and the cluster mean of one or more level-1 variables (a common way to create a level-2 "contextual" predictor, such as a school's mean SES). It returns a new data frame with these columns added to every row.

Basic usage

Example3 = addUnitMeans(data=Example3, unitVariable="schoolID",
                        meanVariables=c("ses"), newNames=c("CMses"))

Arguments

What it adds: a column named "Nper<unitVariable>" giving the number of level-1 rows in that row's cluster, plus one new column per entry in meanVariables/newNames giving that cluster's mean on that variable, merged onto every row belonging to that cluster.

lmerSummary()

Prints a complete, course-formatted summary of a model fit with lmer(): the model fit statistics (AIC, BIC, log-likelihood, -2*logLik), the random effects table, an omnibus F-test that all fixed slopes are jointly zero, and the fixed effects table. Several additional pieces of output—intraclass correlations, a design effect and effective sample size, random slope reliability, random-effect confidence intervals, standardized effect sizes, and a likelihood-ratio test for the random effects—can be turned on as needed. It also automatically creates new objects in your R environment holding each variance/covariance component from the model, so you do not need to look them up by hand for further calculations. If a random-effect variance is estimated at exactly 0 (a singular fit), every correlation involving it is undefined and prints as NaN, with a note under the random effects table; glmerSummary() and clmmSummary() do the same.

Basic usage

lmerSummary(Model1, name="Model 1: Empty Means Model", ICC1=TRUE, LRT=TRUE)

Arguments

What it saves automatically

When saveObjects=TRUE, using the model's own name (e.g. a model called Model1): Model1_VarCorr (the full variance-component table) and Model1_Fixed (the full fixed-effects table), plus one object per variance/covariance component, e.g. Model1_tau0 (random intercept variance), Model1_sigma2 (residual variance), and Model1_tau_<predictor> for any random slope variance. If the model has more than one grouping factor (e.g. students in classrooms in schools), each of those names also gets its grouping factor appended so the factors do not overwrite each other, e.g. Model1_tau0_school and Model1_tau0_class (a factor such as class:school appears as class.school).

glsSummary()

Same idea as lmerSummary(), but for a gls() model—a single-level regression model with no random effects at all. This is most often the comparison/baseline model without any clustering. Prints the model fit statistics, the residual variance, an omnibus F-test that all fixed slopes are jointly zero, and the fixed effects table. For a model fit with weights=varIdent(form=~1|group) (separate residual variances by group, e.g. by treatment arm or zygosity), it also automatically prints a Residual Variance by Group table—the single residual variance above is then only the reference group's. In a no-intercept model with a factor entered as cell means (e.g. y ~ 0 + group + x), that factor's coefficients are the group means, not slopes, so they are left out of the omnibus test that all slopes are zero. By default each fixed effect is given its own Satterthwaite df and the omnibus test a Satterthwaite denominator df (as lmeSummary() does), which matters for a gls model with an estimated residual covariance matrix (e.g., an unstructured R matrix from corSymm() with varIdent(), or an AR1 structure); use ddf="residual" for the older approach of N minus the number of fixed effects for every effect. The Fixed Effects table is followed by a line stating which denominator df method was used. For a model with a residual correlation structure (for example corCompSymm, corAR1, or corSymm), it also prints a Residual Correlation table—the structure, each parameter's name, and its estimate—right after the residual variance; for compound symmetry, Rho is the correlation among all of a person's residuals (the ICC of the equivalent random-intercept model).

Basic usage

glsSummary(Model0, name="Model 0: No Clustering", CI=TRUE)

# Residual df (N - p) for every effect, instead of the default Satterthwaite df
glsSummary(Model2b, ddf="residual")

Arguments

What it saves automatically

When saveObjects=TRUE: Model0_Fixed (the fixed-effects table) and Model0_sigma2 (the residual variance). For a varIdent() model, also Model0_ResidByGroup (the by-group table) and Model0_sigma2_<group> for each group's own residual variance.

lmeSummary()

Serves the same purpose as lmerSummary(), for a multilevel model fit with nlme's lme() function instead of lme4's lmer(). Supports the same ddf/ICC1/Neff/ICC2/slopeRel/randCI/randCIfixed/randCIrandom/LRT/effectSizes/CI options as lmerSummary(), except LRTdrop/LRTdropGroup—nlme's random-effects structure has no ranova()-style per-term breakdown for a custom row to fit into, so LRT here is a set of separate comparisons, described below: just one (the random effects) for a model with only a random intercept, and otherwise one for each piece (each random slope, the set of random slopes, and the residual correlation). One difference under the hood: by default it replaces nlme's own degrees-of-freedom with a genuine Satterthwaite approximation (computed via the emmeans package) so the significance tests match what lmerSummary() reports for an equivalent lmer model. This is done for every fixed effect separatelythe intercept, continuous slopes, factor- or contrast-coded terms, and interactions alikeso each row of the Fixed Effects table has its own df. The overall test that all fixed slopes are zero uses the same Satterthwaite approximation, extended to several df at once, so it is an F-test even when the slopes are at different levels of the model (e.g. a twin-level and a pair-level predictor, where nlme's own anova(Terms=) cannot combine them); with ddf="asymptotic", or if the Satterthwaite df cannot be computed, it is a multivariate Wald chi-square test instead (labeled accordingly). In a no-intercept model with a factor entered as cell means (e.g. y ~ 0 + group + x), that factor's coefficients are the group means, not slopes, so they are left out of this test. For a model fit with weights=varIdent(form=~1|group) (separate residual variances by group, e.g. twins by zygosity), it also automatically prints a Residual Variance by Group table, and the ICC1/Neff/ICC2/slopeRel/randCI sections below automatically switch to a by-group breakdown instead of a single pooled value. For a model with a residual correlation structure (for example corAR1, corCompSymm, or corSymm), its parameter(s) are added to the Model for the Variance table as extra Residual rows in the Corr column (e.g., AR1 Phi), because nlme keeps them outside the random-effects variances; they are for display only and are not saved by saveObjects.

Basic usage

lmeSummary(Model2c, ICC2=TRUE, randCI=TRUE)

Arguments

What it saves automatically

When saveObjects=TRUE, the same naming scheme as lmerSummary(): Model2c_VarCorr, Model2c_Fixed, Model2c_tau0, Model2c_sigma2, and so on for each variance component. For a varIdent() model, also Model2c_ResidByGroup (the by-group residual variance table) and Model2c_sigma2_<group> for each group's own residual variance. With nested grouping factors, each variance/covariance name also gets its factor appended (e.g. Model2c_tau0_school and Model2c_tau0_class), and the saved variance table gains a grp column.

glmerSummary()

Parallel in structure to lmerSummary(), but for a model fit with glmer()—lme4's version for a binary, count, or other non-Gaussian outcome (e.g. a logistic multilevel model for a yes/no outcome). Several things that make sense for lmerSummary()'s Gaussian models do not translate directly to a GLMM, so glmerSummary() adapts them:

• No Satterthwaite/Kenward-Roger df — glmer() fixed effects are tested with Wald z-tests, so there is no ddf argument here; p-values and CIs use the normal distribution.
• Omnibus test — the joint test that all fixed slopes are zero is a multivariate Wald chi-square test (labeled accordingly) rather than an F-test, since glmer models have no residual-variance-based denominator df.
• Effect sizes — Cohen's d/partial r do not apply to a Wald z-test. Instead, for a logit-link model this prints the odds ratio (OR = exp(Estimate)) with its own delta-method SE (OR.SE = OR×SE) and, when CI=TRUE, a 95% CI—under the header "Effect Sizes for Fixed Slopes in Odds Ratios (OR)". For a log-link model (e.g. Poisson), the rate ratio (RR/RR.SE) instead, under the parallel "...in Rate Ratios (RR)" header. For any other link (e.g. probit), a note prints instead of a table, since an exponentiated effect size is not meaningful there. The CI comes from exponentiating the logit/log-scale CI's own endpoints, not from Estimate±1.96×SE, since that symmetric interval is not valid on a ratio scale.
• Predicted probability (or mean) — the Fixed Effects table also gets a Prob/Prob.SE column pair (plus Prob.Lower/Prob.Upper when CI=TRUE), converting the intercept's own linear-predictor-scale estimate to its actual predicted value via the model's own inverse link function, with a delta-method SE—blank for every slope row, since a slope has no single value of its own to convert. For a binomial model this is a genuine probability (bounded 0–1); for any other family (e.g. Poisson), family$linkinv() instead gives a predicted mean count/rate, so the columns are labeled Pred.Count/Pred.Count.SE there instead.
• ICC2 / Neff (design effect) — both require a level-1 ("residual") variance, which a GLMM does not estimate the way a Gaussian model does. These use the standard fixed constant instead—pi^2/3 for a logit link, 1 for a probit link (the same convention performance::icc() uses internally for ICC1), pi^2/6 for a cloglog link—shown as its own row in the Model for the Variance table, labeled with which constant applies and with its Variance/Std.Dev. in parentheses to flag it as assumed rather than estimated. Unavailable, with a clear note, for any other family/link.
• LRT — delegates directly to glmerRanova() (see below), which already does exactly this for glmer models.

Basic usage

EmptyNon = glmer(data=Example3, family=binomial(link="logit"), nAGQ=7,
                 formula=nonwhite~1+(1|schoolID))
glmerSummary(EmptyNon, ICC1=TRUE, effectSizes=TRUE, CI=TRUE)

Arguments

What it saves automatically

When saveObjects=TRUE, the same naming scheme as lmerSummary(): <model>_VarCorr, <model>_Fixed, <model>_tau0, and so on for each variance component, plus <model>_sigma2 for the assumed level-1 constant when the link supports one (logit, probit, or cloglog). With more than one grouping factor, each variance/covariance name also gets its factor appended (e.g. <model>_tau0_school, <model>_tau0_class).

clmmSummary()

Parallel in structure to glmerSummary(), but for a model fit with clmm()—the ordinal package's cumulative link mixed model, the multilevel extension of an ordinal/proportional-odds logistic regression, for an outcome with 3 or more ordered categories (e.g. a Likert-type item). Like glmer models, clmm models are tested with Wald z-tests (no Satterthwaite/Kenward-Roger df) and have no true residual variance, so most of glmerSummary()'s same adaptations apply here too (see that entry above): a Wald chi-square omnibus test, an odds ratio (logit link only) for effect sizes, and an assumed level-1 variance constant standing in for ICC1/ICC2/Neff/slope reliability. The one real structural difference from glmer: an ordinal outcome with k categories has k−1 threshold (cutpoint) parameters instead of a single intercept. clmm() fits proportional-odds models only (it has no nominal= or scale= options and no test of proportional odds); for a model without random effects, or to relax proportional odds, use clmSummary() below.

• Intercepts / Thresholds — printed as two small tables, since the k−1 cutpoints replace the role an intercept plays elsewhere (both are excluded from the omnibus test, the Fixed Effects table, effect sizes, and reference-level detection, the same way an intercept already is everywhere else). Thresholds are clmm's own parameterization, predicting the chosen link of the probability of being in the LOWER category (at or below each cutpoint). Intercepts are the same cutpoints sign-flipped, predicting the probability of being in the HIGHER category—pointing the same direction as the slopes, which is why Intercepts prints first, immediately before Fixed Effects. Thresholds prints last, after Effect Sizes, and only when thresholds=TRUE (off by default, since Intercepts already carries the same information in the more intuitive direction). Both tables label each row with every category on each side of that cutpoint, e.g. a 3-category outcome's two cutpoints are labeled "0|1,2" and "0,1|2", and both add a Prob/Prob.SE column pair (plus Prob.Lower/Prob.Upper when CI=TRUE) converting that row's estimate to an actual probability via the model's own inverse link function, with a delta-method SE.
• Effect sizes — the odds ratio (OR = exp(Estimate)) with its own delta-method SE (OR.SE), and a 95% CI when CI=TRUE, for each fixed slope—logit link only, since an exponentiated effect size is not meaningful for clmm's other four links, where a note prints instead. Header reads "Effect Sizes for Fixed Slopes in Odds Ratios (OR)" when shown.
• Link scope — the assumed level-1 variance constant covers all five links clmm() supports: pi^2/3 for logit, 1 for probit, pi^2/6 for cloglog or loglog, and unavailable for cauchit (the Cauchy distribution has no finite variance).
• LRT — unlike glmerSummary() (which delegates to glmerRanova()), this is a ranova()-style table built directly for clmm models: one row per random-effects term in the model, each refit with that whole term dropped (falling back to a plain clm() model when no random effects remain) and compared to the full model; a term that has random slopes gets one row per slope instead (that slope and its covariances dropped, the rest of the term kept) and no separate intercept row, as in lmerSummary()'s ranova() table, with the same p.mix mixture-corrected p-value used everywhere else in this guide. Each row's "Random Effect Dropped" label is condensed to "Term | Group" form, e.g. "(1 | schoolID)" is shown as "Intercept | schoolID". LRTdrop/LRTdropGroup add one extra custom row, the same as lmerSummary()'s (see that entry above).

Basic usage

EmptyLunch = clmm(data=Example4b, link="logit", formula=as.factor(lunch)~1+(1|schoolID))
clmmSummary(EmptyLunch, ICC1=TRUE, effectSizes=TRUE, CI=TRUE, thresholds=TRUE)

Arguments

What it saves automatically

When saveObjects=TRUE, the same naming scheme as glmerSummary(): <model>_VarCorr, <model>_Fixed, <model>_Thresholds, <model>_Intercepts, <model>_tau0, and so on for each variance component, plus <model>_sigma2 for the assumed level-1 constant when the link supports one (logit, probit, cloglog, or loglog). With more than one grouping factor, each variance/covariance name also gets its factor appended (e.g. <model>_tau0_school, <model>_tau0_class).

clmSummary()

clmmSummary()'s counterpart for an ordinal outcome modeled without random effects, with clm()—the ordinal package's cumulative link (proportional-odds) model, for an outcome with 3 or more ordered categories. It prints the same cutpoint and slope tables as clmmSummary(): an Intercepts table (the k-1 cutpoints sign-flipped, giving the link of the probability of being in the HIGHER category after each cutpoint, so they point the same direction as the slopes), the Fixed Effects table with Wald z-tests, an overall Wald chi-square test that all slopes are zero, and optionally odds ratios (logit link only), model-predicted probabilities, and a Thresholds table. It has no variance, ICC, or random-effect output, since there are no random effects.

It also reports what only clm() can do. A predictor given in clm()'s nominal= argument has its own effect at every cutpoint instead of one slope (non-proportional odds); these get a Nominal Effects table, one row per predictor and cutpoint, signed like the slopes (positive = more likely above that cutpoint). A predictor given in scale= changes the spread of the latent variable; these get a Scale Effects table on the log scale, with SD.Ratio = exp(Est), the ratio of the latent SD per unit of the predictor. The options LRT, R2, and parallelTest add likelihood ratio tests, pseudo-R2, and the test of proportional odds (see below). Only clm()'s default threshold="flexible" and the five links logit, probit, cloglog, loglog, and cauchit are supported. (clmm() has no nominal= or scale= options and no test of proportional odds, so these are available only here.)

Basic usage

From the wine data in the ordinal package (ratings of 1 to 5 by temperature and contact). Model2 gives contact its own effect at each cutpoint.

Model1 = clm(rating ~ temp + contact, data=wine)
clmSummary(Model1)
clmSummary(Model1, CI=TRUE, effectSizes=TRUE, predMeans="temp",
           LRT=TRUE, R2=TRUE, parallelTest=TRUE)
Model2 = clm(rating ~ temp, nominal=~contact, data=wine)
clmSummary(Model2)

Output of the first call:

Model Fit for Predicting rating:
      AIC       BIC    logLik -2*logLik
  184.984   198.644   -86.492   172.984
Number of obs: 72

Overall Wald Chi-Square Test for All Fixed Slopes:
 Model  NumDF  Chisq p.value
 Model1     2 27.912   <.001

Intercepts predicting the logit of the probability
of being in the HIGHER category after | :
 Term          Est    SE z.value p.value  Prob Prob.SE
 1|2,3,4,5   1.344 0.517   2.600   0.009 0.793   0.085
 1,2|3,4,5  -1.251 0.438  -2.857   0.004 0.223   0.076
 1,2,3|4,5  -3.467 0.598  -5.800   <.001 0.030   0.018
 1,2,3,4|5  -5.006 0.731  -6.850   <.001 0.007   0.005

Fixed Effects:
 Term                  Est    SE z.value p.value
 temp (ref = cold)
 tempwarm            2.503 0.529   4.735   <.001
 contact (ref = no)
 contactyes          1.528 0.477   3.205   0.001

Some of the extra tables: the next three come from LRT=TRUE, R2=TRUE, and parallelTest=TRUE for Model1, and the last is the Nominal Effects table printed for Model2:

LRT for each fixed slope (full model vs. that term dropped):
 Term Dropped  Npar   LL -2LL AIC BIC    LRT DF p.value
 <none>           6  -86  173 185 199
 temp             5 -100  200 210 221 26.928  1   <.001
 contact          5  -92  184 194 205 11.043  1   0.001
 All slopes       4 -104  207 215 225 34.454  2   <.001

Pseudo-R2 (relative to a model with only the cutpoints):
 McFadden Cox.Snell Nagelkerke
    0.166     0.380      0.403

Test of Proportional Odds (a separate effect at each cutpoint):
 Term    DF   LRT p.value
 temp     3 3.175   0.365
 contact  3 0.567   0.904
  (significant = that predictor's effect differs across cutpoints)

Nominal Effects (a separate effect at each cutpoint):
 Term       Cutpoint     Est    SE z.value p.value
 contactyes 1|2,3,4,5  1.615 1.162   1.390   0.164
 contactyes 1,2|3,4,5  1.512 0.591   2.559   0.010
 contactyes 1,2,3|4,5  1.675 0.649   2.581   0.010
 contactyes 1,2,3,4|5  1.051 0.897   1.172   0.241

Arguments

The results are also returned, invisibly, as a list (Fixed, Thresholds, Intercepts, Nominal, Scale, omnibusWald, effectSizes, fixedCI, predMeans, LRT, R2, and parallelTest, for whichever were requested or apply), e.g. x = clmSummary(Model1).

con1DSummary()

Formats the result of a single linear-combination ("1 df contrast") test—for example, testing whether one group's slope equals another's, or whether a specific combination of fixed effects equals zero—computed first with lmerTest's contest1D() function. It prints the result in the same table layout as the Fixed Effects table in lmerSummary(), so it looks consistent with the rest of your output.

Basic usage

The examples for con1DSummary(), con1DTable(), and conMDSummary() all use the same model, Model4c from Example 4a: math ~ 1 + WCses + WCnon + CMses + CMnon30 + (1|schoolID). Each level-1 predictor enters as a cluster-mean-centered within-cluster slope (WCses, WCnon) and a cluster-mean between-cluster slope (CMses, CMnon30), so its fixed effects, in order, are Intercept, WCses, WCnon, CMses, and CMnon30. Its contextual effects (between minus within) are the contrasts tested below.

Test1 = contest1D(Model4c, L=c(0, -1, 0, 1, 0))   # CMses - WCses
con1DSummary(Test1, name="Contextual SES", CI=TRUE, effectSizes=TRUE)

Fixed Effect: Contextual SES (with 95% CI)
   Est    SE      DF t.value p.value Lower Upper
 3.375 0.417 172.761   8.090   <.001 2.551 4.198

Effect Size:
   Est    SE p.value Cohen.d Partial.r
 3.375 0.417   <.001   1.231     0.524

Arguments

Example: contextual effects from a within-between model (Example 4a)

In a model that has each predictor as a within-cluster (cluster-mean-centered) slope and a cluster-mean slope, the cluster-mean slope is the between-cluster effect, and the contextual effect is the difference between the two: between minus within. Each number in L goes with one fixed effect, in the order of the Fixed Effects table (Example 4a's Model 4c has Intercept, WCses, WCnon, CMses, CMnon30), so the contextual effect of SES, CMses − WCses, is c(0, -1, 0, 1, 0). With the fixed effects of Model 4c (WCses = 1.952 and CMses = 5.327) that is 5.327 − 1.952 = 3.375, along with its SE and test:

# Model 4c: contextual effect of SES = CMses - WCses
Test1 = contest1D(Model4c, L=c(0, -1, 0, 1, 0))
con1DSummary(Test1, name="Contextual SES")

Fixed Effect: Contextual SES
   Est    SE      DF t.value p.value
 3.375 0.417 172.761   8.090   <.001

contest1D() accepts L only by position, as above; con1DTable() (next) accepts the same contrast written by name, and can request several at once.

con1DTable()

Runs several single-df contrasts (tests of one linear combination of the fixed effects each) against the same model, an lmer, lme, or gls model, and prints them as a single table—one row per contrast—instead of one contest1D()/con1DSummary() call and printed block per contrast. Unlike con1DSummary(), it does the calculation for you, so you do not need to run anything first: for an lmer model it uses lmerTest's contest1D() (which works for lmer models only), and for an lme or gls model, which contest1D() cannot take, it computes each contrast directly from the model's fixed effects and covariance matrix, with the same Satterthwaite df method that lmeSummary() uses, and gives the same columns. To test more than one contrast at once (2 or more df), use conMDSummary() instead.

Basic usage

con1DTable(Model4c, L=list("Contextual SES"=c(WCses=-1, CMses=1),
                           "Contextual Nonwhite"=c(WCnon=-1, CMnon30=1)),
           CI=TRUE)

Requested Fixed Effects (with 95% CI):
 Name                  Est    SE      DF t.value p.value Lower Upper
 Contextual SES      3.375 0.417 172.761   8.090   <.001 2.551 4.198
 Contextual Nonwhite 1.351 0.594 199.136   2.274   0.024 0.179 2.523

Arguments

The printed table (headed "Requested Fixed Effects") has one row per contrast, with columns Name, Est, SE, df, t.value, p.value, plus Lower/Upper (if CI=TRUE) and Cohen.d/Partial.r (if effectSizes=TRUE)—everything con1DSummary() would print for one contrast, combined onto a single row per contrast instead of separate blocks.

Examples: contextual effects, written by position or by name (Example 4a)

The contextual effect of SES can be written by position, c(0, -1, 0, 1, 0), or by name, c(WCses=-1, CMses=1), with every effect you do not name set to 0. A named vector does not depend on the order of the effects and is easier to check. The first two rows below are the same contrast; the third is the contextual effect of the nonwhite predictor, written by name:

con1DTable(Model4c, L=list("Contextual SES (by position)"  = c(0, -1, 0, 1, 0),
                          "Contextual SES (by name)"      = c(WCses=-1, CMses=1),
                          "Contextual Nonwhite (by name)" = c(WCnon=-1, CMnon30=1)))

Requested Fixed Effects:
 Name                            Est    SE      DF t.value p.value
 Contextual SES (by position)  3.375 0.417 172.761   8.090   <.001
 Contextual SES (by name)      3.375 0.417 172.761   8.090   <.001
 Contextual Nonwhite (by name) 1.351 0.594 199.136   2.274   0.024

The same call works unchanged for the same model fit with lme() (lme(math ~ 1 + WCses + WCnon + CMses + CMnon30, random = ~ 1 | schoolID)); the estimates and SEs are identical and only the df differ slightly (174.787 and 201.468), since nlme’s variance estimates are approximate.

If the level-1 predictors were instead grand-mean-centered (rather than cluster-mean-centered), the cluster-mean slope would already be the contextual effect, and no contrast would be needed.

emmSat()

A thin wrapper for emmeans::emmeans(), for getting model-predicted means or pairwise comparisons directly from an lmer, lme, or gls model. emmeans' own default denominator df differ by model type and cannot be changed with an option: containment df for an lme model, and Kenward-Roger df (if the pbkrtest package is installed) for an lmer model. emmSat() sets the df method for you, so the means, standard errors, and tests match the df that lmerSummary(), lmeSummary(), glsSummary(), con1DTable(), and conMDSummary() report. Everything else is called exactly as in emmeans().

Basic usage

The first call prints the model-predicted mean at each occasion (from Example 5b, a gls model with an unstructured R matrix); the second adds all pairwise differences between occasions, unadjusted; the third gives predicted means at chosen values of a numeric predictor (here from an lme model with random linear change):

emmSat(Model2b, ~occasion)
emmSat(Model2b, pairwise~occasion, adjust="none")
emmSat(Model3b, ~time, at=list(time=c(0, 3)))

Output of the first call (the means), and of the second call (which prints the same means under "emmeans:" and then the pairwise differences under "contrasts:"):

 occasion emmean    SE      df lower.CL upper.CL
 1         8.225 0.237 206.548    7.758    8.692
 2         9.080 0.230 194.243    8.626    9.534
 3         9.950 0.239 199.421    9.479   10.421
 4        11.070 0.271 200.661   10.536   11.604
Degrees-of-freedom method: appx-satterthwaite
Confidence level used: 0.95

contrasts:
 contrast              estimate    SE      df t.ratio p.value
 occasion1 - occasion2   -0.855 0.207 199.035  -4.126   <.001
 occasion1 - occasion3   -1.725 0.239 199.128  -7.230   <.001
 occasion1 - occasion4   -2.845 0.277 198.292 -10.257   <.001
 occasion2 - occasion3   -0.870 0.222 197.915  -3.922   <.001
 occasion2 - occasion4   -1.990 0.256 197.932  -7.765   <.001
 occasion3 - occasion4   -1.120 0.225 198.138  -4.973   <.001
Degrees-of-freedom method: appx-satterthwaite

Arguments

The tables are always printed, in the same layout and with the same number formatting as the *Summary() functions, followed by emmeans' own notes about the df method, the confidence level, and any p-value adjustment; a pairwise~ specification prints two tables, labeled "emmeans:" and "contrasts:". The usual emmeans object is returned invisibly, so assign it (e.g. Means = emmSat(Model2b, ~occasion)) if you want to use summary(), confint(), pairs(), contrast(), or test() on it as with emmeans(). For a plain linear combination of fixed effects with Satterthwaite df, con1DTable() is an alternative; for a joint test, use conMDSummary().

conMDSummary()

Tests a set of a model's fixed effects jointly, in one of two ways. With terms, it tests whether a chosen subset of the fixed slopes are all zero—for example, just the two SES-related slopes in a model that also has other predictors—using the same mechanism each *Summary() function's own "Overall F-test/Wald test for All Fixed Slopes" already uses for the whole model, just restricted to the slopes named in terms. With L, it tests any set of linear combinations of the fixed effects equal to zero—for example, that a within-cluster slope and a between-cluster slope are equal, that two group means differ, or several such questions at once. Works with a model fit by lmer(), lme(), gls(), or glmer(): lmerTest's contestMD() for an lmer model (Satterthwaite or Kenward-Roger df); for an lme model, an F-test with Satterthwaite denominator df (the same per-coefficient df lmeSummary()'s Fixed Effects table uses, extended to several df at once the way contestMD() does it) or, with ddf="asymptotic", a multivariate Wald chi-square test; for a gls model, the same Satterthwaite F-test used for lme (which accounts for an estimated residual covariance matrix such as an unstructured R matrix) or, with ddf="residual", an F-test on its residual df, or, with ddf="asymptotic", the Wald chi-square; and always a multivariate Wald chi-square test for a glmer model (since neither contestMD() nor anova(Terms=) supports glmer). Because it does not depend on anova(Terms=), it can test slopes at different levels of a model together (e.g. a twin-level and a pair-level predictor), which nlme's own anova() cannot. Passing every one of a model's slope names reproduces the same result as that model's own *Summary() omnibus test, since that is just the special case of testing all of them at once.

Basic usage

The examples for con1DSummary(), con1DTable(), and conMDSummary() all use the same model, Model4c from Example 4a: math ~ 1 + WCses + WCnon + CMses + CMnon30 + (1|schoolID). Each level-1 predictor enters as a cluster-mean-centered within-cluster slope (WCses, WCnon) and a cluster-mean between-cluster slope (CMses, CMnon30), so its fixed effects, in order, are Intercept, WCses, WCnon, CMses, and CMnon30. Its contextual effects (between minus within) are the contrasts tested below.

conMDSummary(Model4c, terms=c("WCses","CMses"),
             name="Joint Test: SES (Within + Between)")
conMDSummary(Model4c, L=c(WCses=1, CMses=-1), name="Within = Between SES slope")
conMDSummary(Model4c, L=list("Contextual SES"=c(WCses=-1, CMses=1),
                             "Contextual Nonwhite"=c(WCnon=-1, CMnon30=1)),
             name="Both contextual effects zero")
conMDSummary(Model4c, L=list("Contextual SES"=c(WCses=-1, CMses=1), "CMnon30"),
             name="Contextual SES and CMnon30 slope")

Joint Test: SES (Within + Between)
Joint Test for:
WCses, CMses
 Model   NumDF   DenDF F.value p.value
 Model4c     2 291.926 248.279   <.001

Within = Between SES slope
Joint Test of 1 Linear Combination of the Fixed Effects
 Model   NumDF   DenDF F.value p.value
 Model4c     1 172.761  65.453   <.001

Both contextual effects zero
Joint Test of 2 Linear Combinations of the Fixed Effects:
Contextual SES, Contextual Nonwhite
 Model   NumDF   DenDF F.value p.value
 Model4c     2 193.242  33.816   <.001

Contextual SES and CMnon30 slope
Joint Test of 2 Linear Combinations of the Fixed Effects:
Contextual SES, CMnon30
 Model   NumDF   DenDF F.value p.value
 Model4c     2 163.541  59.790   <.001

The first test asks whether SES matters at all (the within and between slopes both zero). The second, a single combination, tests that the within and between SES slopes are equal; this is the same question as the contextual effect of SES being zero, so its F (65.453) is the square of that contrast’s t (8.090) in con1DSummary() and con1DTable(), with the same denominator df (172.761). The third tests both contextual effects at once (2 df), and the fourth shows how a combination and a plain fixed-effect name can be listed together. The same L works unchanged for the lme() version of this model (for example, the two contextual effects give F = 33.816 with denominator df 195.506).

Arguments

The printed output starts with a header line saying what was tested ("Joint Test for:" with terms, or how many linear combinations with L), followed on its own line by the list of terms (or the combinations' labels, if L has any), then a single-row table: Model, NumDF, DenDF, F.value, p.value for an F-test, or Model, NumDF, Chisq, p.value for a Wald chi-square test (which has no DenDF), followed by a short note whenever a Wald chi-square was used instead of an F-test (and, for a gls model, a note stating whether Satterthwaite or residual denominator df were used). A typo'd or invalid term name, or a term not found in the model, stops with a plain-English error listing the model's valid term names.

PseudoR2()

Computes pseudo-R2 for each variance component (e.g., the random intercept variance, a random slope variance, the residual variance)—the percent reduction in that variance when moving from a baseline model to one or more larger models with more predictors. This shows how much of each source of variability your predictors explain. Up to four alternative models can be compared at once, e.g., baseline → add level-1 predictors → add level-2 predictors → add a cross-level interaction—with the change in pseudo-R2 shown at each additional step. Works with models fit by lmer() or by lme() (every model in one call must be the same type, not a mix). For an lme model with weights=varIdent() (separate residual variances by group), the residual gets one row per group, labeled "Residual (<group>)", each compared to the same group in the baseline model.

Basic usage

PseudoR2(baseModel=Model1, altModel1=Model2, altModel2=Model3)

For a model with more than one grouping factor (crossed random effects, as in Example 6a, or three or more nested levels), each term is labeled with its grouping factor, so that rows with the same name, such as the two random intercepts, can be told apart. Nested factors are listed from the highest level down.

PseudoR2(baseModel=Model1b, altModel1=Model2)

Pseudo-R2
 Term                   Model   Variance Pseudo.R2
 Intercept | PschoolID  Model1b    0.175
 Intercept | PschoolID  Model2     0.172     0.015
 Intercept | SschoolID  Model1b    0.079
 Intercept | SschoolID  Model2     0.073     0.079
 Residual               Model1b    0.505
 Residual               Model2     0.504     0.002

Arguments

TotalR2()

Computes total-R2—the proportion of the total outcome variance explained by a model's fixed effects only—for a single model, or the change in total-R2 across up to four models with progressively more fixed effects. Works with models fit by lmer() or lme().

Basic usage

TotalR2(Model2, Model3, Model4)

Arguments

randLRT()

Computes a likelihood ratio test (LRT) between two models you have already fit separately, when they cannot be compared with ranova() or the LRT= option in lmerSummary()/lmeSummary()—most commonly, a single-level model (glm or gls, with no random effects) versus a nested multilevel model (glmer, lmer, or lme) that adds one or more random effects. This is the manual comparison used for GLMMs (binary/count outcomes), such as glm vs. glmer for a logistic model. Unlike a fixed df=1 test, randLRT() computes the test's degrees of freedom automatically as the difference in the number of estimated parameters between the two models, so it works correctly even if the more complex model adds more than one random-effect parameter at once. The output table includes both the ordinary p-value and the p.mix mixture-corrected p-value (see the glossary above), which is the more appropriate one when fullModel adds a variance component to baseModel.

Basic usage

SingleNon = glm(data=Example3, family=binomial(link="logit"), formula=nonwhite~1)
EmptyNon = glmer(data=Example3, family=binomial(link="logit"), nAGQ=7,
                 formula=nonwhite~1+(1|schoolID))
randLRT(SingleNon, EmptyNon, name="LRT for Random Intercept (Nonwhite)")

Arguments

glmerRanova()

Mimics what ranova() does for lmer models—testing each random effect and reporting a likelihood ratio test—but for glmer models, which ranova() does not support. You give it just the one fitted glmer model; it automatically figures out every random effect in it and refits a whole sequence of reduced models (reusing the same data, family, link, and nAGQ as the original model), so no manual nested-model fitting is required. Each row's own "Random Effect Dropped" label is condensed to "Term | Group" form, e.g. "x1 + correlations | schoolID".

For a random-effects term with more than one variance (e.g. a random intercept + slope), the slopes are tested one at a time, in the order they appear in the model's formula, each one compared to a cumulative model built from the intercept and every slope already tested: the first slope's variance plus its correlation with the intercept (2 df), the second slope's variance plus its correlations with the intercept and the first slope (3 df), and so on. As in lmerSummary()'s ranova() table, the intercept of a term that has slopes gets no row of its own. A term with only one variance (e.g. a simple random intercept with no slope) is tested by dropping it entirely, the same way ranova() does.

Example

For a model with a random intercept and two correlated random slopes, formula(1+x1+x2|schoolID), glmerRanova() prints two rows (the intercept of a term with slopes gets no row of its own):

Basic usage

EmptyNon = glmer(data=Example3, family=binomial(link="logit"), nAGQ=7,
                 formula=nonwhite~1+(1|schoolID))
glmerRanova(EmptyNon)

Arguments

The output table includes both the ordinary p-value and the p.mix mixture-corrected p-value (see the glossary above) for every row, since every row removes at least a variance component.

plotResiduals()

Draws plots of a model's level-1 (within-cluster) residuals so you can check the two assumptions made about them: that they are normally distributed and that they have constant variance. It is the stand-alone version of the residPlot option in lmerSummary() and lmeSummary(), with extra options for saving the plots to a file and for choosing which kind of residual to plot. It works for lmer and lme models. There are three plots: a histogram of the residuals with a normal curve laid over it; a normal Q-Q plot (points close to the line = approximately normal); and the residuals against the fitted values, with a zero line and a smoothed trend line (an even band around zero = constant variance; a funnel or curve = a problem). The residuals are the model's own conditional level-1 residuals, so the fitted values include the random effects.

If you give by, a fourth plot is added—boxplots of the residuals at each level of that variable (similar box heights = constant variance across the levels)—and the residuals' N, mean, and SD at each level are printed. This is the natural check for repeated-measures data, where the question is whether the variance changes over time.

Basic usage

plotResiduals(Model2)                                    # draws in the plot window
plotResiduals(Model2, by="time_F")                       # adds boxplots of the residuals at each time
plotResiduals(Model2, by="time_F", file="Residuals.png") # saves a png instead of drawing it

Arguments

The function also returns, invisibly, a data frame of each residual, its fitted value, and its by value, in case you want to look at them further.

marginalV()

Builds the model-implied marginal variance-covariance matrix (V) of one cluster's repeated outcomes, and the corresponding marginal correlation matrix, at time values you choose, from a fitted lmer or lme model. V = Z G Z' + R, where G is the random effects variance-covariance matrix, Z holds the random effects' predictors at the chosen times (a column of 1s for the intercept, and the time values for a random time slope), and R is the residual covariance matrix: the residual variance with no correlation for an lmer model, and, for an lme model, also its residual correlation structure (corAR1, corCAR1, corCompSymm, corSymm, or a spatial structure) and any variance function (varIdent, varPower, and so on), evaluated at the chosen times.

The times do not have to be times that anyone was actually observed (unlike nlme's getVarCov(type="marginal"), which is for an individual who was), so you can see how the model implies that the variances and correlations change across, or between, any times you like. Models with more than one clustering variable (for example, occasions nested in persons nested in schools) are not supported, and neither are glmer models.

Basic usage

From Example 5b: random linear change over time, with a continuous-time AR1 (corCAR1) residual correlation.

Model3c = lme(data=Example5b, method="REML", fixed=dep~1+time,
              random=~1+time|PersonID,
              correlation=corCAR1(form=~time|PersonID))
marginalV(Model3c, times=c(0, 1, 2, 3))
marginalV(Model3c, times=c(0, 0.5, 3), digits=2)

Output of the first call (the second prints the same two matrices for those three times):

Model-Implied Marginal Variance-Covariance Matrix (V) for Model3c:
       0      1      2      3
0 11.229  6.340  5.372  4.811
1  6.340 10.834  6.578  6.242
2  5.372  6.578 11.705  8.081
3  4.811  6.242  8.081 13.840

Model-Implied Marginal Correlation Matrix:
      0     1     2     3
0 1.000 0.575 0.469 0.386
1 0.575 1.000 0.584 0.510
2 0.469 0.584 1.000 0.635
3 0.386 0.510 0.635 1.000

(Rows and columns are values of time)

Arguments

The matrices are also returned, invisibly, as V and Corr in a list that includes the pieces V is built from: G (the random effects variance-covariance matrix), R (the residual covariance matrix), Z (the random effects' predictors at the chosen times), and times (the occasions used). For example, x = marginalV(Model3c, times=0:3) lets you use x$V or x$Corr.