cenpoisson {VGAM} | R Documentation |
Family function for a censored Poisson response.
cenpoisson(link = "loge", earg = list(), imu = NULL)
link |
Link function applied to the mean.
See Links for more choices.
|
earg |
Extra argument optionally used by the link function.
See Links for more information.
|
imu |
Optional initial value.
See CommonVGAMffArguments for more information.
|
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 Surv
.
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 Surv
,
nor is type="counting"
.
Thomas W. Yee
See survival for background.
# Example 1: right censored data set.seed(123) y = rpois(n <- 100, exp(3)) U = 20 cy = pmin(U, y) rcensored = y >= U table(cy) table(rcensored) status = ifelse(rcensored, 0, 1) table(i <- print(Surv(cy, status))) # Check; U+ means >= U fit = vglm(Surv(cy, status) ~ 1, cenpoisson, trace=TRUE) coef(fit, mat=TRUE) table(print(fit@y)) # Another check; U+ means >= U # Example 2: left censored data L = 15 cy = pmax(L, y) lcensored = y < L # Note y < L, not cy == L or y <= L table(cy) table(lcensored) status = ifelse(lcensored, 0, 1) table(i <- print(Surv(cy, status, type="left"))) # Check fit = vglm(Surv(cy, status, type="left") ~ 1, cenpoisson, trace=TRUE) coef(fit, mat=TRUE) # Example 3: interval censored data Lvec = rep(L, len=n) Uvec = rep(U, len=n) icensored = Lvec <= y & y < Uvec # Neither lcensored or rcensored table(icensored) status = rep(3, n) # 3 means interval censored status = ifelse(rcensored, 0, status) # 0 means right censored status = ifelse(lcensored, 2, status) # 2 means left censored # Have to adjust Lvec and Uvec because of the (start, end] format: Lvec[icensored] = Lvec[icensored] - 1 Uvec[icensored] = Uvec[icensored] - 1 Lvec[lcensored] = Lvec[lcensored] # Remains unchanged Lvec[rcensored] = Uvec[rcensored] # Remains unchanged table(i <- print(Surv(Lvec, Uvec, status, type="interval"))) # Check fit = vglm(Surv(Lvec, Uvec, status, type="interval") ~ 1, cenpoisson, trace=TRUE) coef(fit, mat=TRUE) table(print(fit@y)) # Another check # Example 4: Add in some uncensored observations index = (1:n)[icensored] index = index[1:min(4,length(index))] status[index] = 1 # actual or uncensored value Lvec[index] = y[index] table(i <- print(Surv(Lvec, Uvec, status, type="interval"))) # Check fit = vglm(Surv(Lvec, Uvec, status, type="interval") ~ 1, cenpoisson, trace=TRUE, crit="c") coef(fit, mat=TRUE) table(print(fit@y)) # Another check