cenpoisson               package:VGAM               R Documentation

_C_e_n_s_o_r_e_d _P_o_i_s_s_o_n _F_a_m_i_l_y _F_u_n_c_t_i_o_n

_D_e_s_c_r_i_p_t_i_o_n:

     Family function for a censored Poisson response.

_U_s_a_g_e:

     cenpoisson(link = "loge", earg = list(), imu = NULL)

_A_r_g_u_m_e_n_t_s:

    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.

_D_e_t_a_i_l_s:

     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.

_V_a_l_u_e:

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

_W_a_r_n_i_n_g:

     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.

_N_o_t_e:

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

_A_u_t_h_o_r(_s):

     Thomas W. Yee

_R_e_f_e_r_e_n_c_e_s:

     See 'survival' for background.

_S_e_e _A_l_s_o:

     'Surv', 'poissonff', 'Links'.

_E_x_a_m_p_l_e_s:

     # 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 = head(index, 4)
     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

