cmm_reg: Conway Maxwell Multinomial Regression

In andrewraim/COMMultReg: Conway-Maxwell Multinomial Regression

Conway Maxwell Multinomial Regression

Description

Functions for CMM Regression.

Usage

cmm_reg(
  formula_x,
  formula_w = ~1,
  data = NULL,
  beta_init = NULL,
  gamma_init = NULL,
  control = NULL,
  ...
)

cmm_reg_raw(y, X, W, beta_init = NULL, gamma_init = NULL, control = NULL)

## S3 method for class 'cmm_reg'
logLik(object, ...)

## S3 method for class 'cmm_reg'
AIC(object, ..., k = 2)

## S3 method for class 'cmm_reg'
coef(object, extended = FALSE, ...)

## S3 method for class 'cmm_reg'
vcov(object, extended = FALSE, ...)

## S3 method for class 'cmm_reg'
print(x, ...)

Arguments

formula_x	Regression formula to determine \bm{X} design matrix. See details. The outcome must be specified here.
formula_w	Regression formula to determine \bm{W} design matrix.
data	An optional data.frame with variables to be used with regression formulas. Variables not found here are read from the environment.
beta_init	A p_x \times (k-1) matrix whose columns correspond to \bm{β_1}, …, \bm{β_{k-1}}. See details. If provided, will be used for the starting value in Newton-Raphson. If NULL, a default value will be used.
gamma_init	A vector which serves the starting value for γ, if provided. Otherwise, if NULL, a default value will be used.
control	a list of control parameters; see control.
...	Additional optional arguments
y	An n \times k matrix of outcomes, where the ith row \bm{x}_i^\top represents the ith observation.
X	The \bm{X} design matrix; see details.
W	The \bm{W} design matrix; see details.
object	An object output from cmm_reg or cmm_reg_raw.
k	the penalty per parameter to use in AIC; default is 2.
extended	boolean; if FALSE, drop terms associated with the extended parameters μ. See details.
x	Used with the print function; same meaning as object argument.

Details

Fit the model

\bm{Y}_i \sim \textrm{CMM}_k(m_i, \bm{p}_i, ν_i), \quad i = 1, …, n

assuming multinomial logit link

\log≤ft(\frac{p_{i2}}{p_{i1}} \right) = \bm{x}_i^\top \bm{β}_1, \; …, \; \log≤ft(\frac{p_{i,k} }{ p_{i1}} \right) = \bm{x}_i^\top \bm{β}_{k-1},

and identity link

ν_i = \bm{w}_i^\top \bm{γ}.

The first category is assumed to be the baseline by default, but this can be changed to category b by specifying the base = b argument.

The function cmm_reg provides a formula interface, while cmm_reg_raw provides a "raw" interface where variables y, X, and W can be provided directly.

Fitting is carried out with a Newton-Raphson algorithm on the extended parameter \bm{\vartheta} = (\bm{μ}, \bm{ψ}), where \bm{ψ} = (\bm{γ}, \bm{β}_1, …, \bm{β}_{k-1}) are the regression coefficients of interest and \bm{μ} contains elements of the form -\log C(\bm{p}, ν; m) + m \log p_b which are typically not of direct interest to the analyst. See Morris, Raim, and Sellers (2020+) for further details.

Let \bm{\vartheta}^{(g)} denote the estimate from the gth iteration. The algorithm is considered to have converged when ∑_{j} |\vartheta_j^{(g)} - \vartheta_j^{(g-1)}| < ε is sufficiently small, or failed to converge when a maximum number of iterations has been reached. These values can be specified via the control argument.

Several auxiliary functions are provided for convenience:

• cmm_reg_control provides a convenient way to construct the control argument.

• logLik returns the log-likelihood at the solution \hat{\bm{ψ}}.

• AIC returns the Akaike information criterion at the solution.

• coef returns a list with elements beta and gamma at the solution \hat{\bm{ψ}}. If extended = TRUE, an element with μ is also returned.

• vcov returns an estimate of the covariance matrix of \hat{\bm{ψ}} based on the information matrix. If extended = TRUE, elements corresponding to μ are also included.

• print displays estimates and standard errors.

Value

An object of class cmm_reg containing the result.

Examples

# Generate data from CMM regression model
set.seed(1234)

n = 200
m = rep(10, n)
k = 3

x = rnorm(n)
X = model.matrix(~ x)
beta = matrix(NA, 2, k-1)
beta[1,] = -1
beta[2,] = 1
P = t(apply(X %*% beta, 1, inv_mlogit))

w = rnorm(n)
W = model.matrix(~ w)
gamma = c(1, -0.1)
nu = X %*% gamma

y = matrix(NA, n, k)
for (i in 1:n) {
    y[i,] = r_cmm(1, m[i], P[i,], nu[i], burn = 200)
}

dat = data.frame(y1 = y[,1], y2 = y[,2], y3 = y[,3], x = x, w = w)

# Fit CMM regression with formula interface
cmm_out = cmm_reg(formula_x = y ~ x, formula_w = ~w, data = dat)
print(cmm_out)
logLik(cmm_out)
AIC(cmm_out)
coef(cmm_out)
vcov(cmm_out)

# Alteratively, use the "raw" interface
cmm_raw_out = cmm_reg_raw(y, X, W)
print(cmm_raw_out)

andrewraim/COMMultReg documentation built on April 2, 2022, 11:04 p.m.