- Operating System: Mac
- brms Version: 2.10.3
I’ve inherited some count data of self-reported behaviors. The response format was limited in that the lowest count one could endorse was 1 and the highest count one could endorse was 6. The way I see it, this implies each 1 is a mixture of 1s and 0s and each 6 is a mixture of 6s and above. I believe these could be considered censored count data and I’d like to model them with the Poisson likelihood.
Here’s a simplified version of the data. The
y values are the true counts. The
y_cens values are the counts after censoring as described, above.
# load packages library(tidyverse) library(brms) # simulate data set.seed(1) s <- tibble(y = rpois(n = 1e3, lambda = 3)) %>% mutate(y_cens = ifelse(y < 1, 1, ifelse(y > 6, 6, y)))
cens() argument in
brms::brm() requires an additional column in the data indicating the censoring status of each observation. Here’s my attempt.
s <- mutate(censored = ifelse(y == 1, "left", ifelse(y > 6, "right", "none")))
Here’s my model setup.
fit <- brm(data = s, family = poisson, y_cens | cens(censored) ~ 1)
The model output seems to looks good.
Family: poisson Links: mu = log Formula: y_cens | cens(censored) ~ 1 Data: s (Number of observations: 1000) Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1; total post-warmup samples = 4000 Population-Level Effects: Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS Intercept 1.11 0.02 1.07 1.14 1.00 1649 2219 Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence, Rhat = 1).
Does it appear I’ve correctly interpreted how to work with censored count data?