pois4mlogit: Fit multinomial logit models

Description

pois4mlogit is a procedure that fits a multinomial logit model by calling fitmv to fit multivariate poisson(log) surrogate models, according to the following logic. In the (mixed or not) multinomial logit model, the probabilities $p_{ic}$ of the different categories (or types) $c=1,...,C$ for the $i$th multinomial draw $(n_{i1},...,n_{iC})$ are of the form $$p_{ic}=\frac{e^{\eta_{ic}}}{\sum_{c=1}^C e^{\eta_{ic}}}$$ where each $\eta_{ic}$ is a linear predictor. By contrast, the denominator makes $\log(p_{ic})$ non-linear. This type of model can can be fitted as a Poisson non-linear mixed-effect model (e.g., Chen & Kuo, 2001). spaMM does not have general procedures for fitting non-linear mixed-effect models, but its current procedures for GLMMs have been hacked to fit this specific model, and this includes the use of an iteratively updated offset term (see Details for more information). The correctness of the procedure can be verified for binomial data by comparison to a standard binomial GLMM fit (see Examples).

In a pois4mlogit call, the formula for each type of the multinomial response must contain a term offset(.dynoffset) that declares the above-defined offset term, in addition to the terms specifying each type-specific $\eta_{.c}$. The values of this offset for all $i$ are iteratively updated between successive calls to surrogate Poisson multivariate-response GLMMs, each fitted using fitmv.

This procedure is experimental, in the sense that not all post-fit procedures may run or return meaningful results from objects returned by pois4mlogit, and some caveats should be kept in mind when handling the outputs (see Examples). But fits and basic post-fit inference procedures are operational and appear reliable. There is a specific predict method for such fits: this method returns by default the predicted frequencies, summing to 1 for each multinomial draw (see Examples). plot_effects and LRT can also handle pois4mlogit fits. simulate and bootstrap procedures based on it may be expected to work, except possibly for some esoteric options. confint has a method for these fits, that can be used to obtain intervals for fixed-effect coefficients, but so far it appears frustratingly slow. numInfo can be used to compute a numerical information matrix (at least when its argument which is NULL, retaining all random-effect parameters in the matrix dimensions). This computation may be distinctly useful to provide standard errors of parameters (see Details and Examples).

Since fitmv is called, its various specific features can be used, such as random effects correlated across the different response types (mv syntax), or the additional relationships between random effects that can be specified by the aliases argument, or shared fixed-effect coefficients across response types (X2X argument). Further, additional submodels not belonging to the multinomial model can also be included, but most post-fit procedures have not been checked in this case, and are likely to fail. However, pois4mlogit cannot fit a model including several multinomial models.

Users should quickly learn to avoid non-identifiable models, in particular the following two simple cases: (i) models where all submodels have their own intercept; and (ii) models where a fixed-effect term (which may be simply an intercept) has the same value in the different submodels, as occurs if its fitted coefficient is shared among submodels, and the values of its regressor are also identical among submodels.

Usage

pois4mlogit(submodels, data, to.long=FALSE, init=list(),  
            control=list(wdfac=2, p4m="", grad=FALSE), ...,  
            next_inits=c("ranPars","v_h","fixef"),
            types, n_iter = 1000L, tol=c(1e-3,1e-5), 
            initfn=get_inits_from_fit, progress = FALSE)
# S3 method for pois4mlogit
predict(object, newdata=NULL, verbose=NULL, 
                              na.action=na.omit, ...)
reshape2long(data, types)

Value

A list also inheriting from classes pois4mlogit and HLfit.

Arguments

submodels: Passed to fitmv: see its submodels argument. This must contain as many poisson(log) submodels as there are response types in the multinomial model.
data: Data frame; each line contains a multinomial draw and, as usual, the required predictor variables. The counts for the different response types must form different columns of the data (this is convenient as, e.g., exactly the same data format can be used in binomial fits). The data are passed unchanged to fitmv, unless to.long is set to TRUE. An initial .dynoffset can be provided in the data. It will be internally overwritten after the first iteration.
init: list; passed to fitmv. Beyond initiating the first fitmv call, it controls which parameters have explicit initial values in further iterations (though the initial values themselves are then distinct from those of the first iteration).
next_inits: Character vector. Controls initiation of the next fitmv call from the result of the previous one. If "ranPars" is included, the function specified by the initfn argument is applied to this previous fit to provide initial values for (“outer-optimized”) random-effect parameters in the next fitmv call. See Details for possible alternative to the next-init default. If "v_h" is included, the previous fitmv result is used to provide initial values for the random effects. If "fixef" is included, the previous fitmv result is used to provide initial values for the fixed effects.
to.long: Boolean; for optional reformatting of the data (see Details and Examples).
control: list; passed to fitmv, but may also contain extra elements intrinsic to pois4mlogit, as shown in its default value. These elements are currently for development purposes.
object, newdata, verbose, na.action: these arguments are respectively a pois4mlogit fit result, a data frame, and objects of any of the formats accepted for the verbose and na.action arguments of predict.HLfit. The first three are indeed passed to predict.HLfit after a local change of the .dynoffset in the data.
...: Further arguments passed to fitmv or to predict.HLfit.
types: Character vector: type labels, i.e., column names of the counts for the different mutlinomial response categories in data.
n_iter: Integer: maximum number of iteration of iterative algorithm.
tol: Numeric: tolerance thresholds for determining convergence. For development purposes, not documented.

initfn: function: a function that returns initial values. Currently for development purposes.
progress: Boolean or numeric: whether to print information about number of iterations, and (if progress>1L) about each iteration. Negative values suppress convergence warnings.

Details

The fitting procedure combines (1) the iterative adjustment of GLMMs where the $\log(\sum_{c=1}^C e^{\eta_{ic}})$ values are represented by an offset term, modified over iterations; and (2) ad-hoc modifications of the gradient vector and information matrix of the GLMMs (the latter modifications being ignored in a first step of the fitting procedure). In mixed-effect models, the procedure iteratively updating the ad-hoc offset has to be called for each combination of random-effect parameters considered by the overall fitting procedure. At convergence, the dynamic offset should not change over iterations, and the Poisson model predictions (offset included) for each multinomial draw should sum to the sample size of the draw. Two convergence criteria, the Ocrit and the Scrit respectively, assess these properties. Non-convergence may signal a problem with the fitting procedure, *but* may also signal that the model is not identifiable, in which case it should be modified.

To fit fixed-effect models, by default only the first step of the above procedure is applied, involving iterative calls to fitmv, with iterative adjustments of the offset term, as this is sufficient to get correct final resutls.

An alternative fitting procedure that performs iterative adjustment of the offset within a surrogate Poisson fit has been implemented and may be faster in some cases (particularly with not too strict tol). However, preliminary tests show that convergence of iterations is not well controlled, in particular in mixed-effect models with given large variance of random effects (which is not really avoidable in practice). Fixing this convergence issue is not a top priority.

The standard errors for fixed-effect coefficients shown in the summary of the fit are not those for the multinomial model (see Examples), and may be quite misleading. This stems from the fact that these SEs do not take into account the effect of the fixed-effect coefficients on the likelihood through the offset. By contrast, SEs which take this effect into account can be deduced from the numInfo computation (see Examples, and what the numInfo documentation says about SEs). Further, likelihood ratio tests will provide more reliable tests of fixed effects.

Missing data are handled. This includes the detection and correct handling of cases where (i) a multinomial draw is uninformative because only one type has full information (i.e., response value and all required predictor variables); and (ii) a multinomial draw is informative about some but not all types. In that case, the basic form of the model for a multinomial draw still holds for the relative counts of these types.

When to.long=TRUE, the surrogate Poisson model is fitted to data specified in a long form where each multinomial draw of size $s_i$ is described as $s_i$ multinomial draws of size 1 (so that each Poisson submodel is fitted to a response vector of 0s and 1s). In this long format, counts for the different response types still form different columns (here containing only 0 or 1) of the data. This format is not needed (and not memory-efficient) but is available for comparison, as it was used by Chen & Kuo (2001). The utility function reshape2long can also be used to convert the original data frame to the long format.

References

Chen, Z. and Kuo, L. (2001) A note on the estimation of the multinomial logit model with random effects. The American Statistician 55, 89-95. https://www.jstor.org/stable/2685993

Examples

Run this code

#### Fitting a binomial(logit) model by a bivariate poisson(log) surrogate:
## Toy data: Let us say we observe a color polymorphism of irises in 10 populations...
set.seed(123)
ssize <- 10L
shape <- 0.35
toydata <- data.frame(
  yellow=rbinom(ssize, 16, prob=rbeta(ssize,shape,shape)),
  purple=rbinom(ssize, 16, prob=rbeta(ssize,shape,shape)), # (purple ignored below)
  blue=rbinom(ssize, 16, prob=rbeta(ssize,shape,shape)),
  phenotype=rnorm(ssize)
)
## Standard binomial fit (purple flowers are ignored here)
(byB <- fitme(cbind(yellow,blue) ~ phenotype, family = binomial(), 
      data=toydata))

## Surrogate fit
(byP2 <- pois4mlogit(submodels = list(
  list(yellow ~ offset(.dynoffset) + phenotype, family = poisson()),
  list(blue ~ offset(.dynoffset) + 0, family = poisson())),
  data = toydata, types=c("yellow","blue")))
  
# => Note the different standard errors: the surrogate Poisson fit does not
# provide the "correct" standard errors of the binomial fit.
# However, numInfo(), from which standard errors can be deduced,
# returns correct results:
if (spaMM.getOption("example_maxtime")>0.65) {
  (infoP <- numInfo(byP2)) # practically identical to numInfo(byB)
  sqrt(diag(solve(infoP))) # SEs as in byB
}

# The 'over-parametrized model' below gives standard errors
# closer to the ones from the standard binomial fit.


# Add a trivial gaussian submodel just to show that this can be done:   
(just4fun <- pois4mlogit(submodels = list(
  list(yellow ~ offset(.dynoffset) + phenotype, family = poisson()),
  list(blue ~ offset(.dynoffset) + 0, family = poisson()),
  list(phenotype ~1)),
  data = toydata, types=c("yellow","blue")))

#### Over-parametrized model:

(byP2over <- pois4mlogit(submodels = list(
  list(yellow ~ offset(.dynoffset) + phenotype, family = poisson()),
  list(blue ~ offset(.dynoffset) + phenotype, family = poisson())),
  data = toydata, types=c("yellow","blue")))
  
## Coefficients of binomial model recovered as:
  fixef(byP2over)[1:2]-fixef(byP2over)[3:4]
  
## SEs of coefficients of binomial model recovered as:  
  {
    P2B <- rbind(c(1,0,-1,0),c(0,1,0,-1))
    (vcovP2B <- P2B %*% vcov(byP2over) %*% t(P2B))
  } 
  # practically equivalent to 
  vcov(byB)
  # and the original SEs:
  sqrt(diag(vcovP2B))

## Probabilities of each type:
  matrix(predict(byP2),ncol=2,
         dimnames=list(NULL, c("yellow","blue")))           
  
## logLiks
# The likelihood values contained in 'pois4mlogit' objects are now those for 
# the multinomial model rather than for an underlying surrogate Poisson model,
# so examples previous given here interpreting differences between these 
# likelihoods are no longer relevant.

# Have a look at the 'See also' section for useful features 
# not illustrated in the present examples.

Run the code above in your browser using DataLab