| sigma {stats} | R Documentation |
Extract the estimated standard deviation of the errors, the
“residual standard deviation” (misnomed 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(.) extracts the estimated parameter from a fitted
model, i.e., sigma^.
sigma(object, ...) ## Default S3 method: sigma(object, use.fallback = TRUE, ...)
object |
an R object, typically resulting from a model fitting
function such as |
use.fallback |
logical, passed to |
... |
potentially further arguments passed to and from
methods. Passed to |
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.
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”) is actually √(hat(σ^2)).
For multivariate linear models (class "mlm"), a vector
of sigmas is returned, each corresponding to one column of Y.
The misnomer “Residual standard error” has been part of too many R (and S) outputs to be easily changed there.
## -- 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))
## -- 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)))