cens.poisson: Censored Poisson Family Function

In VGAM: Vector Generalized Linear and Additive Models

Censored Poisson Family Function

Description

Family function for a censored Poisson response.

Usage

cens.poisson(link = "loglink", imu = NULL,
             biglambda = 10, smallno = 1e-10)

Arguments

link	Link function applied to the mean; see Links for more choices.
imu	Optional initial value; see CommonVGAMffArguments for more information.
biglambda, smallno	Used to help robustify the code when lambda is very large and the ppois value is so close to 0 that the first derivative is computed to be a NA or NaN. When this occurs mills.ratio is called.

Details

Often a table of Poisson counts has an entry J+ meaning \ge J. This family function is similar to poissonff but handles such censored data. The input requires SurvS4. Only a univariate response is allowed. The Newton-Raphson algorithm is used.

Value

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

Warning

As the response is discrete, care is required with Surv (the old class because of setOldClass(c("SurvS4", "Surv")); see setOldClass), especially with "interval" censored data because of the (start, end] format. See the examples below. The examples have y < L as left censored and y >= U (formatted as U+) as right censored observations, therefore L <= y < U is for uncensored and/or interval censored observations. Consequently the input must be tweaked to conform to the (start, end] format.

A bit of attention has been directed to try robustify the code when lambda is very large, however this currently works for left and right censored data only, not interval censored data. Sometime the fix involves an approximation, hence it is a good idea to set trace = TRUE.

Note

The function poissonff should be used when there are no censored observations. Also, NAs are not permitted with SurvS4, nor is type = "counting".

Author(s)

Thomas W. Yee

References

See survival for background.

See Also

SurvS4, poissonff, Links, mills.ratio.

Examples

# Example 1: right censored data
set.seed(123); U <- 20
cdata <- data.frame(y = rpois(N <- 100, exp(3)))
cdata <- transform(cdata, cy = pmin(U, y),
                          rcensored = (y >= U))
cdata <- transform(cdata, status = ifelse(rcensored, 0, 1))
with(cdata, table(cy))
with(cdata, table(rcensored))
with(cdata, table(print(SurvS4(cy, status))))  # Check; U+ means >= U
fit <- vglm(SurvS4(cy, status) ~ 1, cens.poisson, data = cdata,
            trace = TRUE)
coef(fit, matrix = TRUE)
table(print(depvar(fit)))  # Another check; U+ means >= U

# Example 2: left censored data
L <- 15
cdata <- transform(cdata,
               cY = pmax(L, y),
               lcensored = y <  L)  # Note y < L, not cY == L or y <= L
cdata <- transform(cdata, status = ifelse(lcensored, 0, 1))
with(cdata, table(cY))
with(cdata, table(lcensored))
with(cdata, table(print(SurvS4(cY, status, type = "left"))))  # Check
fit <- vglm(SurvS4(cY, status, type = "left") ~ 1, cens.poisson,
            data = cdata, trace = TRUE)
coef(fit, matrix = TRUE)

# Example 3: interval censored data
cdata <- transform(cdata, Lvec = rep(L, len = N),
                          Uvec = rep(U, len = N))
cdata <-
  transform(cdata,
        icensored = Lvec <= y & y < Uvec)  # Not lcensored or rcensored
with(cdata, table(icensored))
cdata <- transform(cdata, status = rep(3, N))  # 3 == interval censored
cdata <- transform(cdata,
         status = ifelse(rcensored, 0, status))  # 0 means right censored
cdata <- transform(cdata,
         status = ifelse(lcensored, 2, status))  # 2 means left  censored
# Have to adjust Lvec and Uvec because of the (start, end] format:
cdata$Lvec[with(cdata,icensored)] <- cdata$Lvec[with(cdata,icensored)]-1
cdata$Uvec[with(cdata,icensored)] <- cdata$Uvec[with(cdata,icensored)]-1
# Unchanged:
cdata$Lvec[with(cdata, lcensored)] <- cdata$Lvec[with(cdata, lcensored)]
cdata$Lvec[with(cdata, rcensored)] <- cdata$Uvec[with(cdata, rcensored)]
with(cdata,  # Check
 table(ii <- print(SurvS4(Lvec, Uvec, status, type = "interval"))))
fit <- vglm(SurvS4(Lvec, Uvec, status, type = "interval") ~ 1,
            cens.poisson, data = cdata, trace = TRUE)
coef(fit, matrix = TRUE)
table(print(depvar(fit)))  # Another check

# Example 4: Add in some uncensored observations
index <- (1:N)[with(cdata, icensored)]
index <- head(index, 4)
cdata$status[index] <- 1  # actual or uncensored value
cdata$Lvec[index] <- cdata$y[index]
with(cdata, table(ii <- print(SurvS4(Lvec, Uvec, status,
                                     type = "interval"))))  # Check
fit <- vglm(SurvS4(Lvec, Uvec, status, type = "interval") ~ 1,
            cens.poisson, data = cdata, trace = TRUE, crit = "c")
coef(fit, matrix = TRUE)
table(print(depvar(fit)))  # Another check

VGAM documentation built on Sept. 18, 2024, 9:09 a.m.