MixedModelsExtras.jl API

r2(model::StatisticalModel)
r²(model::StatisticalModel)

Coefficient of determination (R-squared).

For a linear model, the R² is defined as $ESS/TSS$, with $ESS$ the explained sum of squares and $TSS$ the total sum of squares.

r2(model::StatisticalModel, variant::Symbol)
r²(model::StatisticalModel, variant::Symbol)

Pseudo-coefficient of determination (pseudo R-squared).

For nonlinear models, one of several pseudo R² definitions must be chosen via variant. Supported variants are:

• :MacFadden (a.k.a. likelihood ratio index), defined as $1 - \log (L)/\log (L_0)$;
• :CoxSnell, defined as $1 - (L_0/L)^{2/n}$;
• :Nagelkerke, defined as $(1 - (L_0/L)^{2/n})/(1 - L_0^{2/n})$.
• :devianceratio, defined as $1 - D/D_0$.

In the above formulas, $L$ is the likelihood of the model, $L_0$ is the likelihood of the null model (the model with only an intercept), $D$ is the deviance of the model (from the saturated model), $D_0$ is the deviance of the null model, $n$ is the number of observations (given by nobs).

The Cox-Snell and the deviance ratio variants both match the classical definition of R² for linear models.

r2(model::LinearMixedModel; conditional=false)
r²(model::LinearMixedModel; conditional=false)

Coefficient of determination (R-squared).

R² is very non trivial for mixed models for a number of reasons and the measures here are particularly simplistic. For conditional=true, the Pearson correlation between the fitted and observed values is simply squared, in line with the relationship between the coefficient of determination and the Pearson correlation for classical OLS models. For conditional=false, new predicted values are generated based on only the fixed-effects estimates and the correlation between these predictions and the observed values is again squared.

In both cases, it is important to note that despite the usual naming convention (R² for the coefficient of determination and r for Pearson's correlation coefficient), these quantitaties are not defined in terms of the other. Instead, the correlation is a standardized measure of covariance and the coefficient of determination is measured relative to a null model. (The total sum of squares is in some sense the squared residual error of the null model.) Even generalizations of the coefficient of determination that define this value in terms of the likelihoods of an intercept-only model (i.e. the sample mean) and the likelihood of the fitted model fail to generalize to the mixed-model case: should "intercept-only" mean truly just the fixed effect (i.e. comparing the likelihoods of a classical OLS and a mixed-effects model)? or should it also include an random intercept for each grouping variable? What is the null model in this case, when determining whether each grouping variable contributes to overall fit? There is no single, clear, agreed-upon answer for these questions. The bottom line is that there are many potential ways to define a coefficient of determination for linear mixed models (and even more for the generalzed case which combines all the problems of R² for GLM and for LMM) and none of them have all the properties of the coefficient of determination for classical OLS regression.

There are more advanced approximations (see e.g. Nakagawa and colleagues' 2013 2017 papers) but the fundamental philosophical issues above remain. In the future, additional approximations may be supported, but this is not a high priority. Pull requests are welcome.

For more information, see the GLMM FAQ

adjr2(model::StatisticalModel)
adjr²(model::StatisticalModel)

Adjusted coefficient of determination (adjusted R-squared).

For linear models, the adjusted R² is defined as $1 - (1 - (1-R^2)(n-1)/(n-p))$, with $R^2$ the coefficient of determination, $n$ the number of observations, and $p$ the number of coefficients (including the intercept). This definition is generally known as the Wherry Formula I.

adjr2(model::StatisticalModel, variant::Symbol)
adjr²(model::StatisticalModel, variant::Symbol)

Adjusted pseudo-coefficient of determination (adjusted pseudo R-squared). For nonlinear models, one of the several pseudo R² definitions must be chosen via variant. The only currently supported variants are :MacFadden, defined as $1 - (\log (L) - k)/\log (L0)$ and :devianceratio, defined as $1 - (D/(n-k))/(D_0/(n-1))$. In these formulas, $L$ is the likelihood of the model, $L0$ that of the null model (the model including only the intercept), $D$ is the deviance of the model, $D_0$ is the deviance of the null model, $n$ is the number of observations (given by nobs) and $k$ is the number of consumed degrees of freedom of the model (as returned by dof).

adjr2(model::LinearMixedModel; conditional=false)
adjr²(model::LinearMixedModel; conditional=false)

Adjusted coefficient of determination (adjusted R-squared).

For linear models, the adjusted R² is defined as $(1 - (1-R^2)(n-1)/(n-p))$, with $R^2$ the coefficient of determination, $n$ the number of observations, and $p$ the number of coefficients (including the intercept). This definition is generally known as the Wherry Formula I.

For mixed-models, we use this same adjustment for the R² value computed with r2. conditional whether or not the random effects are taken into account when computing the quality of the model fit. They are nonetheless included in the adjustment (i.e. as number $p$ of parameters).

R² is not a clearly defined concept for mixed models, nor does it have all the properties typically expected of the coefficient of determination. See r2 for more information.

icc(model::MixedModel, [groups])
icc(boot::MixedModelBootstrap, [family], [groups])

Compute the intra-class correlation coefficient (ICC) for a mixed model.

The ICC is defined as the variance attributable to the groups divided by the total variance from all groups and the observation-level (residual) variance. In other words, the ICC can be interpreted as the proportion of the variance explainable by the grouping/nesting structure.

A single group can be specified as a Symbol, e.g. :subj or a number of groups can be specified as an array: [:subj, :item]. If no groups are specified, then all grouping variables are used.

When a MixedModelBootstrap is passed, a vector of ICC values for each bootstrap iteration is returned. Because MixedModelBootstrap does not store the associated model family for generalized linear mixed models, the family must be specified (e.g., Bernoulli(), Poisson()). A shortest coverage interval can be computed with MixedModels.shortestcovint.

vif(m::RegressionModel)

Compute the variance inflation factor (VIF).

The VIF measures the increase in the variance of a parameter's estimate in a model with multiple parameters relative to the variance of a parameter's estimate in a model containing only that parameter.

See also gvif.

termnames(model::StatisticalModel)

Return the names of terms used in the formula of model.

This is a convenience method for termnames(formula(model)), which returns a two-tuple of termnames applied to the left and right hand sides of the formula.

For RegressionModels with only continuous predictors, this is the same as (responsename(model), coefnames(model)) and coefnames(formula(model)).

For models with categorical predictors, the returned names reflect the variable name and not the coefficients resulting from the choice of contrast coding.

See also coefnames.

termnames(t::FormulaTerm)

Return a two-tuple of termnames applied to the left and right hand sides of the formula.

julia> termnames(@formula(y ~ 1 + x * y + (1+x|g)))
("y", ["1", "x", "y", "x & y", "(1 + x) | g"])

Similarly, formulae with an implicit intercept will not have a "1" in their variable names, because the implicit intercept does not exist until apply_schema is called (and may not exist for certain model contexts).

julia> termnames(@formula(y ~ x * y + (1+x|g)))
("y", ["x", "y", "x & y", "(1 + x) | g"])

termnames(term::AbstractTerm)

Return the name of the statistical variable associated with a term.

Return value is either a String, an iterable of Strings or the empty vector if there is no associated variable (e.g. termnames(InterceptTerm{false}())).

gvif(m::RegressionModel; scale=false)

Compute the generalized variance inflation factor (GVIF).

If scale=true, then each GVIF is scaled by the degrees of freedom for (number of coefficients associated with) the predictor: $GVIF^(1 / (2*df))$.

The GVIF measures the increase in the variance of a (group of) parameter's estimate in a model with multiple parameters relative to the variance of a parameter's estimate in a model containing only that parameter. For continuous, numerical predictors, the GVIF is the same as the VIF, but for categorical predictors, the GVIF provides a single number for the entire group of contrast-coded coefficients associated with a categorical predictor.

See also vif.

References

Fox, J., & Monette, G. (1992). Generalized Collinearity Diagnostics. Journal of the American Statistical Association, 87(417), 178. doi:10.2307/2290467

partial_fitted(model::MixedModel,
               fe::AbstractVector{<:AbstractString},
               re::Dict{Symbol}=Dict(fn => fe for fn in fnames(model));
               mode=:include,
               type=:linpred)

Compute "partial" fitted values.

Partial fitted values are useful for computing partial residuals. They are the fitted values obtained by setting selected model terms to zero, while preserving the other values at their original estimates.

The fixed effects coefficients to use (fe) are specified as vector of strings. For specifying no coefficients, use the empty string vector String[].

The random effects re are specified as a dictionary, with the grouping variables as keys and the vector of group-level coefficients specified as vectors. For example, Dict(:subj => ["(Intercept)"]) specifies that (1|subj) should be kept. The default is to match the specified fixed effects for all grouping variables, but note that this will fail when the fixed effects specification is incompatible with any grouping variable.

The keyword argument mode specifies whether the fixed and random effects specifications are treated as coefficients to :include or :exclude.

For GeneralizedLinearMixedModel, the keyword argument type specifies whether the predictions should be returned on the scale of linear predictor (:linpred) or on the response scale (:response).

This functionality is similar to the remef package in R.

shrinkagenorm(m::MixedModel{T},
              θref::AbstractVector{T}=(isa(m, LinearMixedModel) ? 1e4 : 1) .*
                                      m.optsum.initial;
              uscale=false, p=2)

Returns a NamedTuple of Tables.jl-tables norm of the change from OLS estimates (across all relevant coefficients) to BLUPs from the mixed model.

p corresponds to the $L_p$ norms, i.e. $p=2$ is the Euclidean metric.

Each entry in the named tuple corresponds to a single grouping term.

shrinkagetables(m::MixedModel{T},
                θref::AbstractVector{T}=(isa(m, LinearMixedModel) ? 1e4 : 1) .*
                                        m.optsum.initial;
                uscale=false) where {T}

Returns a NamedTuple of Tables.jl-tables of the change from OLS estimates to BLUPs from the mixed model.

Each entry in the named tuple corresponds to a single grouping term.