pordlink: Poisson-Ordinal Link Function

Description

Computes the Poisson-ordinal transformation, including its inverse and the first two derivatives.

Usage

pordlink(theta, cutpoint = NULL,
         inverse = FALSE, deriv = 0, short = TRUE, tag = FALSE)

Arguments

theta

Numeric or character. See below for further details.

cutpoint

The cutpoints should be non-negative integers. If pordlink() is used as the link function in cumulative then one should choose reverse = TRUE, parallel = TRUE.

inverse, deriv, short, tag

Details at Links.

Value

See Yee (2012) for details.

Warning

Prediction may not work on vglm or vgam etc. objects if this link function is used.

Details

The Poisson-ordinal link function (POLF) can be applied to a parameter lying in the unit interval. Its purpose is to link cumulative probabilities associated with an ordinal response coming from an underlying Poisson distribution. If the cutpoint is zero then a complementary log-log link is used.

See Links for general information about VGAM link functions.

References

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

Examples

Run this code

# NOT RUN {
pordlink("p", cutpoint = 2, short = FALSE)
pordlink("p", cutpoint = 2, tag = TRUE)

p <- seq(0.01, 0.99, by = 0.01)
y <- pordlink(p, cutpoint = 2)
y. <- pordlink(p, cutpoint = 2, deriv = 1)
max(abs(pordlink(y, cutpoint = 2, inv = TRUE) - p))  # Should be 0

#\ dontrun{ par(mfrow = c(2, 1), las = 1)
#plot(p, y, type = "l", col = "blue", main = "pordlink()")
#abline(h = 0, v = 0.5, col = "orange", lty = "dashed")
#
#plot(p, y., type = "l", col = "blue",
#     main = "(Reciprocal of) first POLF derivative")
#}


# Rutherford and Geiger data
ruge <- data.frame(yy = rep(0:14,
      times = c(57,203,383,525,532,408,273,139,45,27,10,4,0,1,1)))
with(ruge, length(yy))  # 2608 1/8-minute intervals
cutpoint <- 5
ruge <- transform(ruge, yy01 = ifelse(yy <= cutpoint, 0, 1))
fit <- vglm(yy01 ~ 1, binomialff(link=pordlink(cutpoint=cutpoint)), ruge)
coef(fit, matrix = TRUE)
exp(coef(fit))


# Another example
pdata <- data.frame(x2 = sort(runif(nn <- 1000)))
pdata <- transform(pdata, x3 = runif(nn))
pdata <- transform(pdata, mymu = exp( 3 + 1 * x2 - 2 * x3))
pdata <- transform(pdata, y1 = rpois(nn, lambda = mymu))
cutpoints <- c(-Inf, 10, 20, Inf)
pdata <- transform(pdata, cuty = Cut(y1, breaks = cutpoints))
#\ dontrun{ with(pdata, plot(x2, x3, col = cuty, pch = as.character(cuty))) }
with(pdata, table(cuty) / sum(table(cuty)))
fit <- vglm(cuty ~ x2 + x3, data = pdata, trace = TRUE,
            cumulative(reverse = TRUE,
                       parallel = TRUE,
                       link = pordlink(cutpoint = cutpoints[2:3]),
                       multiple.responses = TRUE))
head(depvar(fit))
head(fitted(fit))
head(predict(fit))
coef(fit)
coef(fit, matrix = TRUE)
constraints(fit)
fit@misc$earg
# }

Run the code above in your browser using DataLab

Last chance! 50% off unlimited learning