Ordinal Poisson Family Function
Fits a Poisson regression where the response is ordinal (the Poisson counts are grouped between known cutpoints).
ordpoisson(cutpoints, countdata = FALSE, NOS = NULL,
Levels = NULL, init.mu = NULL, parallel = FALSE,
zero = NULL, link = "loglink")cutpoints |
Numeric. The cutpoints, K_l.
These must be non-negative integers.
|
countdata |
Logical. Is the response (LHS of formula) in count-data format?
If not then the response is a matrix or vector with values |
NOS |
Integer. The number of species, or more generally, the number of
response random variates.
This argument must be specified when |
Levels |
Integer vector, recycled to length |
init.mu |
Numeric. Initial values for the means of the Poisson regressions.
Recycled to length |
parallel, zero, link |
See |
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,…,L where L (L >= 2) is the number of levels. In more detail, let Y^*=l if K_{l-1} < Y <= K_{l} where the K_l are the cutpoints. We have K_0=-Inf and K_L=Inf. 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,…,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.
An object of class "vglmff" (see vglmff-class).
The object is used by modelling functions such as vglm
and vgam.
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.
Sometimes there are no observations between two cutpoints. If so,
the arguments Levels and NOS need to be specified too.
See below for an example.
Thomas W. Yee
Yee, T. W. (2020). Ordinal ordination with normalizing link functions for count data, (in preparation).
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))
## End(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 zerosPlease choose more modern alternatives, such as Google Chrome or Mozilla Firefox.