Become an expert in R — Interactive courses, Cheat Sheets, certificates and more!
Get Started for Free

gaitpoisson

Generally–Altered, –Inflated and –Truncated Poisson Regression


Description

Fits a generally–altered, –inflated and –truncated Poisson regression by MLE. The GAIT combo model having 5 types of special values is implemented. This allows mixtures of Poissons on nested and/or partitioned support as well as a multinomial logit model for altered and inflated values. Truncation may include the upper tail.

Usage

gaitpoisson(alt.mix = NULL, inf.mix = NULL, alt.mlm = NULL,
     inf.mlm = NULL, truncate = NULL, max.support = Inf,
     zero = c("pobs", "pstr"), eq.ap = FALSE, eq.ip = FALSE,
     parallel.ap = FALSE, parallel.ip = FALSE,
     llambda.p = "loglink", llambda.a = llambda.p, llambda.i = llambda.p,
     type.fitted = c("mean", "lambdas", "pobs.mlm", "pstr.mlm",
     "pobs.mix", "pstr.mix", "Pobs.mix", "Pstr.mix", "nonspecial",
     "Numer", "Denom.p", "sum.mlm.i", "sum.mix.i", "ptrunc.p",
     "cdf.max.s"), gpstr.mix = ppoints(7) / 3,
     gpstr.mlm = ppoints(7) / (3 + length(inf.mlm)), imethod = 1,
     imux = 0.5, ilambda.p = NULL, ilambda.a = ilambda.p,
     ilambda.i = ilambda.p, ipobs.mix = NULL, ipstr.mix = NULL,
     ipobs.mlm = NULL, ipstr.mlm = NULL, byrow.ai = FALSE,
     ishrinkage = 0.95, probs.y = 0.35)

Arguments

truncate, max.support

Vector of truncated values, i.e., nonnegative integers. For the first five arguments (for the special values) a NULL stands for an empty set, and the five sets must be mutually disjoint. Argument max.support enables RHS-truncation, i.e., something equivalent to truncate = (U+1):Inf for some upper support point U specified by max.support.

alt.mix, inf.mix

Vector of altered and inflated values corresponding to finite mixture models. These are described as parametric or structured.

The parameter lambda.p is always estimated. If length(alt.mix) is 1 or more then the parameter pobs.mix is estimated. If length(inf.mix) is 1 or more then the parameter pstr.mix is estimated. If length(alt.mix) is 2 or more then the parameter lambda.a is estimated. If length(inf.mix) is 2 or more then the parameter lambda.i is estimated.

If length(alt.mix) == 1 or length(inf.mix) == 1 then lambda.a and lambda.i are unidentifiable and therefore ignored. In such cases it would be equivalent to moving alt.mix into alt.mlm, and similarly, moving inf.mix into inf.mlm.

Due to its flexibility, it is easy to misuse this function and ideally the values of the above arguments should be well justified by the application on hand. Adding inappropriate or unnecessary values to these arguments willy-nilly is a recipe for disaster, especially for inf.mix. Using alt.mix effectively removes a subset of the data from the main analysis, therefore may result in a substantial loss of efficiency. For seeped values, alt.mix and alt.mlm should be used only. Heaped values can be handled by inf.mlm and inf.mix, as well as alt.mix and alt.mlm.

alt.mlm, inf.mlm

Vector of altered and inflated values corresponding to the multinomial logit model (MLM) probabilities of observing those values—see multinomial. These are described as nonparametric or unstructured.

llambda.p, llambda.a, llambda.i

Link functions; the suffixes .p, .a and .i refer to the parent, altered and inflated distributions respectively. See Links for more choices and information.

eq.ap, eq.ip

Single logical each. Constrain the rate parameters to be equal? See CommonVGAMffArguments for information.

For the GIT–Pois–Pois submodel, after plotting the responses, if the distribution of the spikes above the nominal probabilities has roughly the same shape as the ordinary values then setting eq.ip = TRUE would be a good idea so that lambda.i == lambda.p. And if inf.mix is of length 2 or a bit more, then TRUE should definitely be entertained. Likewise, for heaped or seeped data, setting eq.ap = TRUE (so that lambda.p == lambda.p) would be a good idea for the GAT–Pois–Pois if the shape of the altered probabilities is roughly the same as the parent distribution.

parallel.ap, parallel.ip

Single logical each. Constrain the MLM probabilities to be equal? If so then this applies to all length(alt.mlm) pobs.mlm probabilities or all length(inf.mlm) pstr.mlm probabilities. See CommonVGAMffArguments for information. The default means that the probabilities are generally unconstrained and unstructured and will follow the shape of the data. See constraints.

type.fitted

See CommonVGAMffArguments and below for information. The first value is the default, and this is usually the unconditional mean. Choosing an irrelevant value may result in an NA being returned and a warning, e.g., "pstr.mlm" for a GAT–MLM model.

The choice "lambdas" returns a matrix with at least 1 column and up to 3 of them, corresponding to all those estimated. In order, their colnames are "lambda.p", "lambda.a" and "lambda.i". For other distributions such as gaitlog type.fitted = "shapes" is permitted and the colnames are "shape.p", "shape.a" and "shape.i", etc.

Option "Pobs.mix" provides more detail about "pobs.mix" by returning a matrix whose columns correspond to each altered value; the row sums (rowSums) of this matrix is "pobs.mix". Likewise "Pstr.mix" about "pstr.mix".

The choice "cdf.max.s" is the CDF evaluated at max.support using the parent distribution, e.g., ppois(max.support, lambda.p) for gaitpoisson. The value should be 1 if max.support = Inf (the default). The choice "nonspecial" is the probability of a nonspecial value. The choices "Denom.p" and "Numer" are quantities found in the GAIT combo PMF and are for convenience only.

gpstr.mix, gpstr.mlm

See CommonVGAMffArguments for information. Gridsearch values for the two parameters. If failure occurs try a finer grid, especially closer to 0, and/or experiment with imux.

imethod, ipobs.mix, ipstr.mix, ipobs.mlm, ipstr.mlm

See CommonVGAMffArguments for information.

imux

Numeric, general downward multiplier for initial values for the sample proportions (MLEs actually). The value 1 makes no adjustment, and in general it should lie in (0, 1] with a value near 0.5 recommended. If too high then grid.search() tends to fail. If this occurs another course of action is to set gpstr.mix and/or gpstr.mlm to be a finer grid closer to 0, e.g., gpstr.mix = seq(5) / 100.

ilambda.p, ilambda.a, ilambda.i

Initial values for the rate parameters; see CommonVGAMffArguments for information.

probs.y, ishrinkage

See CommonVGAMffArguments for information.

byrow.ai

Details are at Gaitpois.

zero

See CommonVGAMffArguments for information. By default, all the MLM probabilities are modelled as simple as possible (intercept-only) to help avoid numerical problems, especially when there are many covariates. The Poisson means are modelled by the covariates, and the default vector is pruned of any irrelevant values. To model all the MLM probabilities with covariates set zero = NULL.

For the MLM probabilities, to model pobs.mix only with covariates set zero = c('pstr', 'pobs.mlm'). Likewise, to model pstr.mix only with covariates set zero = c('pobs', 'pstr.mlm').

It is noted that, amongst other things, zipoisson and zipoissonff differ with respect to zero, and ditto for zapoisson and zapoissonff.

Details

The full GAIT–Pois–Pois–MLM–Pois-MLM combo model may be fitted with this family function. There are five types of special values and all arguments for these may be used in a single model. Here, the MLM represents the nonparametric while the Pois refers to the Poisson mixtures. The defaults for this function correspond to an ordinary Poisson regression so that poissonff is called instead. A MLM with only one probability to model is equivalent to logistic regression (binomialff and logitlink).

The order of the linear/additive predictors is best explained by an example. Suppose a combo model has length(a.mix) > 1 and length(i.mix) > 1, a.mlm = 3:5 and i.mlm = 6:9, say. Then loglink(lambda.p) is the first. The second is multilogitlink(pobs.mix) followed by loglink(lambda.a) because a.mix is long enough. The fourth is multilogitlink(pstr.mix) followed by loglink(lambda.i) because i.mix is long enough. Next are the probabilities for the alt.mlm values. Lastly are the probabilities for the inf.mlm values. All the probabilities are estimated by one big MLM and effectively the "(Others)" column of left over probabilities is associated with the nonspecial values. The dimension of the vector of linear/additive predictors is M=12.

Two mixture submodels that may be fitted can be abbreviated GAT–Pois–Pois or GIT–Pois–Pois. For the GAT model the distribution being fitted is a (spliced) mixture of two Poissons with differing (partitioned) support. Likewise, for the GIT model the distribution being fitted is a mixture of two Poissons with nested support. The two rate parameters may be constrained to be equal using eq.ap and eq.ip.

A good first step is to apply spikeplot for selecting candidate values for altering and inflating. Deciding between parametrically or nonparametrically can also be determined from examining the spike plot. Misspecified alt.mix/alt.mlm/inf.mix/inf.mlm will result in convergence problems (setting trace = TRUE is a very good idea.) This function currently does not handle multiple responses. Further details are at Gaitpois.

A well-conditioned data–model combination should pose no difficulties for the automatic starting value selection being successful. Failure to obtain initial values from this self-starting family function indicates the degree of inflation may be marginal and/or a misspecified model. If this problem is worth surmounting the arguments to focus on especially are imux, gpstr.mix and gpstr.mlm. See below for the stepping-stone trick.

Apart from the order of the linear/additive predictors, the following are (or should be) equivalent: gaitpoisson() and poissonff(), gaitpoisson(alt.mix = 0) and zapoisson(zero = "pobs0"), gaitpoisson(inf.mix = 0) and zipoisson(zero = "pstr0"), gaitpoisson(truncate = 0) and pospoisson(). Likewise, if alt.mix and inf.mix are assigned a scalar then it effectively moves that scalar to alt.mlm and inf.mlm because there is no lambda.a or lambda.i being estimated. Thus gaitpoisson(alt.mix = 0) and gaitpoisson(alt.mlm = 0) are the effectively same, and ditto for gaitpoisson(inf.mix = 0) and gaitpoisson(inf.mlm = 0).

A nonparametric special case submodel is the GAIT–Pois–MLM–MLM, which is where the ordinary values have a Poisson distribution, and there are altered and inflated values having unstructured probabilities. Thus the distribution being fitted is a mixture of a Poisson and two MLMs with the support of one of the MLMs being equal to the set of altered values and the other for the inflated values. Hence the probability for each inflated value comes from two sources: the parent distribution and a MLM.

Value

An object of class "vglmff" (see vglmff-class). The object is used by modelling functions such as vglm, rrvglm and vgam.

The fitted.values slot of the fitted object, which should be extracted by the generic function fitted, returns the mean mu by default. The choice type.fitted = "pobs.mlm" returns a matrix whose columns are the altered probabilities (Greek symbol omega). The choice "pstr.mlm" returns a matrix whose columns are the inflated probabilities (Greek symbol phi).

The choice "ptrunc.p" returns the probability of having a truncated value with respect to the parent distribution. It includes any truncated values in the upper tail beyond max.support. The probability of a value less than or equal to max.support with respect to the parent distribution is "cdf.max.s".

The choice "sum.mlm.i" adds two terms. This gives the probability of an inflated value, and the formula can be loosely written down as something like "pstr.mlm" + "Numer" * dpois(inf.mlm, lambda.p) / "Denom".

Warning

Amateurs tend to be overzealous fitting zero-inflated models when the fitted mean is low—the warning of ziP should be heeded. For GAIT regression the warning applies here to all inf.mix and inf.mlm values, not just 0.

Default values for this and similar family functions may change in the future, e.g., eq.ap and eq.ip. Important internal changes might occur too, such as the ordering of the linear/additive predictors and the quantities returned as the fitted values.

Using inf.mlm requires more caution than alt.mlm because gross inflation is ideally needed for it to work safely. Ditto for inf.mix versus alt.mix. Data exhibiting deflation or no inflation will produce numerical problems, hence set trace = TRUE to monitor convergence. More than c.10 IRLS iterations should raise suspicion.

Parameter estimates close to the boundary of the parameter space indicate model misspecification and this can be detected by hdeff.

This function is quite memory-hungry with respect to length(c(alt.mix, inf.mix, alt.mlm, inf.mlm)).

It is often a good idea to set eq.ip = TRUE, especially when length(inf.mix) is not much more than 2 or the values of inf.mix are not spread over the range of the response. This way the estimation can borrow strength from both the inflated and non-inflated values. If the inf.mix values form a single small cluster then this can easily create estimation difficulties—the idea is somewhat similar to multicollinearity.

Note

Numerical problems can easily arise because of the flexibility of this distribution and/or the lack of sizeable inflation; it is a good idea to gain experience with simulated data first before applying it to real data. Numerical problems may arise if any of the special values are in remote places of the support, e.g., a value y such that dpois(y, lambda.p) is very close to 0. This is because the ratio of two tiny values can be unstable.

Good initial values may be difficult to obtain using self-starting procedures, especially when there are covariates. If so, then it is advisable to use a trick: fit an intercept-only model first and then use etastart = predict(int.only.model) to fit the model with covariates. This uses the simpler model as a stepping-stone.

The labelling of the linear/additive predictors has been abbreviated to reduce space. For example, multilogitlink(pobs.mix) and multilogitlink(pstr.mix) would be more accurately multilogitlink(cbind(pobs.mix, pstr.mix)) because one grand MLM is fitted. This shortening may result in modifications needed in other parts of VGAM to compensate.

Author(s)

T. W. Yee

References

Yee, T. W. and Ma, C. (2021). Generally–altered, –inflated and –truncated regression, with application to heaped and seeped counts. In preparation.

See Also

Examples

a.mix <- c(5, 10)  # Alter these values parametrically
i.mlm <- c(4, 14)  # Inflate these values
i.mix <- c(3, 15)  # Inflate these values
tvec <- c(6, 7)   # Truncate these values
pobs.mix <- pstr.mix <- pstr.mlm <- 0.1  # So parallel.ip = TRUE, etc.
max.support <- 20; set.seed(1)
gdata <- data.frame(x2 = runif(nn <- 1000))
gdata <- transform(gdata, lambda.p = exp(2 + 0.5 * x2))
gdata <- transform(gdata,
  y1 = rgaitpois(nn, lambda.p, alt.mix = a.mix, inf.mix = i.mix,
                 pobs.mix = pobs.mix, pstr.mix = pstr.mix,
                 inf.mlm = i.mlm, pstr.mlm = pstr.mlm,
                 truncate = tvec, max.support = max.support))
gaitpoisson(alt.mix = a.mix, inf.mix = i.mix, inf.mlm = i.mlm)
with(gdata, table(y1))
## Not run:  spikeplot(with(gdata, y1)) 
gaitpfit <- vglm(y1 ~ x2, crit = "coef", trace = TRUE, data = gdata,
                 gaitpoisson(alt.mix = a.mix, inf.mix = i.mix,
                             parallel.ip = TRUE,
                             inf.mlm = i.mlm, eq.ap = TRUE, eq.ip = TRUE,
                             truncate = tvec, max.support = max.support))
head(fitted(gaitpfit, type.fitted = "Pstr.mix"))
head(predict(gaitpfit))
t(coef(gaitpfit, matrix = TRUE))  # Easier to see with t()
summary(gaitpfit, HDEtest = FALSE)  # summary(gaitpfit) is better

VGAM

Vector Generalized Linear and Additive Models

v1.1-5
GPL-3
Authors
Thomas Yee [aut, cre], Cleve Moler [ctb] (author of several LINPACK routines)
Initial release
2021-01-13

We don't support your browser anymore

Please choose more modern alternatives, such as Google Chrome or Mozilla Firefox.