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\).

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

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.

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\).

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.

# NOT RUN { ## -- lm() ------------------------------ lm1 <- lm(Fertility ~ . , data = swiss) sigma(lm1) # ~= 7.165 = "Residual standard error" printed from summary(lm1) 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))) SF <- cbind(numdead, numalive = 20-numdead) 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))) # }