touchard (version 2.0.1)

toureg: Touchard 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.

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.

%\item{offset}{a list of offset vectors used (if any),}

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.

Details

Touchard regression with either \(log(\mu)\) or \(log(\lambda)\) 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 \(\delta\) 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 \(\delta\) given current estimate of \(\beta\). Details are given Andrade et al (submitted).

Finally, estimation can be performed by Poisson Quasi-MLE (or Poisson pseudo-MLE): the estimator for is \(\beta\) 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.

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

# NOT RUN {
### 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")


# }