Version:	0.5
Title:	Analysis of Overdispersed Data using S3 Methods
Depends:	R (≥ 3.0.0), stats
Suggests:	boot, lme4
Description:	Provides functions to analyse overdispersed counts or proportions. These functions should be considered as complements to more sophisticated methods such as generalized estimating equations (GEE) or generalized linear mixed effect models (GLMM). aods3 is an S3 re-implementation of the deprecated S4 package aod.
License:	GPL-2 | GPL-3 [expanded from: GPL (≥ 2)]
Encoding:	UTF-8
NeedsCompilation:	no
Packaged:	2024-11-04 16:10:48 UTC; siberchicot
Author:	Matthieu Lesnoff [aut], Renaud Lancelot [aut], Aurélie Siberchicot [ctb, cre]
Maintainer:	Aurélie Siberchicot <aurelie.siberchicot@univ-lyon1.fr>
Repository:	CRAN
Date/Publication:	2024-11-06 15:00:02 UTC

Antibiotics against Shipping Fever in Calves

Description

Hypothetical drug trial to compare the effect of four antibiotics against Shipping fever in calves (Shoukri and Pause, 1999, Table 3.11).

Usage

data(antibio)

Format

A data frame with 24 observations on the following 3 variables.

treatment

A factor with levels 1, 2, 3 and 4

n

A numeric vector: the number of treated animals within a two-week period.

m

A numeric vector: the number of deaths at the end of the two weeks.

References

Shoukri, M.M., Pause, C.A., 1999, 2nd ed. Statistical methods for health sciences. CRC Press, London.

ML Estimation of Generalized Linear Models for Overdispersed Count Data

Description

The function fits a beta-binomial (BB) or a negative binomial (NB) generalized linear model from clustered data.

For the BB model, data have the form {(n_1, m_1), (n_2, m_2), ..., (n_N, m_N)}, where n_i is the size of cluster i, m_i the number of “successes”, and N the number of clusters. The response is the proportion y = m/n.

For the NB model, data can be of two types. When modeling simple counts, data have the form {m_1, m_2, ..., m_N}, where m_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 n_i is the denominator of the rate for cluster i (considered as an “offset”, i.e. a constant known value) and m_i the number of occurences of the event. For both types of data, the response is the count y = m.

Usage

  aodml(formula,
  data,
  family = c("bb", "nb"),
  link = c("logit", "cloglog", "probit"), 
	phi.formula = ~ 1,
  phi.scale = c("identity", "log", "inverse"),
  phi.start = NULL,
  fixpar = list(),
  hessian = TRUE,
  method = c("BFGS", "Nelder-Mead", "CG", "SANN"),
  control = list(maxit = 3000, trace = 0), ...)
  
  ## S3 method for class 'aodml'
print(x, ...)
  ## S3 method for class 'aodml'
summary(object, ...)
  ## S3 method for class 'aodml'
anova(object, ...)
  ## S3 method for class 'anova.aodml'
print(x, digits, ...)
  ## S3 method for class 'aodml'
fitted(object, ..., what = c("mu", "nu", "eta", "phi"))
  ## S3 method for class 'aodml'
residuals(object, ..., type = c("deviance", "pearson", "response"))
  ## S3 method for class 'aodml'
coef(object, ...)
  ## S3 method for class 'aodml'
df.residual(object, ...)
  ## S3 method for class 'aodml'
logLik(object, ...)
  ## S3 method for class 'aodml'
deviance(object, ...)
  ## S3 method for class 'aodml'
AIC(object, ..., k = 2)
  ## S3 method for class 'aodml'
vcov(object, ...)
  ## S3 method for class 'aodml'
predict(object, ..., type = c("link", "response"), se.fit = FALSE, newdata = NULL)
  

Arguments

Details

Beta-binomial model (BB)

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

m | \lambda,n ~ Binomial(n, \lambda)

where \lambda follows a Beta distribution Beta(a_1, a_2). Noting B the beta function, the marginal (beta-binomial) distribution of m is

P(m | n) = C(n, m) * B(a_1 + m, a_2 + n - m) / B(a_1, a_2)

Using the re-parameterization \mu = a_1 / (a_1 + a_2) and \Phi = 1 / (a_1 + a_2 + 1), the marginal mean and variance are

E[m] = n * \mu

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

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

Negative binomial model (NB)

—— Simple counts (model with no offset)

For a given cluster (m), the model is

y | \lambda ~ Poisson(\lambda)

where \lambda follows a Gamma distribution of mean \mu and shape k (Var[\lambda] = \mu^2 / k). Noting G the gamma function, the marginal (negative binomial) distribution of m is

P(m) = {G(m+k) / (m! * G(k))} * (k / (k + \mu))^k * (\mu / (k + \mu))^m

Using the re-parameterization \Phi = 1 / k, the marginal mean and variance are

E[m] = \mu

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

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

—— Rates (model with offset)

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

m | \lambda,n ~ Poisson(\lambda)

The marginal (negative binomial) distribution P(m|n) is the same as for the case with no offset (= P(m)). The response in aodml is y = m. The mean is E[y] = \mu, defined such as \mu = exp(X * b + log(n)) = exp(\nu + log(n)) = exp(\eta), where log(n) is the offset. The variance is Var[y] = \mu + \Phi * \mu^2.

Other details

Argument phi.scale of function aodml enables to estimate the over-dispersion parameter under different scales.

If phi.scale = "identity" (Default), the function estimates \Phi.

If phi.scale = "log", the function estimates log(\Phi).

If phi.scale = "inverse", the function estimates 1/\Phi.

The full parameter vector returned by aodml, say param, is equal to (b, \Phi). This vector is estimated by maximizing the log-likelihood of the marginal model using function optim. The estimated variances-covariances matrix of param is calculated by the inverse of the observed hessian matrix returned by optim, and is referred to as varparam.

Value

An object of class aodml, printed and summarized by various functions. Function deviance.aodml returns the value -2 * (logL - logL_max). The “deviance” used in function AIC.aodml to calculate AIC and AICc is -2 * logL.

References

Crowder, M.J., 1978. Beta-binomial anova for proportions. Appl. Statist. 27, 34-37.
Griffiths, D.A., 1973. Maximum likelihood estimation for the beta-binomial distribution and an application to the household distribution of the total number of cases of disease. Biometrics 29, 637-648.
Lawless, J.F., 1987. Negative binomial and mixed Poisson regression. The Canadian Journal of Statistics, 15(3): 209-225.
McCullagh, P., Nelder, J. A., 1989, 2nd ed. Generalized linear models. New York, USA: Chapman and Hall.
Prentice, R.L., 1986. Binary regression using an extended beta-binomial distribution, with discussion of correlation induced by covariate measurement errors. J.A.S.A. 81, 321-327.
Williams, D.A., 1975. The analysis of binary responses from toxicological experiments involving reproduction and teratogenicity. Biometrics 31, 949-952.

See Also

glm and optim

Examples

#------ Beta-binomial model

data(orob2)

fm1 <- aodml(cbind(m, n - m) ~ seed, data = orob2, family = "bb")
names(fm1)

# summaries
fm1
summary(fm1)
coef(fm1)
vcov(fm1)
logLik(fm1)
deviance(fm1)
AIC(fm1)
gof(fm1)

# predictions
cbind(
  fitted(fm1),
  fitted(fm1, what = "nu"),
  fitted(fm1, what = "eta"),
  fitted(fm1, what = "phi")
)
predict(fm1, type = "response", se.fit = TRUE)
newdat <- data.frame(seed = c("O73", "O75"))
predict(fm1, type = "response", se.fit = TRUE, newdata = newdat)

# model with heterogeneity in phi
fm <- aodml(cbind(m, n - m) ~ seed, data = orob2,
  family = "bb", phi.formula = ~ seed)
summary(fm)
AIC(fm1, fm)

# various phi scales
fm <- aodml(cbind(m, n - m) ~ seed, data = orob2, family = "bb")
fm$phi
fm$phi.scale
fm <- aodml(cbind(m, n - m) ~ seed, data = orob2, family = "bb",
            phi.scale = "log")
fm$phi
fm$phi.scale
fm <- aodml(cbind(m, n - m) ~ seed, data = orob2, family = "bb",
            phi.scale = "inverse")
fm$phi
fm$phi.scale

### Models with coefficient(s) set as constant

# model with 1 phi coefficient, set as constant "0.02"
fm <- aodml(formula = cbind(m, n - m) ~ seed * root, data = orob2,
  family = "bb", fixpar = list(5, 0.02))
fm$param
fm$varparam

# model with 2 phi coefficients, with the first set as constant ~ "0"
fm <- aodml(formula = cbind(m, n - m) ~ seed * root, data = orob2,
  family = "bb", phi.formula = ~ seed, fixpar = list(5, 1e-8))
fm$param
fm$varparam

# model with 2 phi coefficients, with the first set as constant ~ "0",
# and the mu intercept (1rst coef of vector b) set as as constant "-0.5"
fm <- aodml(formula = cbind(m, n - m) ~ seed * root, data = orob2,
            family = "bb", phi.formula = ~ seed,
            fixpar = list(c(1, 5), c(-0.5, 1e-8)))
fm$param
fm$varparam
  
### Model tests

# LR tests - non-constant phi
fm0 <- aodml(cbind(m, n - m) ~ 1, data = orob2, family = "bb")
fm2 <- aodml(cbind(m, n - m) ~ seed + root, data = orob2, family = "bb")
fm3 <- aodml(cbind(m, n - m) ~ seed * root, data = orob2, family = "bb")
anova(fm0, fm1, fm2, fm3)

# LR tests - constant phi
# phi is assumed to be estimated from fm3
fm2.bis <- aodml(cbind(m, n - m) ~ seed  + root, data = orob2,
                 family = "bb", fixpar = list(4, fm3$phi))
LRstat <- 2 * (logLik(fm3) - logLik(fm2.bis))  
pchisq(LRstat, df = 1, lower.tail = FALSE)  
  
# Wald test of the seed factor in fm1
wald.test(b = coef(fm3), varb = vcov(fm3), Terms = 4)

#------ Negative binomial model

### Modelling counts

data(salmonella)
fm1 <- aodml(m ~ log(dose + 10) + dose, data = salmonella, family = "nb")
## fm1 <- aodml(m ~ log(dose + 10) + dose, data = salmonella, family = "nb",
##              method = "Nelder-Mead")

# summaries
fm1
summary(fm1)
coef(fm1)
vcov(fm1)
logLik(fm1)
deviance(fm1)
AIC(fm1)
gof(fm1)

# predictions
cbind(
  fitted(fm1),
  fitted(fm1, what = "nu"),
  fitted(fm1, what = "eta"),
  fitted(fm1, what = "phi")
)
predict(fm1, type = "response", se.fit = TRUE)
newdat <- data.frame(dose = c(20, 40))
predict(fm1, type = "response", se.fit = TRUE, newdata = newdat)

# various phi scales
fm <- aodml(m ~ log(dose + 10) + dose, data = salmonella, family = "nb")
fm$phi
fm$phi.scale
fm <- aodml(m ~ log(dose + 10) + dose, data = salmonella,
            family = "nb", phi.scale = "log")
fm$phi
fm$phi.scale
fm <- aodml(m ~ log(dose + 10) + dose, data = salmonella,
            family = "nb", phi.scale = "inverse")
fm$phi
fm$phi.scale
  
# LR and Wald tests of the "log(dose + 10) + dose" factors
fm0 <- aodml(m ~ 1, data = salmonella, family = "nb")
anova(fm0, fm1)
fm0.bis <- aodml(m ~ 1, data = salmonella, family = "nb",
                 fixpar = list(2, fm1$phi))
LRstat <- 2 * (logLik(fm1) - logLik(fm0.bis))  
pchisq(LRstat, df = 2, lower.tail = FALSE)  
wald.test(b = coef(fm1), varb = vcov(fm1), Terms = 2:3)

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

fm <- aodml(formula = m ~ group + offset(log(trisk)),
            phi.formula = ~ group, data = dja, family = "nb",
            phi.scale = "log")
summary(fm)
  
# model with 1 phi coefficient, set as constant "0.8"
fm <- aodml(formula = m ~ group + offset(log(trisk)), data = dja,
            family = "nb", phi.formula = ~1, fixpar = list(3, 0.8))
fm$param
fm$varparam

# model with 2 phi coefficients, with the first set as constant ~ "0" in the identity scale
fm <- aodml(formula = m ~ group + offset(log(trisk)), data = dja,
            family = "nb", phi.formula = ~ group, phi.scale = "log",
            fixpar = list(4, -15))
fm$param
fm$varparam

# model with 2 phi coefficients, with the first set as constant "0" in the identity scale,
# and the mu intercept coefficient (1rst coef of vector b) set as as constant "-0.5"
fm <- aodml(formula = m ~ group + offset(log(trisk)), data = dja,
  family = "nb", phi.formula = ~ group, phi.scale = "log",
  fixpar = list(c(1, 4), c(-0.5, -15)))
fm$param
fm$varparam

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 {(n_1, m_1), (n_2, m_2), ..., (n_N, m_N)}, where n_i is the size of cluster i, m_i the number of “successes”, and N the number of clusters. The response is the proportion y = m/n.

For Poisson-type models, data can be of two forms. When modeling “simple counts”, data have the form {m_1, m_2, ..., m_N}, where m_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 n_i is the denominator of the rate for cluster i (considered as an “offset”, i.e. a constant known value) and m_i the number of occurences of the event. For both forms of data, the response is the count y = 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

Details

Binomial-type models

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

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

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

E[m] = n * \mu

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

The response in aodql is y = m/n. The mean is E[y] = \mu, defined such as \mu = g^{-1}(X * b) = g^{-1}(\nu), where g is the link function, X is a design-matrix, b a vector of fixed effects and \nu = X * b is the corresponding linear predictor. The variance is 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), the model is

y | \lambda \sim Poisson(\lambda)

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

E[m] = \mu

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

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

—— Rates (model with offset)

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

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[m | n] = \mu

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

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

Other details

Vector b 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 b (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 b 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)

Analysis of Overdispersed Data

Description

This package provides functions to analyse overdispersed counts or proportions. The functions should be considered as complements to more sophisticated methods such as generalized estimating equations (GEE) or generalized linear mixed effect models (GLMM).

Details

Functions Index :

Data sets Index :

Author(s)

Matthieu Lesnoff <matthieu.lesnoff@cirad.fr> and Renaud Lancelot <renaud.lancelot@cirad.fr>
Maintainer: Renaud Lancelot <renaud.lancelot@cirad.fr>

Age, Period and Cohort Effects for Vital Rates

Description

Number of prostate cancer deaths and midperiod population for nonwhites in the USA by age and period. The cohort index k is related to age and period indices (i and j, respectively) by k = j + I - i, where I = max(i) (Holford, 1983, Table 2).

Usage

data(cohorts)

Format

A data frame with 49 observations on the following 4 variables.

period

A factor with levels 1935-, 1940-, ..., 1965-.

age

A factor with levels 50-, 55-, ..., 80-.

m

Numeric: the number of prostate cancer deaths.

n

Numeric: the midperiod population size.

References

Holford, T.R., 1983. The estimation of age, period and cohort effects for vital rates. Biometrics 39, 311-324.

Mortality of Djallonke Lambs in Senegal

Description

Field trial to assess the effect of ewes deworming (prevention of gastro-intestinal parasitism) on the mortality of their offspring (age < 1 year). This data set is extracted from a large database on small ruminants production and health in Senegal (Lancelot et al., 1998). Data were collected in a sample of herds in Kolda (Upper Casamance, Senegal) during a multi-site survey (Faug?re et al., 1992). See also the references below for a presentation of the follow-up survey (Faug?re and Faug?re, 1986) and a description of the farming systems (Faug?re et al., 1990).

Usage

data(dja)

Format

A data frame with 21 observations on the following 4 variables.

group

a factor with 2 levels: CTRL and TREAT, indicating the treatment.

village

a factor indicating the village of the herd.

herd

a factor indicating the herd.

n

a numeric vector: the number of animals exposed to mortality.

trisk

a numeric vector: the exposition time to mortality (in year).

m

a numeric vector: the number of deaths.

References

Faug?re, O., Faug?re, B., 1986. Suivi de troupeaux et contr?le des performances individuelles des petits ruminants en milieu traditionnel africain. Aspects m?thodologiques. Rev. Elev. M?d. v?t. Pays trop., 39 (1): 29-40.
Faug?re, O., Dock?s, A.-C., Perrot, C., Faug?re, B., 1990. L'?levage traditionnel des petits ruminants au S?n?gal. I. Pratiques de conduite et d'exploitation des animaux chez les ?leveurs de la r?gion de Kolda. Revue Elev. M?d. v?t. Pays trop. 43: 249-259.
Faug?re, O., Tillard, E., Faug?re, B., 1992. Prophylaxie chez les petits ruminants au S?n?gal : r?gionalisation d'une politique nationale de protection sanitaire. In: B. Rey, S. H. B. Lebbie, L. Reynolds (Ed.), First biennial conference of the African Small Ruminant Research Network, ILCA, 1990, ILRAD, Nairobi, pp. 307-314.
Lancelot, R., Faye, B., Juan?s, X., Ndiaye, M., P?rochon, L., Tillard, E., 1998. La base de donn?es BAOBAB: un outil pour mod?liser la production et la sant? des petits ruminants dans les syst?mes d'?levage traditionnels au S?n?gal. Revue Elev. M?d. v?t. Pays trop., 51 (2): 135-146.

Test of Proportion Homogeneity between Groups using Donner's and Rao-Scott's Adjustments

Description

The function tests the homogeneity of probabilities between J groups (H_0: \mu_1 = \mu_2 = ... = \mu_J) from clustered binomial data {(n_1, m_1), (n_2, m_2), ..., (n_N, m_N)}, where n_i is the size of cluster i, m_i the number of “successes” (proportions are y = m/n), and N the number of clusters. The function uses adjusted chi-squared statistics, with either the correction proposed by proposed by Donner (1989) or the correction proposed by Rao and Scott (1993).

Usage

  drs(formula, data, type = c("d", "rs"), C = NULL, pooled = FALSE)
  
  ## S3 method for class 'drs'
print(x, ...)
  

Arguments

Details

Donner's test
The chi-squared statistic is adjusted with the correction factor C_j computed in each group j. The test statistic is given by:

X^2 = \sum_{j} ( (m_j - n_j * \mu)^2 / (C_j * n_j * \mu * (1 - \mu)) )

where \mu = \sum_{j} (m_j) / \sum_{j} (n_j) and C_j = 1 + (n_{A,j} - 1) * \rho. n_{A,j} is a scalar depending on the cluster sizes, and \rho is the ANOVA estimate of the intra-cluster correlation assumed common across groups (see Donner, 1989 or Donner et al., 1994). The statistic is compared to a chi-squared distribution with J - 1 degrees of freedom. Fixed correction factors C_j can be specified with the argument C.

Rao ans Scott's test
The method uses design effects and “effective” sample sizes. The design effect C_j in each group j is estimated by C_j = v_{ratio,j} / v_{bin,j}, where v_{ratio,j} is the variance of the ratio estimate of the probability in group i (Cochran, 1999, p. 32 and p. 66) and v_{bin,j} is the standard binomial variance. The C_j are used to compute the effective sample sizes n_{adj,j} = n_j / C_j, the effective numbers of successes m_{adj,j} = m_j / C_j in each group j, and the overall effective proportion mu_adj = \sum_{j} m_{adj,j} / \sum_{j} C_j. The test statistic is obtained by substituting these quantities in the usual chi-squared statistic, yielding:

X^2 = \sum_{j} ( (m_{adj,j} - n_{adj,j} * muadj)^2 / (n_{adj,j} * muadj * (1 - muadj)) )

which is compared to a chi-squared distribution with J - 1 degrees of freedom.
A pooled design effect over the J groups is estimated if argument pooled = TRUE (see Rao and Scott, 1993, Eq. 6). Fixed design effects C_j can be specified with the argument C.

Value

An object of class drs, printed with print.drs.

References

Donner, A., 1989. Statistical methods in ophthalmology: an adjusted chi-squared approach. Biometrics 45, 605-611.
Donner, A., 1993. The comparison of proportions in the presence of litter effects. Prev. Vet. Med. 18, 17-26.
Donner, A., Eliasziw, M., Klar, N., 1994. A comparison of methods for testing homogeneity of proportions in teratologic studies. Stat. Med. 13, 1253-1264.

See Also

chisq.test

Examples

data(dja)
# Donner
drs(formula = cbind(m, n - m) ~ group, data = dja, type = "d")
# Rao and Scott
drs(formula = cbind(m, n - m) ~ group, type = "rs", data = dja)
drs(formula = cbind(m, n - m) ~ group, type = "rs", data = dja, pooled = TRUE)
# standard chi2 test
drs(formula = cbind(m, n - m) ~ group, data = dja, type = "d", C = c(1:1))
drs(formula = cbind(m, n - m) ~ group, data = dja, type = "rs", C = c(1:1))

Test of Goodness-of-Fit of Models for Count data

Description

The function returns a chi-squared test of goodness of fit for models of class glm, aodml or aodql.

Usage

  gof(object)
  gof.default(object)
  
  ## S3 method for class 'gof'
print(x, ..., digits =  max(3, getOption("digits") - 3))
  

Arguments

Details

Function gof calculates the deviance D and the Pearson chi-squared X^2 statistics for the model under consideration. Let y be the observed response, and E[y] = \mu and Var[y] its mean and variance estimated from the model, statistic X^2 is calculated by:

X^2 = \sum_{i}( (y_i - \mu)^2/Var[y_i] )

Assuming that the data length is N and the number of the parameters in the model is p, D and X^2 are compared to a chi-squared distribution with N-p degrees of freedom.

Value

An object of class gof, printed with print.gof.

References

Agresti, A. Categorical data analysis. Wiley, 1990.

See Also

residuals, chisq.test

Examples

data(orob2)
fm1 <- glm(cbind(m, n - m) ~ seed, data = orob2, family = binomial)
fm2 <- aodml(cbind(m, n - m) ~ seed, data = orob2, family = "bb")
gof(fm1)
gof(fm2)

Intra-Cluster Correlation for Clustered Binomial data

Description

The function estimates the intraclass correlation \rho from clustered binomial data:

{(n_1, m_1), (n_2, m_2), ..., (n_N, m_N)},

where n_i is the size of cluster i, m_i the number of “successes” (proportions are y = m/n), and N the number of clusters. The function uses a one-way random effect model. Three estimates, corresponding to methods referred to as “A”, “B” and “C” in Goldstein et al. (2002), can be returned.

Usage

  iccbin(n, m, method = c("A", "B", "C"), nAGQ = 1, M = 1000)
  
  ## S3 method for class 'iccbin'
print(x, ...)
  

Arguments

Details

Before computations, the clustered data are split to binary “0/1” observations y_{ij} (observation j in cluster i). The methods of calculation are described in Goldstein et al. (2002).

Methods "A" and "B" use the 1-way logistic binomial-Gaussian model

y_{ij} | \mu_{ij} \sim Bernoulli(\mu_{ij})

logit(\mu_{ij}) = b_0 + u_i,

where b_0 is a constant and u_i a cluster random effect with u_i \sim N(0, s^2_u). The ML estimate of the variance component s^2_u is calculated with the function glmer of package lme4. The intra-class correlation \rho = Corr[y_{ij}, y_{ik}] is then calculated from a first-order model linearization around E[u_i]=0 in method “A”, and with Monte Carlo simulations in method “B”.

Method "C" provides the common ANOVA (moment) estimate of \rho. For details, see for instance Donner (1986), Searle et al. (1992) or Ukoumunne (2002).

Value

An object of class iccbin, printed with print.iccbin.

References

Donner A., 1986, A review of inference procedures for the intraclass correlation coefficient in the one-way random effects model. International Statistical Review 54, 67-82.
Searle, S.R., Casella, G., McCulloch, C.E., 1992. Variance components. Wiley, New York.
Ukoumunne, O. C., 2002. A comparison of confidence interval methods for the intraclass correlation coefficient in cluster randomized trials. Statistics in Medicine 21, 3757-3774.
Golstein, H., Browne, H., Rasbash, J., 2002. Partitioning variation in multilevel models. Understanding Statistics 1(4), 223-231.

See Also

glmer

Examples

data(rats)
z <- rats[rats$group == "TREAT", ]
# A: glmm (model linearization)
iccbin(z$n, z$m, method = "A")
iccbin(z$n, z$m, method = "A", nAGQ = 10)
# B: glmm (Monte Carlo)
iccbin(z$n, z$m, method = "B")
iccbin(z$n, z$m, method = "B", nAGQ = 10, M = 1500)
# C: lmm (ANOVA moments)
iccbin(z$n, z$m, method = "C")

## Not run: 
  # Example of CI calculation with nonparametric bootstrap
  require(boot)
  foo <- function(X, ind) {
    n <- X$n[ind]
    m <- X$m[ind]
    iccbin(n = n, m = m, method = "C")$rho
    }
  res <- boot(data = z[, c("n", "m")], statistic = foo, R = 500, sim = "ordinary", stype = "i")
  res
  boot.ci(res, conf = 0.95, type = "basic")
  
## End(Not run)
  

Transformation from the Link Scale to the Observation Scale

Description

The function transforms a variable from the link scale to the observation scale (probability or count).

Usage

invlink(x, type = c("cloglog", "log", "logit", "probit"))

Arguments

Value

clog-log: 1 - exp(-exp(x))
log: exp(x)
logit: exp(x) / (1 + exp(x))
probit: pnorm(x)

See Also

link

Examples

x <- seq(-5, 5, length = 100)
plot(x, invlink(x, type = "logit"), type = "l", lwd = 2, ylab = "Probability")
lines(x, invlink(x, type = "cloglog"), lty = 2, lwd = 2)
grid(col = "black")
legend(-5, 1, legend = c("alogit(x)", "acloglog(x)"), lty = c(1, 2), bg = "white")

Transformation from the Observation Scale to the Link Scale

Description

The function transforms a variable from the observation scale (probability or count) to the link scale.

Usage

link(x, type = c("cloglog", "log", "logit", "probit"))

Arguments

Value

clog-log: log(-log(1 - x ))
log: log(x)
logit: log(x / (1 - x))
probit: qnorm(x)

See Also

invlink

Examples

x <- seq(.001, .999, length = 100)
plot(x, link(x, type = "logit"), type = "l", lwd = 2, ylab = "link(proba.)")
lines(x, link(x, type = "cloglog"), lty = 2, lwd = 2)
grid(col = "black")
legend(0, 6, legend = c("logit(x)", "cloglog(x)"), lty = c(1, 2), bg = "white")

A Comparison of Site Preferences of Two Species of Lizard

Description

“These data describe the daytime habits of two species of lizards, grahami and opalinus. They were collected by observing occupied sites or perches and recording the appropriate description, namely species involved, time of the day, height and diameter of the perch and whether the site was sunny or shaded. Time of the day is recorded as early, mid-day or late.” (McCullagh and Nelder, 1989, p.129).

Usage

data(lizards)

Format

A data frame with 24 observations on the following 6 variables.

Site

A factor with levels Sun and Shade.

Diameter

A factor with levels D <= 2 and D > 2 (inches).

Height

A factor with levels H < 5 and H >= 5 (feet).

Time

A factor with levels Early, Mid-day and Late.

grahami

A numeric vector giving the observed sample size for grahami lizards.

opalinus

A numeric vector giving the observed sample size for opalinus lizards.

Details

The data were originally published in Fienberg (1970).

Source

McCullagh, P., Nelder, J. A., 1989, 2nd ed. Generalized linear models. New York, USA: Chapman and Hall.

References

Fienberg, S.E., 1970. The analysis of multidimensional contingency tables. Ecology 51: 419-433.

Pregnant Female Mice Experiment

Description

Unpublished laboratory data on the proportion of affected foetuses in two groups (control and treatment) of 10 pregnant female mice (Kupper and Haseman, 1978, p. 75).

Usage

data(mice)

Format

A data frame with 20 observations on the following 3 variables.

group

a factor with levels CTRL and TREAT

n

a numeric vector: the total number of foetuses.

m

a numeric vector: the number of affected foetuses.

References

Kupper, L.L., Haseman, J.K., 1978. The use of a correlated binomial model for the analysis of a certain toxicological experiments. Biometrics 34, 69-76.

Germination Data

Description

Data describing the germination for seed Orobanche cernua cultivated in three dilutions of a bean root extract (Crowder, 1978, Table 1). The mean proportions of the three groups are 0.142, 0.872 and 0.842, and the overall mean is 0.614.

Usage

data(orob1)

Format

A data frame with 16 observations on the following 3 variables.

dilution

a factor with 3 levels: 1/1, 1/25 and 1/625

.

n

a numeric vector: the number of seeds exposed to germination.

m

a numeric vector: the number of seeds which actually germinated.

References

Crowder, M.J., 1978. Beta-binomial anova for proportions. Appl. Statist. 27, 34-37.

Germination Data

Description

A 2 x 2 factorial experiment comparing 2 types of seed and 2 root extracts (Crowder, 1978, Table 3). There are 5 or 6 replicates in each of the 4 treatment groups, and each replicate comprises a number of seeds varying between 4 and 81. The response variable is the proportion of seeds germinating in each replicate.

Usage

data(orob2)

Format

A data frame with 21 observations on the following 4 variables.

seed

a factor with 2 levels: O73 and O75.

root

a factor with 2 levels BEAN and CUCUMBER.

n

a numeric vector: the number of seeds exposed to germination.

m

a numeric vector: the number of seeds which actually germinated.

References

Crowder, M.J., 1978. Beta-binomial anova for proportions. Appl. Statist. 27, 34-37.

Rabbits Foetuses Survival Experiment

Description

Experimental data for analyzing the effect of an increasing dose of a compound on the proportion of live foetuses affected (Paul, 1982, Table 1). Four treatment-groups were considered: control “C”, low dose “L”, medium dose “M” and high dose “H”. The animal species used in the experiment was banded Dutch rabbit.

Usage

data(rabbits)

Format

A data frame with 84 observations on the following 3 variables.

group

a factor with levels C, H, L and M

n

a numeric vector: the total number of foetuses.

m

a numeric vector: the number of affected foetuses.

References

Paul, S.R., 1982. Analysis of proportions of affected foetuses in teratological experiments. Biometrics 38, 361-370.

Rats Diet Experiment

Description

“Weil (1970) in Table 1 gives the results from an experiment comprising two treatments. One group of 16 pregnant female rats was fed a control diet during pregnancy and lactation, the diet of a second group of 16 pregnant females was treated with a chemical. For each litter the number n of pups alive at 4 days and the number x of pups that survived the 21 day lactation period were recorded.” (Williams, 1975, p. 951).

Usage

data(rats)

Format

A data frame with 32 observations on the following 3 variables.

group

A factor with levels CTRL and TREAT

n

A numeric vector: the number of pups alive at 4 days.

m

A numeric vector: the number of pups that survived the 21 day lactation.

Source

Williams, D.A., 1975. The analysis of binary responses from toxicological experiments involving reproduction and teratogenicity. Biometrics 31, 949-952.

References

Weil, C.S., 1970. Selection of the valid number of sampling units and a consideration of their combination in toxicological studies involving reproduction, teratogenesis or carcinogenesis. Fd. Cosmet. Toxicol. 8, 177-182.

Salmonella Reverse Mutagenicity Assay

Description

“Data for our third example were compiled by Margolin et al. (1981) from an Ames Salmonella reverse mutagenicity assay. Table 1 shows the number of revertant colonies observed on each of 3 replicate plates tested at each of 6 dose levels of quinoline.” (Breslow, 1984, Table 1).

Usage

data(salmonella)

Format

A data frame with 18 observations on the following 2 variables.

dose

a numeric vector: the dose level of quinoline (microgram per plate).

m

a numeric vector: the number of revertant colonies of TA98 Salmonella.

Source

Breslow, N.E., 1984. Extra-Poisson variation in log-linear models. Applied Statistics 33(1), 38-44.

References

Margolin, B.H., Kaplan, N., Zeiger, E., 1981. Statistical analysis of the Ames Salmonella / microsome test. Proc. Natl Acad. Sci. USA 76, 3779-3783.

Split Grouped Data Into Individual Data

Description

The function splits grouped binomial data and optional covariates to individual binary data. Two types of grouped data are managed by splitbin:

• Grouped data with weights;

• Grouped data of form cbind(success, failure).

When weights, successes or failures involve non-integer numbers, these numbers are rounded (using round()) before splitting.

Usage

splitbin(formula, data, id = "idbin")

Arguments

Value

A data frame built according to the formula and function used in the call.

Examples

# grouped data with weights
z <- data.frame(
    m = c(0, 1, 0, 1),
    f1 = c("A", "A", "B", "B"),
    f2 = c("C", "D", "C", "D"),
    n = c(4, 2, 1, 3)
    )
z
splitbin(formula = n ~ f1, data = z)$tab
splitbin(formula = n ~ f1 + f2 + m , data = z)$tab

# grouped data of form "cbind(success, failure)"
z <- data.frame(
    m = c(4, 1),
    n = c(5, 3),
    f1 = c("A", "B"),
    f2 = c("C", "D")
    )
z
splitbin(formula = cbind(m, n - m) ~ 1, data = z)$tab
splitbin(formula = cbind(m, n - m) ~ f1 + f2, data = z)$tab
splitbin(formula = cbind(m, n - m) ~ f1 + f2, data = z)$tab

Estimate of a Probability from Clustered Binomial Data

Description

The function estimates a probability and its variance from clustered binomial data

{(n_1, m_1), (n_2, m_2), ..., (n_N, m_N)},

where n_i is the size of cluster i, m_i the number of “successes” (proportions are y = m/n), and N the number of clusters. Confidence intervals are calculated using a normal approximation, which might be inappropriate when the probability is close to 0 or 1.

Usage

  varbin(n, m, alpha = 0.05, R = 5000)
  
  ## S3 method for class 'varbin'
print(x, ...)
  

Arguments

Details

Five methods are used for the estimations. Let us consider N clusters of sizes n_1, \ldots, n_N with observed count responses m_1, \ldots, m_N. We note y_i = m_i/n_i (i = 1, \ldots, N) the observed proportions. The underlying assumption is that the probability, say mu, is homogeneous across the clusters.

Binomial method: the probability estimate and its variance are calculated by

\mu = (sum_{i} (m_i)) / (sum_{i} (n_i)) (ratio estimate) and

\mu * (1 - \mu) / (sum_{i} (n_i) - 1), respectively.

Ratio method: the probability \mu is estimated as for the binomial method (ratio estimate). The one-stage cluster sampling formula is used to calculate the variance of \mu (see Cochran, 1999, p. 32 and p. 66).

Arithmetic method: the probability is estimated by \mu = sum_{i} (y_i) / N. The variance of \mu is estimated by sum_{i} (y_i - \mu)^2 / (N * (N - 1)).

Jackknife method: the probability is estimated by \mu defined by the arithmetic mean of the pseudovalues y_{v,i}. The variance is estimated by sum_{i} (y_{v,i} - \mu)^2 / (N * (N - 1)) (Gladen, 1977, Paul, 1982).

Bootstrap method: R samples of clusters of size N are drawn with equal probability from the initial sample (y_1, \ldots , y_N) (Efron and Tibshirani, 1993). The bootstrap estimate \mu and its estimated variance are the arithmetic mean and the empirical variance (computed with denominator R - 1) of the R binomial ratio estimates, respectively.

Value

An object of class varbin, printed with print.varbin.

References

Cochran, W.G., 1999, 3th ed. Sampling techniques. Wiley, New York.
Efron, B., Tibshirani, R., 1993. An introduction to the bootstrap. Chapman and Hall, London.
Gladen, B., 1977. The use of the jackknife to estimate proportions from toxicological data in the presence of litter effects. JASA 74(366), 278-283.
Paul, S.R., 1982. Analysis of proportions of affected foetuses in teratological experiments. Biometrics 38, 361-370.

See Also

boot

Examples

data(rabbits)
z <- rabbits[rabbits$group == "M", ]
varbin(z$n, z$m)
by(rabbits,
	list(group = rabbits$group),
  function(x) varbin(n = x$n, m = x$m, R = 1000))

Wald Test for Model Coefficients

Description

The function returns a Wald chi-squared test or a F test for a vector of model coefficients (possibly of length one), given its variance-covariance matrix.

Usage

  wald.test(b, varb, Terms = NULL, L = NULL, H0 = NULL,  df = NULL, verbose = FALSE, ...)
  
  ## S3 method for class 'wald.test'
print(x, ..., digits = max(3, getOption("digits") - 3))
  

Arguments

Details

The assumption is that the coefficients follow asymptotically a multivariate normal distribution with mean equal to the model coefficients b and variance equal to their var-cov matrix varb.

One (and only one) of Terms or L must be given. When L is given, it must have the same number of columns as the length of b, and the same number of rows as the number of linear combinations of coefficients.

When df is given, the chi-squared Wald statistic is divided by m, the number of linear combinations of coefficients to be tested (i.e., length(Terms) or nrow(L)). Under the null hypothesis H_0, this new statistic follows an F(m, df) distribution.

Value

An object of class wald.test, printed with print.wald.test.

References

Diggle, P.J., Liang, K.-Y., Zeger, S.L., 1994. Analysis of longitudinal data. Oxford, Clarendon Press, 253 p.
Draper, N.R., Smith, H., 1998. Applied Regression Analysis. New York, John Wiley & Sons, Inc., 706 p.

Examples

data(orob2)
fm <- aodql(cbind(m, n - m) ~ seed * root, data = orob2, family = "qbin")
# Wald chi2 test for the effect of root
wald.test(b = coef(fm), varb = vcov(fm), Terms = 3:4)
L <- matrix(c(0, 0, 1, 0, 0, 0, 0, 1), nrow = 2, byrow = TRUE)
wald.test(b = coef(fm), varb = vcov(fm), L = L)