VGAM (version 1.1-2)

ordpoisson: Ordinal Poisson Family Function

Description

Fits a Poisson regression where the response is ordinal (the Poisson counts are grouped between known cutpoints).

Usage

ordpoisson(cutpoints, countdata = FALSE, NOS = NULL,
           Levels = NULL, init.mu = NULL, parallel = FALSE,
           zero = NULL, link = "loglink")

Arguments

cutpoints

Numeric. The cutpoints, \(K_l\). These must be non-negative integers. Inf values may be included. See below for further details.

countdata

Logical. Is the response (LHS of formula) in count-data format? If not then the response is a matrix or vector with values 1, 2, …, L, say, where L is the number of levels. Such input can be generated with cut with argument labels = FALSE. If countdata = TRUE then the response is expected to be in the same format as depvar(fit) where fit is a fitted model with ordpoisson as the VGAM family function. That is, the response is matrix of counts with L columns (if NOS = 1).

NOS

Integer. The number of species, or more generally, the number of response random variates. This argument must be specified when countdata = TRUE. Usually NOS = 1.

Levels

Integer vector, recycled to length NOS if necessary. The number of levels for each response random variate. This argument should agree with cutpoints. This argument must be specified when countdata = TRUE.

init.mu

Numeric. Initial values for the means of the Poisson regressions. Recycled to length NOS if necessary. Use this argument if the default initial values fail (the default is to compute an initial value internally).

parallel, zero, link

See poissonff.

Value

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

Warning

The input requires care as little to no checking is done. If fit is the fitted object, have a look at fit@extra and depvar(fit) to check.

Details

This VGAM family function uses maximum likelihood estimation (Fisher scoring) to fit a Poisson regression to each column of a matrix response. The data, however, is ordinal, and is obtained from known integer cutpoints. Here, \(l=1,\ldots,L\) where \(L\) (\(L \geq 2\)) is the number of levels. In more detail, let \(Y^*=l\) if \(K_{l-1} < Y \leq K_{l}\) where the \(K_l\) are the cutpoints. We have \(K_0=-\infty\) and \(K_L=\infty\). The response for this family function corresponds to \(Y^*\) but we are really interested in the Poisson regression of \(Y\).

If NOS=1 then the argument cutpoints is a vector \((K_1,K_2,\ldots,K_L)\) where the last value (Inf) is optional. If NOS>1 then the vector should have NOS-1 Inf values separating the cutpoints. For example, if there are NOS=3 responses, then something like ordpoisson(cut = c(0, 5, 10, Inf, 20, 30, Inf, 0, 10, 40, Inf)) is valid.

References

Yee, T. W. (2012) Ordinal ordination with normalizing link functions for count data, (in preparation).

See Also

poissonff, polf, ordered.

Examples

# NOT RUN {
set.seed(123)  # Example 1
x2 <- runif(n <- 1000); x3 <- runif(n)
mymu <- exp(3 - 1 * x2 + 2 * x3)
y1 <- rpois(n, lambda = mymu)
cutpts <- c(-Inf, 20, 30, Inf)
fcutpts <- cutpts[is.finite(cutpts)]  # finite cutpoints
ystar <- cut(y1, breaks = cutpts, labels = FALSE)
# }
# NOT RUN {
plot(x2, x3, col = ystar, pch = as.character(ystar))
# }
# NOT RUN {
table(ystar) / sum(table(ystar))
fit <- vglm(ystar ~ x2 + x3, fam = ordpoisson(cutpoi = fcutpts))
head(depvar(fit))  # This can be input if countdata = TRUE
head(fitted(fit))
head(predict(fit))
coef(fit, matrix = TRUE)
fit@extra

# Example 2: multivariate and there are no obsns between some cutpoints
cutpts2 <- c(-Inf, 0, 9, 10, 20, 70, 200, 201, Inf)
fcutpts2 <- cutpts2[is.finite(cutpts2)]  # finite cutpoints
y2 <- rpois(n, lambda = mymu)   # Same model as y1
ystar2 <- cut(y2, breaks = cutpts2, labels = FALSE)
table(ystar2) / sum(table(ystar2))
fit <- vglm(cbind(ystar,ystar2) ~ x2 + x3, fam =
            ordpoisson(cutpoi = c(fcutpts,Inf,fcutpts2,Inf),
                       Levels = c(length(fcutpts)+1,length(fcutpts2)+1),
                       parallel = TRUE), trace = TRUE)
coef(fit, matrix = TRUE)
fit@extra
constraints(fit)
summary(depvar(fit))  # Some columns have all zeros
# }