Gaussian dispersion models
This function estimates Gaussian dispersion regression models.
lm.disp(formula, var.formula, data = list(), maxit = 30,
epsilon = glm.control()$epsilon, subset, na.action = na.omit,
contrasts = NULL, offset = NULL)formula |
a symbolic description of the mean function of the model to be fit. For the details of model formula specification see |
var.formula |
a symbolic description of the variance function of the model to be fit. This must be a one-sided formula; if omitted the same terms used for the mean function are used. For the details of model formula specification see |
data |
an optional data frame containing the variables in the model. By default the variables are taken from |
maxit |
integer giving the maximal number of iterations for the model fitting procedure. |
epsilon |
tolerance value for checking convergence. See |
subset |
an optional vector specifying a subset of observations to be used in the fitting process. |
na.action |
a function which indicates what should happen when the data contain |
contrasts |
an optional list as described in the |
offset |
this can be used to specify an a priori known component to be included in the linear predictor during fitting. An |
Gaussian dispersion models allow to model variance heterogeneity in Gaussian regression analysis using a log-linear model for the variance.
Suppose a response y is modelled as a function of a set of p predictors x through the linear model
y_i = β'x_i + e_i
where e_i ~ N(0, σ^2) under homogeneity.
Variance heterogeneity is modelled as
V(e_i) = σ^2 = exp(λ'z_i)
where z_i may contain some or all the variables in x_i and other variables not included in x_i; z_i is however assumed to contain a constant term.
The full model can be expressed as
E(y|x) = β'x
V(y|x) = exp(λ'z)
and it is fitted by maximum likelihood following the algorithm described in Aitkin (1987).
lm.dispmod() returns an object of class "dispmod".
The summary method can be used to obtain and print a summary of the results.
An object of class "dispmod" is a list containing the following components:
call |
the matched call. |
mean |
an object of class |
var |
an object of class |
initial.deviance |
the value of the deviance at the beginning of the iterative procedure, i.e. assuming constant variance. |
deviance |
the value of the deviance at the end of the iterative procedure. |
Based on a similar procedure available in Arc (Cook and Weisberg, http://www.stat.umn.edu/arc)
Aitkin, M. (1987), Modelling variance heterogeneity in normal regression models using GLIM, Applied Statistics, 36, 332–339.
data(minitab)
minitab <- within(minitab, y <- V^(1/3) )
mod <- lm(y ~ H + D, data = minitab)
summary(mod)
mod.disp1 <- lm.disp(y ~ H + D, data = minitab)
summary(mod.disp1)
mod.disp2 <- lm.disp(y ~ H + D, ~ H, data = minitab)
summary(mod.disp2)
# Likelihood ratio test
deviances <- c(mod.disp1$initial.deviance,
mod.disp2$deviance,
mod.disp1$deviance)
lrt <- c(NA, abs(diff(deviances)))
cbind(deviances, lrt, p.value = 1-pchisq(lrt, 1))
# quadratic dispersion model on D (as discussed by Aitkin)
mod.disp4 <- lm.disp(y ~ H + D, ~ D + I(D^2), data = minitab)
summary(mod.disp4)
r <- mod$residuals
phi.est <- mod.disp4$var$fitted.values
plot(minitab$D, log(r^2))
lines(minitab$D, log(phi.est))Please choose more modern alternatives, such as Google Chrome or Mozilla Firefox.