How to fit a negative binomial distribution in R while incorporating censoring

I need to fit $Y_{ij} \sim NegBin(m_{ij},k)$, i.e. a negative binomial distribution to count data. However, the data I have observed are censored - I know the value of $y_{ij}$, but it could be more than that value. The log-likelihood is

\begin{equation} ll = \sum_{i=1}^n w_i (c_i \log(P(Y_{ij}=y_{ij}|X_{ij})) + (1- c_i) \log(1- \sum_{k=1}^32 P(Y_{ij} = k|X_{ij}))) \end{equation}

where $X_{ij}$ represent the design matrix (with the covariates of interest), $w_i$ is the weight for each observation, $y_{ij}$ is the response variable and $P(Y_{ij}=y_{ij}|X_{ij})$ is the negative binomial distribution where the $m_{ij}=exp(X_{ij} \beta)$ and $\alpha$ is the over-dispersion parameter.

Does anyone know of an R package to tackle this problem?

You can try gamlss.cens package.

Another R package that seems to do what you want, is pscal. The associated vignette has lots of examples.