stats (version 3.6.2)

# sigma: Extract Residual Standard Deviation 'Sigma'

## Description

Extract the estimated standard deviation of the errors, the “residual standard deviation” (misnamed also “residual standard error”, e.g., in summary.lm()'s output, from a fitted model).

Many classical statistical models have a scale parameter, typically the standard deviation of a zero-mean normal (or Gaussian) random variable which is denoted as $$\sigma$$. sigma(.) extracts the estimated parameter from a fitted model, i.e., $$\hat\sigma$$.

## Usage

sigma(object, ...)# S3 method for default
sigma(object, use.fallback = TRUE, ...)

## Arguments

object

an R object, typically resulting from a model fitting function such as lm.

use.fallback

logical, passed to nobs.

potentially further arguments passed to and from methods. Passed to deviance(*, ...) for the default method.

## Value

typically a number, the estimated standard deviation of the errors (“residual standard deviation”) for Gaussian models, and---less interpretably---the square root of the residual deviance per degree of freedom in more general models. In some generalized linear modelling (glm) contexts, $$sigma^2$$ (sigma(.)^2) is called “dispersion (parameter)”. Consequently, for well-fitting binomial or Poisson GLMs, sigma is around 1.

Very strictly speaking, $$\hat{\sigma}$$ (“$$\sigma$$ hat”) is actually $$\sqrt{\widehat{\sigma^2}}$$.

For multivariate linear models (class "mlm"), a vector of sigmas is returned, each corresponding to one column of $$Y$$.

## Details

The stats package provides the S3 generic and a default method. The latter is correct typically for (asymptotically / approximately) generalized gaussian (“least squares”) problems, since it is defined as


sigma.default <- function (object, use.fallback = TRUE, ...)
sqrt( deviance(object, ...) / (NN - PP) )


where NN <- nobs(object, use.fallback = use.fallback) and PP <- sum(!is.na(coef(object))) -- where in older R versions this was length(coef(object)) which is too large in case of undetermined coefficients, e.g., for rank deficient model fits.

deviance, nobs, vcov.

## Examples

# NOT RUN {
## -- lm() ------------------------------
lm1 <- lm(Fertility ~ . , data = swiss)
stopifnot(all.equal(sigma(lm1), summary(lm1)$sigma, tol=1e-15)) ## -- nls() ----------------------------- DNase1 <- subset(DNase, Run == 1) fm.DN1 <- nls(density ~ SSlogis(log(conc), Asym, xmid, scal), DNase1) sigma(fm.DN1) # ~= 0.01919 as from summary(..) stopifnot(all.equal(sigma(fm.DN1), summary(fm.DN1)$sigma, tol=1e-15))

# }
# NOT RUN {
<!-- % example from ./predict.glm.R -->
# }
# NOT RUN {
## -- glm() -----------------------------
## -- a) Binomial -- Example from MASS
ldose <- rep(0:5, 2)
numdead <- c(1, 4, 9, 13, 18, 20, 0, 2, 6, 10, 12, 16)
sex <- factor(rep(c("M", "F"), c(6, 6)))
sigma(budworm.lg <- glm(SF ~ sex*ldose, family = binomial))

## -- b) Poisson -- from ?glm :
## Dobson (1990) Page 93: Randomized Controlled Trial :
counts <- c(18,17,15,20,10,20,25,13,12)
outcome <- gl(3,1,9)
treatment <- gl(3,3)
sigma(glm.D93 <- glm(counts ~ outcome + treatment, family = poisson()))
## (currently) *differs* from
summary(glm.D93)$dispersion # == 1 ## and the *Quasi*poisson's dispersion sigma(glm.qD93 <- update(glm.D93, family = quasipoisson())) sigma (glm.qD93)^2 # 1.282285 is close, but not the same summary(glm.qD93)$dispersion # == 1.2933

## -- Multivariate lm() "mlm" -----------
utils::example("SSD", echo=FALSE)
sigma(mlmfit) # is the same as {but more efficient than}
sqrt(diag(estVar(mlmfit)))
# }