# How to derive probability distributions for event count time series using R?

I need direction for a robust approach in R for deriving probability distributions for selected time points in event count time series.

In the below illustration:

1. "Period_1" shows elapsed number of months,
2. "1stStateX" shows the number of times elements in the population reach a state of X (it's binary, either an element reaches X or not and it's a "dead state" in that an element can only hit X once),
3. "cumStateX" runs a running cumulative total of 1stStateX,
4. "totalUnits" is the number of units in the population,
5. "rateCumX" is the mean (or cumStateX/totalUnits), and then
6. variance, standard deviation, and standard error; I show my calculation for these at the bottom.

For example, I highlight below period 18 where the cumulative mean rate of reaching X is 47.89%. I would like to derive a probability distribution around that 47.89% in period 18. Perhaps a Poisson Autoregressive Model or Autoregressive Conditional Poisson Model is the way to go for forecasting (I've been fiddling around with R packages ACP and tscount for forecasting), but at this point I'm more interested in simply drawing probability distributions around selected periods.

Could someone please advise on how to proceed for deriving point-in-time probability distributions, R packages that may help, etc.?

In the above, I calculate as follows using data.table (calculations are obvious even if you don't know data.table):

[rateCumX:= cumStateX/totalUnits]
[, var:=totalUnits * rateCumX * (1-rateCumX)]
[, sdv:=sqrt(var)]
[, serr:=sdv/sqrt(totalUnits)]


Revised standard error of the mean for binomial distribution ("serr")(1st 18 rows only of same data from original post), using data.table R package:

,rateCumX:= cumStateX/totalUnits][
,var:=totalUnits * rateCumX * (1-rateCumX)][
,sdv:=sqrt(var)][
,serr:=sdv/totalUnits]


• 1. I think you mean serr := sdv/totalUnits. 2. It's not clear to me what you're asking. Do you want an extension of the binomial distribution approach you're working on, or is that just context and you want an entirely new approach? 3. It seems like a survival analysis approach is more applicable than a time-series approach. Are you particularly keen on a time-series approach, and could you explain why? Nov 26, 2022 at 14:39
• Hi, thank you for the feedback. I'm not keen on time-series, I felt it may be applicable because counts reaching state X is time dependent and the counts demonstrate serial dependence (?), autoregressive models predict future values based on past values which is often the case here, etc. I'd like a recommendation for approach and will look at survival analysis. The elements move among various states across time "X" which is the terminal dead state; an element's movements across the non-X states are predictive of the probability of hitting X, there's correlation among elements. Nov 26, 2022 at 19:28
• So I'm open to any other suggestions! I'll start researching survival analysis with R packages now. I have the data showing each element's movement across states over time. I have transition matrices for these movements. Eventually I'm going to have someone implement MCMC simulation for this, but want to get my head around the best approach for simply looking at a probability distribution of outcomes on an aggregate basis, like I show in the illustration, for specific time periods, given the mean/sdev/serr/etc. what would the normal or lognormal distribution look like etc. Nov 26, 2022 at 19:36
• I think the new serr is more useful. Now you can calculate rateCumX +/- 1.96 * serr and get an approximate 95% confidence interval for rateCumX, or some other confidence interval. That's not necessarily the ultimate result you want, but it's a useful statistic to have. Nov 29, 2022 at 0:13

## 1 Answer

From your description it sounds like a survival analysis approach would be a reasonable way to attack this problem. In particular, because you have a known number of individuals and each individual can be in exclusively one state at a time, it sounds like a multi-state survival model would work. If you're interested in reading about that I'd recommend the competing risk vignette for the R survival package. https://cran.r-project.org/web/packages/survival/vignettes/compete.pdf

For the question at hand, I'll interpret "robust" as meaning a non-parametric or assumption-free method. The Kaplan-Meier estimate is a natural generalization of what you've done so far. Here's an example of using a Kaplan-Meier estimate of the survival curve to generate confidence intervals for the proportion of individuals who will be in the "death" state at time t. From that it's easy to find confidence intervals for the number of individuals who will be in the "death" state.

library(survival)
n <- 1280
event_time <- c(
rep(8, 10), rep(9, 112), rep(10, 114),
rep(11, 69), rep(12, 59), rep(13, 77),
rep(14, 45), rep(15, 42), rep(16, 32),
rep(17, 26), rep(18, 27), rep(19, 21),
rep(20, 19), rep(21, 22), rep(22, 24),
rep(23, 8), rep(24, 10), rep(25, 16),
rep(26, 17), rep(27, 12), rep(28, 5),
rep(29, 5), rep(30, 4), rep(31, 4),
rep(32, 4), rep(33, 4), rep(34, 4),
rep(35, 8), rep(36, 4), rep(37, 3),
rep(38, 6), rep(39, 5), rep(40, 0),
rep(41, n - 818))
event_ind <- c(rep(1, 818), rep(0, n - 818))
mod <- survfit(Surv(event_time, event_ind) ~ 1)
plot(mod, ylab = "surviving proportion", xlab = "time")


Created on 2022-11-26 with reprex v2.0.2

Because the model predicts survival proportions, we need to take 1 minus these values to find the proportions in the death state.

indx <- which(mod$time == 18) # "Death" Probabilities c(lower = 1 - mod$upper[indx],
est = 1 - mod$surv[indx], upper = 1 - mod$lower[indx])
#>     lower       est     upper
#> 0.4508079 0.4789062 0.5055670
# Mean number of units in state 2
c(lower = (1 - mod$upper[indx])*n, est = (1 - mod$surv[indx])*n,
upper = (1 - mod\$lower[indx])*n)
#>    lower      est    upper
#> 577.0341 613.0000 647.1257


Note that this confidence interval is not symmetric around the estimate. We can also compare with the calculation that I think you were trying to do, and see the result is quite comparable.

# Binomial Variance calculation
p_est <- 613/n
var <- p_est*(1-p_est)*n
sdv <- sqrt(var)
sderr <- sdv/n
c(p_est - 1.96*sderr, p_est, p_est + 1.96*sderr)
#> 0.4515388 0.4789062 0.5062737
c((p_est - 1.96*sderr)*n, p_est*n,
(p_est + 1.96*sderr)*n)
#> 577.9697 613.0000 648.0303


This isn't a probability distribution for the point-in-time estimates. The K-M approach is a non-parametric approach, so it doesn't need to assume a distribution. That's a strength and a weakness. If you have a distribution you think would be good fit, a model using that distribution will be able to make more accurate predictions and forecast into the future. The K-M estimate necessarily stops at the end of the data and can't forecast beyond the observed data.

• Have you given me a lot to think about and study. I'm going through the materials you referenced and your example now. THANK YOU! Nov 27, 2022 at 6:45
• Mr. Thiessen, I revised the serr formula per your recommendation in this post and my prior post. Please let me know if you don't agree with my revision. Nov 28, 2022 at 8:49