R: Fitting Joint Mean and Dispersion Effects Models

jmdem {jmdem}

R Documentation

Fitting Joint Mean and Dispersion Effects Models

Description

jmdem is used to fit joint mean and dispersion effects models, specified by giving a symbolic description of the linear predictors for the mean and dispersion and a description of the error distribution

Usage

jmdem(mformula, dformula, data, mfamily = gaussian, dfamily = Gamma, 
      weights, subset, dev.type = c("deviance", "pearson"), 
      moffset = NULL, doffset = NULL, mustart = NULL, phistart = NULL, 
      betastart = NULL, lambdastart = NULL, hessian = TRUE, na.action, 
      grad.func = TRUE, fit.method = "jmdem.fit", 
      method = c("Nelder-Mead", "BFGS", "CG", "L-BFGS-B", "SANN", "Brent"), 
      df.adj = FALSE, disp.adj = FALSE, full.loglik = FALSE, 
      beta.first = TRUE, prefit = TRUE, mcontrasts = NULL, 
      dcontrasts = NULL, control = list(...), 
      minv.method = c("solve", "chol2inv", "ginv"), ...)

jmdem.fit(x, y, z = NULL, weights, mfamily = gaussian, dfamily = Gamma, 
          mu, phi, beta, lambda, moffset = NULL, doffset = NULL, 
          dev.type = c("deviance", "pearson"), hessian = TRUE, 
          method = c("Nelder-Mead", "BFGS", "CG", "L-BFGS-B", "SANN", "Brent"), 
          disp.adj = FALSE, df.adj = FALSE, full.loglik = FALSE, 
          control = list(), mintercept = TRUE, dintercept = TRUE, 
          grad.func = TRUE, lower = -Inf, upper = Inf, ...)

Arguments

`mformula`	an object of class "`formula`" (or one that can be coerced to that class): a symbolic description of the mean submodel to be fitted. The details of model specification are given under 'Details'.
`dformula`	a symbolic description of the dispersion submodel to be fitted. The details are also given under 'Details'.
`data`	an optional data frame, list or environment (or object coercible by as.data.frame to a 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 `jmdem` is called.
`mfamily`	a description of the error distribution and link function to be used in the mean submodel. This can be a character string naming a family function, a family function or the result of a call to a family function. (See `family` for details of family functions.)
`dfamily`	a description of the error distribution and link function to be used in the dispersion submodel. (Also see `family` for details of family functions.)
`weights`	an optional vector of 'prior weights' to be used in the fitting process. Should be `NULL` or a numeric vector.
`subset`	an optional vector specifying a subset of observations to be used in the fitting process.
`dev.type`	a specification of the type of residuals to be used as the response of the dispersion submodel. The ML estimates of the jmdem are the optima of either the quasi-likelihood function for deviance residuals, or the pseudo-likelihood function for Pearson residuals.
`moffset`	an a priori known component to be included in the linear predictor of the mean submodel during fitting. This should be NULL or a numeric vector of length equal to the number of cases. One or more offset terms can be included in the formula instead or as well, and if more than one is specified their sum is used. See `model.offset`.
`doffset`	an a priori known component to be included in the linear predictor of the dispersion submodel during fitting. See `model.offset`.
`mustart`, `mu`	a vector of starting values of individual means.
`phistart`, `phi`	a vector of starting values of individual dispersion.
`betastart`, `beta`	a vector of starting values for the regression parameters of the mean submodel.
`lambdastart`, `lambda`	a vector of starting values for the regression parameters of the dispersion submodel.
`hessian`	the method used to compute the information matrix. Hessian matrix will be calculated for `"TRUE"`, Fisher matrix for `"FALSE"`.
`na.action`	a function which indicates what should happen when the data contain `NA`s. The default is set by the `na.action` setting of `options`, and is `na.fail` if that is unset. The 'factory-fresh' default is `na.omit`. Another possible value is `NULL`, no action. Value `na.exclude` can be useful.
`grad.func`	the gradient function will be included in the optimisation for the "`BFGS`", "`CG`" and "`L-BFGS-B`" methods for `"TRUE"`. If it is `NULL`, a finite-difference approximation will be used. For the "`SANN`" method it specifies a function to generate a new candidate point. If it is `NULL` a default Gaussian Markov kernel is used.
`fit.method`	the method to be used in fitting the model. The default method "`jmdem.fit`" uses the general-purpose optimisation (`optim`): the alternative "model.frame" returns the model frame and does no fitting. User-supplied fitting functions can be supplied either as a function or a character string naming a function, with a function which takes the same arguments as `jmdem.fit`. If specified as a character string it is looked up from within the stats namespace.
`method`	the method to be used for the optimisation. See `optim` for details.
`df.adj`	an adjustment factor for the degrees of freedom `(n-p)/n`, where `n` is the number of observations and `p` is the number of parameters to be estimated in `jmdem`, will be multiplied to the likelihood function before the optimisation for `"TRUE"`.
`disp.adj`	an adjustment factor for the dispersion weight will be multiplied to the estimated dispersion parameter during the optimisation for `"TRUE"`. For details, please see McCullagh and Nelder (1989, Ch. 10, P. 362).
`full.loglik`	the full likelihood function instead of the quasi- or pseudo-likelihood function will be used for the optimisation for `TRUE`.
`beta.first`	the mean effects will be estimated (assuming constant sample dispersion) at the initial stage for `TRUE`. For `FALSE`, the dispersion effects will be estimated first (assuming constantly zero mean for the whole sample).
`prefit`	a specfication whether `jmdem` uses `glm` to prefit the starting values of the mean and dispersion parameters. For `FALSE`, the initial parameter values of all the regressors are set to zero and the sample mean and sample dispersion will be used as the starting values of the corresponding submodel intercepts instead. If the submodels have no intercept, all parameters will also be set to zero. The sample mean and sample dispersion will then be used as `mustart` and `phistart` in the internal computation (they will not be officially recorded in `mustart` and `phistart` in the output object). Defaule value is `TRUE`.
`mcontrasts`	an optional list for the mean effect constrasts. See the `contrasts.arg` of `model.matrix.default`.
`dcontrasts`	an optional list for the dispersion effect constrasts. See the `contrasts.arg` of `model.matrix.default`.
`control`	a list of parameters for controlling the fitting process. For `jmdem.fit` this is passed to `jmdem.control`.
`minv.method`	the method used to invert matrices during the estimation process. "`solve`" gives the solutions of a system of equations, "`chol2inv`" gives the inverse from Choleski or QR decomposition and "`ginv`" gives the generalized inverse of a matrix. If none of the methods is specified or if they are specified in a vector such as `c("solve", "chol2inv", "ginv")`, the matrices will be inverted by the methods in the sequence as given in the vector until it is found.
`x`, `y`, `z`	`x` is a mean submodel's design matrix of dimension `n * p`, `z` is a dispersion submodel's design matrix of dimension `n * k`, and `y` is a vector of observations of length `n`. If `z` is `NULL`, the dispersion submodel only contains an intercept.
`mintercept`	a specification whether the intercept term is included in the mean submodel.
`dintercept`	a specification whether the intercept term is included in the dispersion submodel.
`lower`, `upper`	bounds on the variables for the "`L-BFGS-B`" optimisation method.
`...`	For `control`: arguments to be used to form the default control argument if it is not supplied directly. For `jmdem` and `jmdem.fit`: further arguments passed to or from other methods.

Details

A typical predictor has the form response ~ terms where response is the (numeric) response vector and terms is a series of terms which specifies a linear predictor for response.

A terms specification of the form first + second indicates all the terms in first together with all the terms in second with any duplicates removed. A specification of the form first:second indicates the set of terms obtained by taking the interactions of all terms in first with all terms in second. The specification first * second indicates the cross of first and second. This is the same as first + second + first:second.

The terms in the formula will be re-ordered so that main effects come first, followed by the interactions, all second-order, all third-order and so on: to avoid this pass a terms object as the formula.

An additional term response ~ terms + eta can be added to dformula if the mean submodel is nested in the dispersion submodel in the form such that

g(E(y_i))=\boldsymbol{x}_i\boldsymbol{\beta}=\eta_i, h(\phi)=\boldsymbol{z}_i\boldsymbol{\lambda}+\eta_i\gamma.

In the contrary, if the dispersion submodel is nested in the mean submodel such that

g(E(y_i))=\boldsymbol{x}_i\boldsymbol{\beta}+\delta_i\kappa, h(\phi_i)=\boldsymbol{z}_i\boldsymbol{\lambda}=\delta_i,

mformula can be specified as response ~ terms + delta.

Non-NULL weights can be used to indicate that different observations have different dispersions (with the values in weights being inversely proportional to the dispersions); or equivalently, when the elements of weights are positive integers w_i, that each response y_i is the mean of w_i unit-weight observations. For a binomial GLM prior weights are used to give the number of trials when the response is the proportion of successes: they would rarely be used for a Poisson GLM.

If more than one of etastart and mustart is specified, the first in the list will be used. It is often advisable to supply starting values for a quasi family, and also for families with unusual links such as gaussian("log").

glm.fit is the workhorse function: it is not normally called directly but can be more efficient where the response vector, design matrix and family have already been calculated.

Value

`coefficients`	a named vector of estimated coefficients of both the mean and dispersion submodel
`beta`	estimated coefficients of the mean submodel
`lambda`	estimated coefficients of the dispersion submodel
`residuals`	the working residuals, that is the residuals in the final iteration of the `optim` fit. Depending on the type of deviance specified by `dev.type`, `residuals` corresponds to `deviance.residuals` or `pearson.residuals`. Since cases with zero weights are omitted, their working residuals are `NA`.
`deviance.residuals`	the deviance residuals resulting from the final iteration of the `optim` fit.
`pearson.residuals`	the pearson residuals resulting from the final iteration of the `optim` fit.
`fitted.values`	the fitted mean values, obtained by transforming the linear predictors by the inverse of the link function.
`dispersion`	the fitted individual dispersion values, obtained by transforming the linear predictors of the dispersion submodel by the corresponding inverse of the link function.
`mean.rank`	the numeric rank of the fitted mean submodel.
`dispersion.rank`	the numeric rank of the fitted dispersion submodel.
`rank`	the total numeric rank of the fitted model. `mean.rank` and `dispersion.rank` are the corresponding ranks of the fitted mean and dispersion submodels.
`mean.family`	the `family` object used for the mean submodel.
`dispersion.family`	the `family` object used for the dispersion submodel.
`mean.linear.predictors`	the linear fit on link scale of the mean submodel.
`dispersion.linear.predictors`	the linear fit on link scale of the dispersion submodel.
`deviance`	the residual sum of squares of the complete fitted model.
`individual.loglik`	individual value of the log-likelihood function given the estimated mean and dispersion.
`aic`	the Akaike Information Criterion, minus twice the maximised log-likelihood plus twice the number of parameters.
`iter`	number of iterations needed for the fit.
`weights`	the working weights, that is the weights in the final iteration of the `optim` fit.
`prior.weights`	the weights initially supplied, a vector of `1`s if none were.
`info.matrix`	the information matrix given the estimated model coefficients. The diagonal elements of its inverse are the standard errors of the model parameters.
`df.residual`	the residual degrees of freedom of the complete fitted model.
`y`	the `y` vector used.
`x`	the mean submodel design matrix.
`z`	the dispersion submodel design matrix.
`log.llh`	the maximised log-likelihood of the entire sample.
`converged`	logical. Was the `optim` algorithm judged to have converged?
`gradient`	logical. Was the gradient function included in the `optim` algorithm?
`deviance.type`	the type of redidual deviance specified, it is either "`deviance`" or "`pearson`".
`information.type`	the type of information matrix specified, it is either "`Hessian`" or "`Fisher`".
`dispersion.adjustment`	logical. Was the dispersion parameter adjusted by an adjustment factor during the optimisation?
`df.adjustment`	logical. Was the likelihood function adjusted by the degrees of freedom adjustment factor?
`optim.method`	the name of the method used in `optim`.
`control`	the value of the control argument used.
`data`	the evaluated dataset specified in the `data` argument.
`mean.model`	the model frame of the mean submodel.
`dispersion.model`	the model frame of the dispersion submodel.
`call`	the matched call.
`mean.formula`	the formula of the mean submodel supplied.
`dispersion.formula`	the formula of the dispersion submodel supplied.
`fit.method`	the name of the fit function used, currently always "`jmdem.fit`".
`mean.offset`	the offset vector used in the mean submodel.
`dispersion.offset`	the offset vector used in the dispersion submodel.
`dispersion.deviance`	the deviance sum of squares of the dispersion submodel.
`dispersion.df.residual`	the residual degrees of freedom of the dispersion submodel.
`null.deviance`	the residual sum of squares of the complete null model.
`df.null`	the residual degrees of freedom for the complete null model.
`dispersion.null.deviance`	the residual sum of squares of the dispersion null submodel.
`dispersion.df.null`	the residual degrees of freedom for the dispersion null submodel.
`beta.null`	the estimated coefficients of the mean null submodel.
`lambda.null`	the estimated coefficients of the dispersion null submodel.
`dispersion.null`	the estimated dispersion of the complete null model.
`residuals.null`	the residuals of the complete null model.
`mustart`	the vector of starting values for individual means used.
`phistart`	the vector of starting values for individual dispersion used.
`betastart`	the vector of starting values for the mean submodel parameters used.
`lambdastart`	the vector of starting values for the dispersion submodel parameters used.
`mean.terms`	the terms object used for the mean submodel.
`dispersion.terms`	the terms object used for the dispersion submodel.
`xlevels`	a record of the levels of the factors used in fitting the mean submodel.
`zlevels`	a record of the levels of the factors used in fitting the dispersion submodel.
`mean.contrasts`	the contrasts used for the mean submodel.
`dispersion.contrasts`	the contrasts used for the dispersion submodel.
`na.action`	information returned by model.frame on the special handling of `NA`s.
`init.mean.fit`	the initial values of the mean submodel coefficients, linear predictors and fitted values.
`init.dispersion.fit`	the initial values of the dispersion submodel coefficients, linear predictors and fitted values.
`matrix.inverse.method`	information returned on the method used for inverting matrices during optimisation.

Author(s)

Karl Wu Ka Yui (karlwuky@suss.edu.sg)

References

Carroll, R.J., Ruppert, D. (1988). Transformation and Weighting in Regression. London: Chapman and Hall.

Cordeiro, M.G., Simas, A.B. (2009). The distribution of pearson residuals in generalized linear models. Comput. Statist. Data Anal., 53, 3397-3411.

McCullagh, P. (1983). Quasi-likelihood functions. Annals of Statistics 11 (1), 59-67.

McCullagh P. and Nelder, J.A. (1989) Generalized Linear Models. London: Chapman and Hall.

Nash, J.C. (1990). Compact Numerical Methods for Computers. Linear Algebra and Function Minimisation. Adam Hilger.

Nelder, J.A., Lee, Y., Bergman, B., Hynen, A., Huele, A.F., Engel, J. (1998). Joint modelling of mean and dispersion. Technometrics, 40 (2), 168-175.

Nelder, J.A., Pregibon, D. (1987). An extended quasi-likelihood function. Biometrika, 74 (2), 221-232.

Nocedal, J., Wright, S.J. (1999). Numerical Optimization. Springer.

Smyth, G.K. (1989). Generalized linear models with varying dispersion. J. R. Statist. Soc. B, 51 (1), 47-60.

Smyth, G.K., Verbyla, A.P. (1996). A conditional likelihood approach to residual maximum linear estimation in generalized linear models. J. R. Statist. Soc. B, 58 (3), 565-572.

Smyth, G.K., Verbyla, A.P. (1999). Adjusted likelihood methods for modelling dispersion in generalized linear models. Environmetrics, 10, 695-709.

Wedderburn, R. (1974). Quasi-likelihood functions, generalized linear models, and the Gauss-Newton method. Biometrika, 61 (3), 439-447.

Wu, K.Y.K., Li, W.K. (2016). On a dispersion model with Pearson residual responses. Comput. Statist. Data Anal., 103, 17-27.

Examples

## Fit poisson counts by unnested mean and dispersion submodels. 
## Use log-links for both submodels. Set dispersion fitting based 
## on deviance residuals. Use conjugate gradient (CG) as 
## optimisation method.
MyData <- simdata.jmdem.sim(mformula = y ~ x, dformula = ~ z, 
                            mfamily = poisson(), 
                            dfamily = Gamma(link = "log"), 
                            beta.true = c(0.5, 4), 
                            lambda.true = c(2.5, 3), n = 100)
                            
fit <- jmdem(mformula = y ~ x, dformula = ~ z, data = MyData, 
             mfamily = poisson, dfamily = Gamma(link = "log"), 
             dev.type = "deviance", method = "CG")
      
## Fit Gaussian responses by nesting dispersion submodel in the mean 
## submodel. Use default link for both submodels. Set dispersion fitting 
## based on pearson residuals. Use quasi-Newton (BFGS) as optimisation 
## method. Adjust degrees of freedom for the likelihood function.
MyData <- simdata.jmdem.sim(mformula = y ~ x + delta, dformula = ~ z, 
                            mfamily = gaussian(), 
                            dfamily = Gamma(link = "log"), 
                            beta.true = c(0.5, 4, 1), 
                            lambda.true = c(2.5, 3), n = 100)

fit <- jmdem(mformula = y ~ x + delta, dformula = ~ z, data = MyData, 
             mfamily = gaussian, dfamily = Gamma, dev.type = "pearson", 
             method = "BFGS", df.adj = TRUE)

[Package jmdem version 1.0.1 Index]