cenpoisson {VGAM} | R Documentation |
Family function for a censored Poisson response.
cenpoisson(link = "loge", earg = list(), imu = NULL)
link, earg |
Link function and its extra argument 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, cenpoisson, cdata, trace = TRUE) coef(fit, matrix = TRUE) table(print(fit@y)) # 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, cenpoisson, 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) # Neither 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)] # Remains unchanged cdata$Lvec[with(cdata, rcensored)] = cdata$Uvec[with(cdata, rcensored)] # Remains unchanged with(cdata, table(ii <- print(SurvS4(Lvec, Uvec, status, type = "interval")))) # Check fit = vglm(SurvS4(Lvec, Uvec, status, type = "interval") ~ 1, cenpoisson, cdata, trace = TRUE) coef(fit, matrix = TRUE) table(print(fit@y)) # 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, cenpoisson, cdata, trace = TRUE, crit = "c") coef(fit, matrix = TRUE) table(print(fit@y)) # Another check