R: Likelihood-Ratio Test Corrected for Overdispersion

anovaOD {AICcmodavg}

R Documentation

Likelihood-Ratio Test Corrected for Overdispersion

Description

Compute likelihood-ratio test between a given model and a simpler model.

Usage

anovaOD(mod.simple, mod.complex, c.hat = 1, 
        nobs = NULL, ...)

## S3 method for class 'glm'
anovaOD(mod.simple, mod.complex, c.hat = 1,  
        nobs = NULL, ...)

## S3 method for class 'unmarkedFitOccu'
anovaOD(mod.simple, mod.complex, c.hat = 1,  
        nobs = NULL, ...)

## S3 method for class 'unmarkedFitColExt'
anovaOD(mod.simple, mod.complex, c.hat = 1,  
        nobs = NULL, ...)

## S3 method for class 'unmarkedFitOccuRN'
anovaOD(mod.simple, mod.complex, c.hat = 1,  
        nobs = NULL, ...)

## S3 method for class 'unmarkedFitPCount'
anovaOD(mod.simple, mod.complex, c.hat = 1,  
        nobs = NULL, ...)

## S3 method for class 'unmarkedFitPCO'
anovaOD(mod.simple, mod.complex, c.hat = 1,  
        nobs = NULL, ...)

## S3 method for class 'unmarkedFitDS'
anovaOD(mod.simple, mod.complex, c.hat = 1,  
        nobs = NULL, ...)

## S3 method for class 'unmarkedFitGDS'
anovaOD(mod.simple, mod.complex, c.hat = 1,  
        nobs = NULL, ...)

## S3 method for class 'unmarkedFitOccuFP'
anovaOD(mod.simple, mod.complex, c.hat = 1,  
        nobs = NULL, ...)

## S3 method for class 'unmarkedFitMPois'
anovaOD(mod.simple, mod.complex, c.hat = 1,  
        nobs = NULL, ...)

## S3 method for class 'unmarkedFitGMM'
anovaOD(mod.simple, mod.complex, c.hat = 1,  
        nobs = NULL, ...)

## S3 method for class 'unmarkedFitGPC'
anovaOD(mod.simple, mod.complex, c.hat = 1,  
        nobs = NULL, ...)

## S3 method for class 'unmarkedFitOccuMS'
anovaOD(mod.simple, mod.complex, c.hat = 1,  
        nobs = NULL, ...)

## S3 method for class 'unmarkedFitOccuTTD'
anovaOD(mod.simple, mod.complex, c.hat = 1,  
        nobs = NULL, ...)

## S3 method for class 'unmarkedFitMMO'
anovaOD(mod.simple, mod.complex, c.hat = 1,  
        nobs = NULL, ...)

## S3 method for class 'unmarkedFitDSO'
anovaOD(mod.simple, mod.complex, c.hat = 1,  
        nobs = NULL, ...)

## S3 method for class 'glmerMod'
anovaOD(mod.simple, mod.complex, c.hat = 1,  
        nobs = NULL, ...)

## S3 method for class 'maxlikeFit'
anovaOD(mod.simple, mod.complex, c.hat = 1,  
        nobs = NULL, ...)

## S3 method for class 'multinom'
anovaOD(mod.simple, mod.complex, c.hat = 1,  
        nobs = NULL, ...)

## S3 method for class 'vglm'
anovaOD(mod.simple, mod.complex, c.hat = 1,  
        nobs = NULL, ...)

Arguments

`mod.simple`	an object of class `glm`, `glmmTMB`, `maxlikeFit`, `mer`, `merMod`, `multinom`, `vglm`, and various `unmarkedFit` classes containing the output of a model. This model should be a simpler version of `mod.complex` resulting from a deletion of certain terms (i.e., nested model).
`mod.complex`	an object of the same class as `mod.simple`.
`c.hat`	value of overdispersion parameter (i.e., variance inflation factor) such as that obtained from `c_hat`, `mb.gof.test`, or `Nmix.gof.test`. Typically, this value should be computed for the most complex model and applied to simpler models.
`nobs`	the number of observations used in the analysis. If `nobs = NULL`, the total number of rows are used as the sample size to compute the residual degrees of freedom as `nobs - K`, where K is the number of estimated parameters. This is relevant only for mixed models or various models of `unmarkedFit` classes where sample size is not straightforward. In such cases, one might use total number of observations or number of independent clusters (e.g., sites) as the value of `nobs`.
`...`	additional arguments passed to the function.

Details

This function applies a correction for overdispersion on the likelihood-ratio test between a model and its simpler counterpart. The simpler model must be nested within the more complex model, typically as the result of deleting terms. You should supply the c.hat value of the most complex of the two models you are comparing.

When 1 < c.hat < 4, the likelihood-ratio test is computed as:

LR = \frac{-2 * (LL.simple - LL.complex)}{(K.complex - K.simple) * c.hat}

where LL.simple and LL.complex are the log-likelihoods of the simple and complex models, respectively, and where K.complex and K.simple are the number of estimated parameters in each model. The test statistic is approximately distributed as F_{K.complex - K.simple, n - K.complex}, where n is the number of observations (i.e., nobs) used in the analysis (Venables and Ripley 2002).

When nobs = NULL, the number of observations is based on the number of rows of the data frame used in the analysis. For mixed models or various models of unmarkedFit, sample size is less straightforward, and nobs could be based on the total number of observations or on the number of independent clusters (e.g., sites), among other choices.

When c.hat = 1, the likelihood-ratio test simplifies to:

LR = -2 * (LL.simple - LL.complex)

where in this case the test statistic is distributed as a \chi^2_{K.complex - K.simple} (McCullagh and Nelder 1989).

The function supports different model types such as Poisson GLM's and GLMM's, single-season and dynamic occupancy models (MacKenzie et al. 2002, 2003), and various N-mixture models (Royle 2004, Dail and Madsen 2011).

Value

anovaOD returns an object of class anovaOD as a list with the following components:

`form.simple`	a character string of the parameters estimated in `mod.simple`.
`form.complex`	a character string of the parameters estimated in `mod.complex`.
`c.hat`	the `c.hat` estimate used to adjust the likelihood-ratio test.
`devMat`	a matrix storing as columns the number of parameters estimated (`K`), the log-likelihood of each model `logLik`, the difference in estimated parameters between the two models (`Kdiff`), minus twice the difference in log-likelihoods between the models (`-2LL`), the test statistic, and the associated P-value.

Author(s)

Marc J. Mazerolle

References

Dail, D., Madsen, L. (2011) Models for estimating abundance from repeated counts of an open population. Biometrics 67, 577–587.

MacKenzie, D. I., Nichols, J. D., Lachman, G. B., Droege, S., Royle, J. A., Langtimm, C. A. (2002) Estimating site occupancy rates when detection probabilities are less than one. Ecology 83, 2248–2255.

MacKenzie, D. I., Nichols, J. D., Hines, J. E., Knutson, M. G., Franklin, A. B. (2003) Estimating site occupancy, colonization, and local extinction when a species is detected imperfectly. Ecology 84, 2200–2207.

Mazerolle, M. J. (2006) Improving data analysis in herpetology: using Akaike's Information Criterion (AIC) to assess the strength of biological hypotheses. Amphibia-Reptilia 27, 169–180.

McCullagh, P., Nelder, J. A. (1989) Generalized Linear Models. Second edition. Chapman and Hall: New York.

Royle, J. A. (2004) N-mixture models for estimating population size from spatially replicated counts. Biometrics 60, 108–115.

Venables, W. N., Ripley, B. D. (2002) Modern Applied Statistics with S. Second edition. Springer-Verlag: New York.

Examples

##anuran larvae example from Mazerolle (2006)
data(min.trap)
##assign "UPLAND" as the reference level as in Mazerolle (2006)          
min.trap$Type <- relevel(min.trap$Type, ref = "UPLAND") 

##run model
m1 <- glm(Num_anura ~ Type + log.Perimeter + Num_ranatra,
          family = poisson, offset = log(Effort),
          data = min.trap) 
##null model
m0 <- glm(Num_anura ~ 1,
          family = poisson, offset = log(Effort),
          data = min.trap) 

##check c-hat for global model
c_hat(m1) #uses Pearson's chi-square/df

##likelihood ratio test corrected for overdispersion
anovaOD(mod.simple = m0, mod.complex = m1, c.hat = c_hat(m1))
##compare without overdispersion correction
anovaOD(mod.simple = m0, mod.complex = m1)


##example with occupancy model
## Not run: 
##load unmarked package
if(require(unmarked)){
   
   data(bullfrog)
     
   ##detection data
   detections <- bullfrog[, 3:9]

   ##assemble in unmarkedFrameOccu
   bfrog <- unmarkedFrameOccu(y = detections)
     
   ##run model
   fm <- occu(~ 1 ~ Reed.presence, data = bfrog)
   ##null model
   fm0 <- occu(~ 1 ~ 1, data = bfrog)

   ##check GOF
   ##GOF <- mb.gof.test(fm, nsim = 1000)
   ##estimate of c-hat:  1.89

   ##display results after overdispersion adjustment
   anovaOD(fm0, fm, c.hat = 1.89)

   detach(package:unmarked)
}

## End(Not run)

[Package AICcmodavg version 2.3-3 Index]