Generally-Altered and -Truncated Negative Binomial Regression Family Function (GAT–NB–MLM Variant)

Description

Fits a generally-altered and -truncated negative binomial regression (using a multinomial logit model for the altered values).

Usage

gatnbinomial.mlm(alter = NULL, truncate = NULL, zero = "size",
   lmunb = "loglink", lsize = "loglink", type.fitted = c("mean",
   "pobs.a", "Pobs.a", "prob.a", "prob.i", "prob.t", "lhs.prob"),
   imethod = 1, imunb = NULL, isize = exp(1), ishrinkage = 0.95,
   probs.y = 0.35, cutoff.prob = 0.999, eps.trig = 1e-7,
   max.chunk.MB = 30)

Arguments

Details

The generally-truncated (GT) negative binomial distribution is an ordinary negative binomial distribution with the probability of certain values (given by the truncate argument) being zero. Thus the other probabilities are scaled up. The (0-truncated) positive-negative binomial distribution is a special case (posnegbinomial).

The generally-altered (GA) negative binomial distribution is an ordinary negative binomial distribution with the probability of certain values (given by the alter argument) being modelled using a multinomial logit model (see multinomial). The 0-altered negative binomial distribution is a special case (zanegbinomial) and it is called a hurdle model by some people.

This function can fit both the GA and GT models simultaneously, called the GAT-NB-MLM. It might be applied to heaped data. That is, each special value can be altered or truncated but not both. The default settings make this family function equivalent to negbinomial.

This function implements Fisher scoring and currently does not handle multiple responses for GT or the GA model. Compared to what could be ginbinomial.mlm this family function handles deflation and inflation, therefore handles a wider range of data. For further details please see Gaitnbinom.mlm.

Value

An object of class "vglmff" (see vglmff-class). The object is used by modelling functions such as vglm, and vgam. See gaitpoisson.mlm for more information.

Warning

Due to its flexibility, it is easy to misuse this function; the truncate vector should ideally be not very long and have values that can be justified by the application on hand. Likewise, the alter vector should be short too, and each value should have good justification for being included. Adding unnecessary values to these two arguments willy-nilly is a recipe for disaster.

Regarding truncation, under- or over-flow may occur if the data is ill-conditioned. The response is checked to see that no values equal any values of the truncate vector.

Compared to gaitpoisson.mlm this family function is even more difficult to fit because it is more flexible and involves approximating an infinite series when computing the expected information matrix.

Note

See gaitpoisson.mlm for more information about robustness.

This family function effectively combines what could have been ganbinomial.mlm() and gtnbinomial.mlm() together. The former would have effectively included multinomial inside it. It is possible in the near future that zanegbinomial, zanegbinomialff, Zanegbin, posnegbinomial, Posnegbin be moved into VGAMdata. Likewise for zapoisson, zapoissonff, Zapois.

An argument max.support is not implemented because the upper RHS tail of the NB distribution does not admit practical formulas that can be computed. However, instead of using something like max.support = A one could try something like truncate = (A+1):(A+large) for some positive integer large that is large enough, but within reason. For example, instead of max.support = 20 one could try truncate = 21:30 or truncate = 21:40.

Author(s)

T. W. Yee

See Also

Gaitnbinom.mlm, gaitpoisson.mlm, gatnbinomial.mix, multinomial, negbinomial, posnegbinomial, zanegbinomial, CommonVGAMffArguments, rootogram4, simulate.vlm.

Examples

## Not run: 
avec <- c(10, 20, 30)    # Alter these values
tvec <- 0  # Truncate this value
pobs.a <- logitlink(-(2:4), inverse = TRUE)  # Between 0.02 and 0.12
size1 <- exp(1)
set.seed(1)
gdata <- data.frame(x2 = runif(nn <- 1000))
gdata <- transform(gdata, lambda1 = exp(2))
gdata <- transform(gdata,
  y1 = rgaitnbinom.mlm(nn, size1, mu = lambda1, pobs.a = pobs.a,
                       truncate = tvec, byrow = TRUE, alter = avec))
gatnbinomial.mlm(alter = avec)
(ty1 <- with(gdata, table(y1)))
propn1 <- c(ty1) / sum(ty1)
plot(as.numeric(names(ty1)), propn1, las = 1, xlab = "y",
     yaxs = "i", ylim = c(0, max(propn1) * 1.1), main = "Heaped data",
     ylab = "Proportion", lwd = 3, type = "h", col = "blue")
fit1 <- vglm(y1 ~ 1, trace = TRUE, data = gdata, crit = "coef",
             gatnbinomial.mlm(alter = avec, truncate = tvec))
head(fitted(fit1))
head(predict(fit1))
coef(fit1, matrix = TRUE)
summary(fit1) 
## End(Not run)