Calculate (Pseudo) Rsquared for a fitted model, defined here as the squared multiple correlation between the observed and fitted values for the response variable. 'Adjusted' and 'Predicted' versions are also calculated (see Details).
R2( mod, data = NULL, adj = TRUE, pred = TRUE, offset = FALSE, re.form = NULL, type = c("pearson", "spearman"), adj.type = c("olkinpratt", "ezekiel"), positive.only = TRUE, env = NULL )
mod  A fitted model object, or a list or nested list of such objects. 

data  An optional dataset, used to first refit the model(s). 
adj, pred  Logical. If 
offset  Logical. If 
re.form  For mixed models of class 
type  The type of correlation coefficient to use. Can be 
adj.type  The type of adjusted Rsquared estimator to use. Can be

positive.only  Logical, whether to return only positive values for Rsquared (negative values replaced with zero). 
env  Environment in which to look for model data (if none supplied).
Defaults to the 
A numeric vector of the Rsquared value(s), or an array, list of vectors/arrays, or nested list.
Various approaches to the calculation of a goodness of fit measure for GLMs analogous to Rsquared in the ordinary linear model have been proposed. Generally termed 'pseudo Rsquared' measures, they include variancebased, likelihoodbased, and distributionspecific approaches. Here however, a more straightforward definition is used, which can be applied to any model for which fitted values of the response variable are generated: Rsquared is calculated as the squared (weighted) correlation between the observed and fitted values of the response (in the original units). This is simply the squared version of the correlation measure advocated by Zheng & Agresti (2000), itself an intuitive measure of goodness of fit describing the predictive power of a model. As the measure does not depend on any specific error distribution or model estimating procedure, it is also generally comparable across many different types of model (Kvalseth, 1985). In the case of the ordinary linear model, the measure is exactly equal to the traditional Rsquared based on sums of squares.
If adj = TRUE
(default), the 'adjusted' Rsquared value is also returned,
which provides an estimate of the population – as opposed to sample –
Rsquared. This is achieved via an analytical formula which adjusts
Rsquared using the 'degrees of freedom' of the model (i.e. the ratio of
observations to parameters), helping to counter multiple Rsquared's
positive bias and guard against overfitting of the model to noise in the
original sample. By default, this is calculated via the exact 'OlkinPratt'
estimator, shown in recent simulations to be the optimal unbiased
population Rsquared estimator across a range of estimators and
specification scenarios (Karch, 2020), and thus a good general first
choice, even for smaller sample sizes. Setting adj.type = "ezekiel"
however will use the simpler and more common 'Ezekiel' formula, which can
be more appropriate where minimising the mean squared error (MSE) of the
estimate is more important than strict unbiasedness (Hittner, 2019; Karch,
2020).
If pred = TRUE
(default), a 'predicted' Rsquared is also returned, which
is calculated via the same formula as for Rsquared but using
crossvalidated, rather than original, fitted values. These are obtained by
dividing the model residuals (in the response scale) by the complement of
the observation leverages (diagonals of the hat matrix), then subtracting
these inflated 'predicted' residuals from the response variable. This is
essentially a short cut to obtaining 'outofsample' predictions, normally
arising via a 'leaveoneout' crossvalidation procedure (they are not
exactly equal for GLMs). The resulting Rsquared is an estimate of the
Rsquared that would result were the model fit to new data, and will be
lower than the original – and likely also the adjusted – Rsquared,
highlighting the loss of explanatory power due to sample noise. Predicted
Rsquared may be a more powerful and general indicator of overfitting than adjusted Rsquared,
as it provides a true outofsample test. This measure is a variant of an
existing one,
calculated by substituting the 'PRESS' statistic, i.e. the sum of squares
of the predicted residuals (Allen, 1974), for the residual sum of squares
in the classic Rsquared formula. It is not calculated here for GLMMs, as
the interpretation of the hat matrix is not reliable (see
hatvalues.merMod()
).
For models fitted with one or more offsets, these will be removed by
default from the response variable and fitted values prior to calculations.
Thus Rsquared will measure goodness of fit only for the predictors with
estimated, rather than fixed, coefficients. This is likely to be the
intended behaviour in most circumstances, though if users wish to include
variation due to the offset(s) they can set offset = TRUE
.
For mixed models, the function will, by default, calculate all Rsquared
measures using fitted values incorporating both the fixed and random
effects, thus encompassing all variation captured by the model. This is
equivalent to the 'conditional' Rsquared of Nakagawa et al. (2017) (though
see that reference for a more advanced approach to Rsquared for mixed
models). To include only some or no random effects, simply set the
appropriate formula using the argument re.form
, which is passed directly
to predict.merMod()
. If re.form = NA
, Rsquared is calculated for the
fixed effects only, i.e. the 'marginal' Rsquared of Nakagawa et al.
(2017).
As Rsquared is calculated here as a squared correlation, the type
of
correlation coefficient can also be specified. Setting this to "spearman"
may be desirable in some cases, such as where the relationship between
response and fitted values is not necessarily bivariate normal or linear,
and a correlation of the ranks may be more informative and/or general. This
purely monotonic Rsquared can also be considered a useful goodness of fit measure,
and may be more appropriate for comparisons between GLMs and ordinary
linear models in some scenarios.
Rsquared values produced by this function will by default be in the range
01, meaning that any negative values arising from calculations will be
converted to zero. Negative values essentially mean that the fit is 'worse'
than the null expectation of no relationship between the variables, which
can be difficult to interpret in practice and in any case usually only
occurs in rare situations, such as where the intercept is suppressed or
where a low value of Rsquared is adjusted downwards via an analytic
estimator. Such values are also 'impossible' in practice, given that
Rsquared is a strictly positive measure (as generally known). Hence, for
simplicity and ease of interpretation, values less than zero are presented
as a complete lack of model fit. This is also recommended by Shieh (2008),
who shows for adjusted Rsquared that such 'positivepart' estimators have
lower MSE in estimating the population Rsquared (though higher bias). To
allow return of negative values however, set positive.only = FALSE
. This
may be desirable for simulation purposes, and/or where strict unbiasedness
is prioritised.
Caution must be exercised in interpreting the values of any pseudo Rsquared measure calculated for a GLM or mixed model (including those produced by this function), as such measures do not hold all the properties of Rsquared in the ordinary linear model and as such may not always behave as expected. Care must also be taken in comparing the measures to their equivalents from ordinary linear models, particularly the adjusted and predicted versions, as assumptions and/or calculations may not generalise well. With that being said, the value of standardised Rsquared measures for even 'rough' model fit assessment and comparison may outweigh such reservations, and the adjusted and predicted versions in particular may aid the user in diagnosing and preventing overfitting. They should NOT, however, replace other measures such as AIC or BIC for formally comparing and/or ranking competing models fit to the same response variable.
Allen, D. M. (1974). The Relationship Between Variable Selection and Data Augmentation and a Method for Prediction. Technometrics, 16(1), 125127. doi: 10/gfgv57
Hittner, J. B. (2019). Ezekiel’s classic estimator of the population squared multiple correlation coefficient: Monte Carlobased extensions and refinements. The Journal of General Psychology, 147(3), 213–227. doi: 10/gk53wb
Karch, J. (2020). Improving on Adjusted RSquared. Collabra: Psychology, 6(1). doi: 10/gkgk2v
Kvalseth, T. O. (1985). Cautionary Note about R2. The American Statistician, 39(4), 279285. doi: 10/b8b782
Nakagawa, S., Johnson, P. C. D., & Schielzeth, H. (2017). The coefficient of determination R2 and intraclass correlation coefficient from generalized linear mixedeffects models revisited and expanded. Journal of the Royal Society Interface, 14(134). doi: 10/gddpnq
Shieh, G. (2008). Improved Shrinkage Estimation of Squared Multiple Correlation Coefficient and Squared CrossValidity Coefficient. Organizational Research Methods, 11(2), 387–407. doi: 10/bcwqf3
Zheng, B., & Agresti, A. (2000). Summarizing the predictive power of a generalized linear model. Statistics in Medicine, 19(13), 17711781. doi: 10/db7rfv
# Pseudo Rsquared for mixed models R2(shipley.sem) # fixed + random ('conditional')#> DD Date Growth Live #> R.squared 0.7080083 0.9855557 0.7938367 0.2668652 #> R.squared.adj 0.7074791 0.9855351 0.7934879 0.2655971 #> R.squared.pred 0.6835636 0.9820241 0.7552854 NAR2(shipley.sem, re.form = ~ (1  tree)) # fixed + 'tree'#> DD Date Growth Live #> R.squared 0.5305163 0.8491536 0.6068561 0.2453019 #> R.squared.adj 0.5295483 0.8489101 0.6060877 0.2439737 #> R.squared.pred 0.4898536 0.6831439 0.5012969 NAR2(shipley.sem, re.form = ~ (1  site)) # fixed + 'site'#> DD Date Growth Live #> R.squared 0.6839237 0.6925661 0.3081469 0.2016339 #> R.squared.adj 0.6833402 0.6920022 0.3065040 0.2001797 #> R.squared.pred 0.6567506 0.2541820 0.1153087 NAR2(shipley.sem, re.form = NA) # fixed only ('marginal')#> DD Date Growth Live #> R.squared 0.5012513 0.4250445 0.04814953 0.1834597 #> R.squared.adj 0.5002024 0.4237737 0.04554043 0.1819516 #> R.squared.pred 0.4578167 0.0801710 0.00000000 NAR2(shipley.sem, re.form = NA, type = "spearman") # using Spearman's Rho#> DD Date Growth Live #> R.squared 0.4043303 0.3919543 0.048572025 0.04724129 #> R.squared.adj 0.4029963 0.3905821 0.045964645 0.04529886 #> R.squared.pred 0.3008997 0.2992427 0.001204466 NA# Predicted Rsquared: compare crossvalidated predictions calculated/ # approximated via the hat matrix to standard method (leaveoneout) # Fit test models using Shipley data – compare lm vs glm d < na.omit(shipley) m < lm(Live ~ Date + DD + lat, d) # m < glm(Live ~ Date + DD + lat, binomial, d) # Manual CV predictions (leaveoneout) cvf1 < sapply(1:nrow(d), function(i) { m.ni < update(m, data = d[i, ]) predict(m.ni, d[i, ], type = "response") }) # Shortcut via the hat matrix y < getY(m) f < fitted(m) cvf2 < y  (y  f) / (1  hatvalues(m)) # Compare predictions (not exactly equal for GLMs) all.equal(cvf1, cvf2)#> [1] TRUE#> [1] 1# lm: 1; glm: 0.9999987 # NOTE: comparison not tested here for mixed models, as hierarchical data can # complicate the choice of an appropriate leaveoneout procedure. However, # there is no obvious reason why use of the leverage values (diagonals of the # hat matrix) to estimate CV predictions shouldn't generalise, roughly, to # the mixed model case (at least for LMMs). In any case, users should # exercise the appropriate caution in interpretation.