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.
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)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.
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.
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.
Link functions;
the suffixes .p, .a and .i refer to the parent,
altered and inflated distributions respectively.
See Links for more choices and information.
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.
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.
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.
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.
See CommonVGAMffArguments for information.
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.
Initial values for the rate parameters;
see CommonVGAMffArguments for information.
See CommonVGAMffArguments for information.
Details are at Gaitpois.
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.
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_s\)).
The choice "pstr.mlm" returns
a matrix whose columns are
the inflated probabilities (Greek symbol \(\phi_s\)).
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".
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.
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.
Yee, T. W. and Ma, C. (2021). Generally--altered, --inflated and --truncated regression, with application to heaped and seeped counts. In preparation.
Gaitpois,
multinomial,
rootogram4,
specials,
plotdgait,
spikeplot,
meangait,
poissonff,
gaitlog,
gaitzeta,
poissonff,
zapoisson,
zipoisson,
pospoisson,
CommonVGAMffArguments,
simulate.vlm.
# NOT RUN {
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))
# }
# NOT RUN {
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
# }
Run the code above in your browser using DataLab