spaMM (version 4.5.0)

residuals.HLfit: Extract model residuals

Description

Extracts several types of residuals from an object of class HLfit. Note that the default type ("deviance") of returned residuals differs from the default (response residuals) of equivalent functions in base R.

Usage

# S3 method for HLfit
residuals(object, 
  type = c("deviance", "pearson", "response", "working", "std_dev_res"), force=FALSE, ...)

Value

A vector of residuals

Arguments

object

An object of class HLfit, as returned by the fitting functions in spaMM.

type

The type of residuals which should be returned. See Details for additional information.

force

Boolean: to force recomputation of the "std_dev_res" residuals even if they are available in the object, for checking purposes.

...

For consistency with the generic.

Details

The four types "deviance" (default), "pearson", "response" are "working" are, for GLM families, the same that are returned by residuals.glm. "working" residuals may be returned only for fixed-effect models. "deviance" residuals are the signed square root of those returned by dev_resids when there are no prior weights.

In the presence of prior weights, what the standard extractors do is often a matter of confusion and spaMM has not always been consistent with them. For a gaussian-response GLM (see Examples) stats::deviance.lm calls weighted.residuals() which returns unscaled deviance residuals weighted by prior weights. Unscaled deviance residuals are defined in McCullagh and Nelder 1989, p. 34 and depend on the response values and fitted values but not on the canonical \(\phi\) parameter, and prior weights are not considered. weighted.residuals() still ignores \(\phi\) but accounts for prior weights. This means that different residuals(<glm>) and deviance(<glm>) will be returned for equivalent fits with different parametrizations of the residual variance (as produced by glm(., family=gaussian, weights=rep(2,nrow<data>)) versus the glm call without weights). residuals(<HLfit object>,"deviance") and deviance(<HLfit object>,"deviance") are consistent with this behavior. By contrast, dev_resids(<HLfit object>) always return the unscaled deviance residuals by default.

Following Lee et al. (2006, p.52), the standardized deviance residuals returned for type="std_dev_res" are defined as the deviance residuals divided by \(\phi\sqrt(1-q)\), where the deviance residuals are defined as for a GLM, \(\phi\) is the dispersion parameter of the response family (a vector of values, for heteroscedastic cases), and \(q\) is a vector of leverages given by hatvalues(., type="std") (see hatvalues for details about these specific standardizing leverages).

Some definitions must be extended for non-GLM response families. In the latter case, the deviance residuals are as defined in Details of llm.fit (there is no concept of unscaled residuals here, nor indeed of scaled ones since the residual dispersion parameter is not generally a scale factor, but the returned deviance residuals for non-GLMs are analogous to the scaled ones for GLMs as they depend on residual dispersion). "std_dev_res" residuals are defined from them as shown above for GLM response families, with the additional convention that \(\phi=1\) (since the family's own residual dispersion parameter already enters in the definition of deviance residuals for non-GLM families). Pearson residuals and response residuals are defined as in stats:::residuals.glm. The "working" residuals are defined for each response as \(- [d \log(clik)/d \eta]/[d^2 \log(clik)/d \eta^2]\) where clik is the conditional likelihood.

References

Lee, Y., Nelder, J. A. and Pawitan, Y. (2006). Generalized linear models with random effects: unified analysis via h-likelihood. Chapman & Hall: London.

Examples

data("wafers")
fit <- fitme(y ~X1+(1|batch) ,data=wafers, init=list(phi=NaN))  # : this 'init' 
#                 implies that standardized deviance residuals are saved in the 
#                 fit result, allowing the following comparison: 

r1 <- residuals(fit, type="std_dev_res") # gets stored value
r2 <- residuals(fit, type="std_dev_res", force=TRUE) # forced recomputation
if (diff(range(r1-r2))>1e-14) stop()

##### 
if (FALSE) {
glmfit <- glm(I(y/1000)~X1, family=gaussian(), data=wafers)
deviance(glmfit) #           3...                       (a)
sum(residuals(glmfit)^2) #   3...                       (b) 

# Same model, with different parametrization of residual variance 
glmfit2 <- glm(I(y/1000)~X1, family=gaussian(), data=wafers, weights=rep(2,198))
deviance(glmfit2) #          6...                       (c)  
sum(residuals(glmfit2)^2) #  6...                       (d)

# Same comparison but for HLfit objects:
spfit <- fitme(I(y/1000)~X1, family=gaussian(), data=wafers)
deviance(spfit) #            3...                       (e)
sum(residuals(spfit)^2) #    3...                       (f)  
sum(dev_resids(spfit)) #     3...                         

spfit2 <- fitme(I(y/1000)~X1, family=gaussian(), data=wafers, prior.weights=rep(2,198))
deviance(spfit2) #           6...                       (g) ~ (c,d) # post v4.2.0
sum(residuals(spfit2)^2) #   6...                       (h) ~ (c,d) 
sum(dev_resids(spfit2)) #    3...                         

# Unscaled residuals should no depend on arbitrarily fixed residual variance:
spfit3 <- fitme(I(y/1000)~X1, family=gaussian(), data=wafers, fixed=list(phi=2),
                prior.weights=rep(2,198))
deviance(spfit3) #           6...                       (i) ~ (g)
sum(residuals(spfit3)^2) #   6...                       (k) ~ (h)
sum(dev_resids(spfit3)) #    3...                         


}