ptmixed {ptmixed}

Poisson-Tweedie generalized linear mixed model

Description

Estimates the Poisson-Tweedie generalized linear mixed model with random intercept. Likelihood approximation for the model is based on the adaptive Gauss-Hermite quadrature rule.

Usage

ptmixed(fixef.formula, id, offset = NULL, data, npoints = 5, hessian = T,
  trace = T, theta.start = NULL, reltol = 1e-08, maxit = c(10000, 100),
  freq.updates = 200, min.var.init = 0.001)

Arguments

fixef.formula

A formula for the fixed effects part of the model. It should be in the form y ~ x1 + x2

id

A variable to distinguish observations from the same subject.

offset

An offset to be added to the linear predictor. Default is NULL.

data

A data frame containing the variables declared in fixef.formula.

npoints

Number of quadrature points employed in the adaptive quadrature. Default is 5.

hessian

Logical value. If TRUE, the hessian matrix is evaluated at the MLE to derive the observed Fisher information matrix. Default is TRUE.

trace

Logical value. If TRUE, additional information is printed during the optimization. Default is TRUE.

theta.start

Numeric vector comprising initial parameter values for the vector of regression coefficients, the dispersion parameter, the power parameter and the variance of the random intercept (to be specified exactlyin this order!). Default is NULL: initial parameter estimates are computed automatically by the function.

reltol

Relative tolerance to be used in optim. Default to 1e-8

maxit

Vector containing the maximum number of iterations used in optim by the Nelder-Mead method and, if this fails, by the BFGS method

freq.updates

Number of iterations after which the quadrature points are updated when the Nelder-Mead algorithm is used for the optimization. Default value is 200. To update the quadrature points at every iteration (note that this may make the computation about 10x slower), set freq.updates = 1 or freq.updates = NA. The function first tries to optimize the loglikelihood using the Nelder-Mead algorithm, updating the quadrature points every freq.updates iterations. If this fails to converge, a second attempt is made using the BFGS algorithm, for which the quadrature points are updated at every iteration.

min.var.init

If the initial estimate of the variance of the random intercept is smaller than this value, estimation is stopped and the user is advided to use the simpler Poisson-Tweedie GLM is used. Default is 1e-3.

Value

A list containing the following elements: function's call (call); maximum likelihood estimate (mle); value of the loglikelihood at the mle (logl); convergence value (if 0, the optimization converged); the observed Fisher information (fisher.info), if hessian = T; the number of quadrature points used (quad.points) and the starting value used in the optimization (theta.init); relevant warnings (warnings).

Author(s)

Mirko Signorelli

References

Signorelli, M., Spitali, P., Tsonaka, R. (2021). Poisson-Tweedie mixed-effects model: a flexible approach for the analysis of longitudinal RNA-seq data. Statistical Modelling, 21 (6), 520-545. URL: https://doi.org/10.1177/1471082X20936017

See Also

summary.ptglmm, ranef

Examples

data(df1, package = 'ptmixed')
head(df1)

# 1) Quick example (just 1 quadrature point, hessian and SEs 
# not computed - NB: we recommend to increase the number of
# quadrature points to obtain much more accurate estimates,
# as shown in example 2 below where we use 5 quadrature points)

# estimate the model
fit1 = ptmixed(fixef.formula = y ~ group + time, id = id,
              offset = offset, data = df1, npoints = 1, 
              freq.updates = 200, hessian = FALSE, trace = TRUE)
              
# print summary:
summary(fit1, wald = FALSE)


# 2) Full computation that uses more quadrature points
# for the likelihood approximation and includes numeric
# evaluation of the hessian matrix

# estimate the model:
fit2 = ptmixed(fixef.formula = y ~ group + time, id = id,
              offset = offset, data = df1, npoints = 5, 
              freq.updates = 200, hessian = TRUE, trace = TRUE)
              
# print summary:
summary(fit2)

# extract summary:
results = summary(fit2)
ls(results)
results$coefficients