toureg: Touchard Regression

In touchard: Touchard Model and Regression

Description

Touchard Regression via either maximum likelihood or quasi-likelihood.

Usage

toureg(formula, data, x = FALSE, y = FALSE, start.beta, start.delta, 
        parscale = rep.int(1, length(start.beta)+1),  maxit, abstol = -Inf, 
        reltol = 1e-6, etol = 1e-6, gtol = 1e-4,
        N=100, eps=1e-6, dm = 10, regress = c("mu", "lambda"), 
        method = c("BFGS", "CG", "Nelder-Mead", "glm", "qp1", "qp2"),  ... )

 

Arguments

formula	an object of class "formula" (or one that can be coerced to that class): a symbolic description of the model to be fitted.
data	an optional data frame containing the variables in the model. If not found in data, the variables are taken from environment(formula), typically the environment from which toureg is called.
start.beta, start.delta	starting values for the parameters in the linear predictor. If missing, the default values explained in the reference below are used.
parscale, maxit, abstol, reltol	arguments passed to control in optim. reltol is alos used as a relative tolerance for termination of iteratively weighted least squares (method="glm") if the relative absolute change in the function (f = loglikelihood) falls below reltol*(reltol + |f|).
gtol	iteratively weighted least squares (method="glm") stops if largest component (in magnitude) of gradient is less than this value.
etol	iteratively weighted least squares (method="glm") stops when the last relative step length is sufficiently small, i.e. below etol*(etol + ||b||_2), where b is current state of the minimizer.
dm	non-zero scalar: with method set to "glm" or "qp2" the estimated delta is obtained by solving some nonlinear equation with root-finding started in the interval start.delta +/- dm.
regress	whether regression is based on log("mu" or log("lambda".
method	optimization method for maximization of loglikelihood: (i) Broyden-Fletcher-Goldfarb-Shanno, Conjugate Gradient or Nelder-Mead as implemented in optim; (ii) iteratively weighted least squares (given delta) combined with optimization over delta (given the regression coefficients) as in GLM-type models; or (iii) quasi-Poisson-Touchard (QPT) method with two variants: "qp1" assumes variance = mu-delta and "qp2" assumes the exact Touchard variance (see also dm).
N, eps	arguments passed to tau.
x, y	logical values indicating whether the response vector and model matrix used in the fitting process should be returned as components of the returned value.
...	not used.

Details

Touchard regression with either log(μ) or log(λ) modeled linearly on the predictors as described in Andrade et al (submitted). Estimation can be performed by maximum likelihood via optim with three available methods ('BFGS', 'CG', 'Nelder-Mead') and analytical gradients. Default starting values for the coefficients are obtained from Poisson GLM. Default starting value for δ is obtained by regressing the metameter on the sufficient statistics Y and log(Y+1). Standard errors are obtained from the diagonal of inverse of observed Fisher information as reported at the final iteration.

Estimation may also be performed by combination of iteratively weighted least squares and maximization over δ given current estimate of β. Details are given Andrade et al (submitted).

Finally, estimation can be performed by Poisson Quasi-MLE (or Poisson pseudo-MLE): the estimator for is β is the same as in the Poisson model (which can be thought of as simply a motivation to the first-order condition defining the estimator); the variance is specified independently without restriction of equidispersion. Two specifications are available: (i) a linear specification variance = mu-delta which corresponds to an approximation to the Touchard variance and (ii) the exact Touchard variance, both allowing for under- and over-dispersion. Details are given Andrade et al (submitted).

Extractor functions for fitted model objects (of class "toureg"): print, summary, plot, residuals, predict, cooks.distance, hatvalues and gleverage.

toureg returns an object of class "toureg", a list with components as described below.

Value

call	the original function call.
coefficients	named vector of estimated regression coefficients.
convergence	integer code from otim indicating either successful completion or faulty termination.
data	the data provided in the function call.
delta	named vector (of length one) of estimated delta parameter.
df	residual degrees of freedom in the fitted model.
fitted.values	a vector of fitted values of lambda.
formula	the formula provided in the function call.
lambda	vector of fitted values of lambda.
loglik	log-likelihood of the fitted model or pseudo-log-likelihood in case os method set to "glm", "qp1" or "qp2".
method	method used.
mu	vector of fitted means.
residuals	vector of raw residuals (y - mu).
se	standard errors of estimated parameters.
start.beta, start.delta	the starting values for the parameters passed to the optimizations routines.
w	weights in the (projection) hat matrix analogous to GLMs.

terms	the 'terms' object used.
var	vector of fitted variances.
vcov	covariance matrix of estimates.
x	if requested, the model matrix.
y	if requested, the response vector.

Author(s)

Bernardo Andrade and Sandro Oliveira

References

Andrade, BB; Matsushita, RY; Oliveira, SB (submitted) Analyzing Count Data with the Touchard Model in R. available upon request.

See Also

glm, formula

Examples

### Horseshoe crab data used by several textbook sources	  
data(Crabs)

### Model Fitting (with different methods) and Plotting
summary( fm <- toureg(y ~ weight + color, data=Crabs) )
# same as
# summary( fm <- toureg(y ~ weight + color, data=Crabs, regress='lambda', method='BFGS') )

# other methods based on log(mu):
# summary( fm2 <- toureg(y ~ weight + color, data=Crabs, regress='mu', method='glm') )
# summary( fm3 <- toureg(y ~ weight + color, data=Crabs, regress='mu', method='qp1') )


plot(fm)
plot(fm , which = 1)
rgram(fm)

### Diagnostics
plot(hvalues(fm))
plot(gleverage(fm))
plot(cooks.dist(fm))


sum(residuals(fm,'response')^2)
sum(residuals(fm,'pearson')^2)
sum(residuals(fm,'deviance')^2)



### Predicted values for 'newdata' ###
 
# Predicted mean values (on the scale of the response variable, i.e. \hat{\mu}):
predict(fm, newdata=data.frame(weight=c(5,6), color=c(2,4)), type="response", se.fit=TRUE)
# Predicted values of lambda:
predict(fm, newdata=data.frame(weight=c(5,6), color=c(2,4)), type="lambda", se.fit=TRUE)
# Predicted values of the linear predictor  x'beta, SEs not yet available:
predict(fm, newdata=data.frame(weight=c(5,6), color=c(2,4)), type="linpred")
# Predicted variances, i.e. \hat{\sigma}^2, SEs not yet available:
predict(fm, newdata=data.frame(weight=c(5,6), color=c(2,4)), type="variance")

touchard documentation built on May 31, 2019, 5:04 p.m.