# Censored Poisson Family Function Family function for a censored Poisson response. ```r 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 $>= 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. ## Returns An object of class `"vglmff"` (see `vglmff-class`). The object is used by modelling functions such as `vglm` and `vgam`. ## References See `survival` for background. ## Author(s) Thomas W. Yee ## Note The function `poissonff` should be used when there are no censored observations. Also, `NA`s are not permitted with `SurvS4`, nor is `type = "counting"`. ## 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`. ## See Also `SurvS4`, `poissonff`, `Links`, `mills.ratio`. ## Examples ```r # 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 ```

cens.poisson function