cens.poisson {VGAM} | R Documentation |
Family function for a censored Poisson response.
cens.poisson(link = "loge", imu = NULL)
link |
Link function applied to the mean;
see |
imu |
Optional initial value;
see |
Often a table of Poisson counts has an entry J+ meaning
>= 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.
An object of class "vglmff"
(see vglmff-class
).
The object is used by modelling functions such as
vglm
and
vgam
.
As the response is discrete,
care is required with Surv
, 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.
The function poissonff
should be used
when there are no censored observations.
Also, NA
s are not permitted with SurvS4
,
nor is type = "counting"
.
Thomas W. Yee
See survival for background.
# 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(ii <- 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(ii <- 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 means 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 cdata$Lvec[with(cdata, lcensored)] <- cdata$Lvec[with(cdata, lcensored)] # Unchanged cdata$Lvec[with(cdata, rcensored)] <- cdata$Uvec[with(cdata, rcensored)] # Unchanged 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) 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