aodql {aods3}R Documentation

QL/MM Estimation of Generalized Linear Models for Overdispersed Count Data

Description

From clustered data, the function fits generalized linear models containing an over-dispersion parameter Φ\Phi using quasi-likelihood estimating equations for the mean μ\mu and a moment estimator for Φ\Phi.

For binomial-type models, data have the form {(n1,m1),(n2,m2),...,(nN,mN)(n_1, m_1), (n_2, m_2), ..., (n_N, m_N)}, where nin_i is the size of cluster ii, mim_i the number of “successes”, and NN the number of clusters. The response is the proportion y=m/ny = m/n.

For Poisson-type models, data can be of two forms. When modeling “simple counts”, data have the form {m1,m2,...,mNm_1, m_2, ..., m_N}, where mim_i is the number of occurences of the event under study. When modeling rates (e.g. hazard rates), data have the same form as for the BB model, where nin_i is the denominator of the rate for cluster ii (considered as an “offset”, i.e. a constant known value) and mim_i the number of occurences of the event. For both forms of data, the response is the count y=my = m.

Usage

  aodql(formula,
        data,
        family = c("qbin", "qpois"),
        link = c("logit", "cloglog", "probit"),
        method = c("chisq", "dev"),
        phi = NULL,
        tol = 1e-5, ...)
        
  ## S3 method for class 'aodql'
anova(object, ...)
  ## S3 method for class 'aodql'
coef(object, ...)
  ## S3 method for class 'aodql'
deviance(object, ...)
  ## S3 method for class 'aodql'
df.residual(object, ...)
  ## S3 method for class 'aodql'
fitted(object, ...)
  ## S3 method for class 'aodql'
logLik(object, ...)
  ## S3 method for class 'aodql'
predict(object, ...)
  ## S3 method for class 'aodql'
print(x, ...)
  ## S3 method for class 'aodql'
residuals(object, ...)
  ## S3 method for class 'aodql'
summary(object, ...)
  ## S3 method for class 'aodql'
vcov(object, ...)  
  

Arguments

formula

A formula for the mean μ\mu, defining the parameter vector bb (see details). For binomial-type models, the left-hand side of the formula must be of the form cbind(m, n - m) ~ ... where the fitted proportion is m/n. For Poisson-type models, the left-hand side of the formula must be of the form m ~ ... where the fitted count is m. To fit a rate, argument offset must be used in the right-hand side of the formula (see examples).

data

A data frame containing the response (m and, optionnally, n) and the explanatory variable(s).

family

Define the model which is fitted: “qbin” for binomial-type models and “qpois” for Poisson-type models.

link

For binomial-type models only. Define the link function gg for the mean μ\mu: “cloglog”, “logit” (default) or “probit”. For Poisson-type models, link is automatically set to “log”.

method

For function aodql, define the statistics used for the moment estimation of phiphi; legal values are “chisq” (default) for the chi-squared statistics or “dev” for the deviance statistics.

phi

An optional value defining the over-dispersion parameter Φ\Phi if it is set as constant. Default to NULL (in that case, Φ\Phi is estimated).

tol

A positive scalar (default to 0.001). The algorithm stops at iteration r+1r + 1 when the condition χ2[r+1]χ2[r]<=tol\chi{^2}[r+1] - \chi{^2}[r] <= tol is met by the χ2\chi^2 statistics .

...

Further arguments to passed to the appropriate functions.

object

An object of class “aodql”.

x

An object of class “aodql”.

Details

Binomial-type models

For a given cluster (n,m)(n, m), the model is

mλ,nBinomial(n,λ)m | \lambda,n \sim Binomial(n, \lambda)

where λ\lambda follows a random variable of mean E[λ]=μE[\lambda] = \mu and variance Var[λ]=Φμ(1μ)Var[\lambda] = \Phi * \mu * (1 - \mu). The marginal mean and variance of mm are

E[m]=nμE[m] = n * \mu

Var[m]=nμ(1μ)(1+(n1)Φ)Var[m] = n * \mu * (1 - \mu) * (1 + (n - 1) * \Phi)

The response in aodql is y=m/ny = m/n. The mean is E[y]=μE[y] = \mu, defined such as μ=g1(Xb)=g1(ν)\mu = g^{-1}(X * b) = g^{-1}(\nu), where gg is the link function, XX is a design-matrix, bb a vector of fixed effects and ν=Xb\nu = X * b is the corresponding linear predictor. The variance is Var[y]=(1/n)μ(1μ)(1+(n1)Φ)Var[y] = (1 / n) * \mu * (1 - \mu) * (1 + (n - 1) * \Phi).

Poisson-type models

—— Simple counts (model with no offset)

For a given cluster (m)(m), the model is

yλPoisson(λ)y | \lambda \sim Poisson(\lambda)

where λ\lambda follows a random distribution of mean μ\mu and variance Φμ2\Phi * \mu^2. The mean and variance of the marginal distribution of mm are

E[m]=μE[m] = \mu

Var[m]=μ+Φμ2Var[m] = \mu + \Phi * \mu^2

The response in aodql is y=my = m. The mean is E[y]=μE[y] = \mu, defined such as μ=exp(Xb)=exp(ν)\mu = exp(X * b) = exp(\nu). The variance is Var[y]=μ+Φμ2Var[y] = \mu + \Phi * \mu^2.

—— Rates (model with offset)

For a given cluster (n,m)(n, m), the model is

mλ,nPoisson(λ)m | \lambda,n \sim Poisson(\lambda)

where λ\lambda follows the same random distribution as for the case with no offset. The marginal mean and variance are

E[mn]=μE[m | n] = \mu

Var[mn]=μ+Φμ2Var[m | n] = \mu + \Phi * \mu^2

The response in aodql is y=my = m. The mean is E[y]=μE[y] = \mu, defined such as μ=exp(Xb+log(n))=exp(ν+log(n))=exp(η)\mu = exp(X * b + log(n)) = exp(\nu + log(n)) = exp(\eta), where log(n)log(n) is the offset. The variance is Var[y]=μ+Φμ2Var[y] = \mu + \Phi * \mu^2.

Other details

Vector bb and parameter Φ\Phi are estimated iteratively, using procedures referred to as "Model I" in Williams (1982) for binomial-type models, and "Procedure II" in Breslow (1984) for Poisson-type models.

Iterations are as follows. Quasi-likelihood estimating equations (McCullagh & Nelder, 1989) are used to estimate bb (using function glm and its weights argument), Φ\Phi being set to a constant value. Then, Φ\Phi is calculated by the moment estimator, obtained by equalizing the goodness-of-fit statistic (chi-squared X2 or deviance D) of the model to its degrees of freedom.
Parameter Φ\Phi can be set as constant, using argument phi. In that case, only bb is estimated.

Value

An object of class aodql, printed and summarized by various functions.

References

Breslow, N.E., 1984. Extra-Poisson variation in log-linear models. Appl. Statist. 33, 38-44.
Moore, D.F., 1987, Modelling the extraneous variance in the presence of extra-binomial variation. Appl. Statist. 36, 8-14.
Moore, D.F., Tsiatis, A., 1991. Robust estimation of the variance in moment methods for extra-binomial and extra-poisson variation. Biometrics 47, 383-401. McCullagh, P., Nelder, J. A., 1989, 2nd ed. Generalized linear models. New York, USA: Chapman and Hall.
Williams, D.A., 1982, Extra-binomial variation in logistic linear models. Appl. Statist. 31, 144-148.

See Also

glm

Examples


#------ Binomial-type models

data(orob2)
fm <- aodql(cbind(m, n - m) ~ seed, data = orob2, family = "qbin")
coef(fm)
vcov(fm)
summary(fm)
# chi2 tests of the seed factor in fm
wald.test(b = coef(fm), varb = vcov(fm), Terms = 2)

# chi-2 vs. deviance statistic to estimate phi
fm1 <- aodql(cbind(m, n - m) ~ seed + root, data = orob2, family = "qbin")
fm2 <- aodql(cbind(m, n - m) ~ seed + root, data = orob2, family = "qbin", method = "dev")
coef(fm1)
coef(fm2)
fm1$phi
fm2$phi
vcov(fm1)
vcov(fm2)
gof(fm1)
gof(fm2)

# estimate with fixed phi
fm <- aodql(cbind(m, n - m) ~ seed, data = orob2, family = "qbin", phi = 0.05)
coef(fm)
vcov(fm)
summary(fm)

#------ Poisson-type models

data(salmonella)
fm <- aodql(m ~ log(dose + 10) + dose, data = salmonella, family = "qpois")
coef(fm)
vcov(fm)
summary(fm)
# chi2 tests of the "log(dose + 10) + dose" factors
wald.test(b = coef(fm), varb = vcov(fm), Terms = 2:3)

# chi-2 vs. deviance statistic to estimate phi
fm1 <- aodql(m ~ log(dose + 10) + dose, data = salmonella, family = "qpois")
fm2 <- aodql(m ~ log(dose + 10) + dose, data = salmonella, family = "qpois", method = "dev")
coef(fm1)
coef(fm2)
fm1$phi
fm2$phi
vcov(fm1)
vcov(fm2)
gof(fm1)
gof(fm2)

# estimate with fixed phi
fm <- aodql(m ~ log(dose + 10) + dose, data = salmonella, family = "qpois", phi = 0.05)
coef(fm)
vcov(fm)
summary(fm)

# modelling a rate
data(dja)
# rate "m / trisk"
fm <- aodql(formula = m ~ group + offset(log(trisk)), data = dja, family = "qpois")
summary(fm)


[Package aods3 version 0.4-1.2 Index]