Approximate Bayesian computation in controlled branching processes: the role of summary statistics

Abstract

A controlled branching process is a stochastic growth population model in which the number of individuals with reproductive capacity in each generation is given by a random control function. The purpose of the present work was to examine the Approximate Bayesian Computation sequential Monte Carlo method, and to propose appropriate summary statistics for these processes. The method’s success is shown to rely on this latter issue, and its accuracy is illustrated and compared with a “likelihood free” Markov chain Monte Carlo technique by means of a simulated example. How to extend the method to a controlled multitype branching process is also illustrated, and an application is made to model real data from the cell kinetics field. Both illustrations are developed using the R statistical software environment.

Appendices

Supplementary material

Simulated data

The data for the simulated example presented in Sect. 4 are provided in Table 3. Recall that for the simulated CBP, which starts with $Z_0=1$ individual, the reproduction law is a geometric distribution with parameter $q=0.4$, and for each $k\in \mathbb {N}_0$, the probability distribution of the control variable $\phi _n(k)$ is a binomial distribution with parameters $\xi (k)$ and $\gamma =0.75$, with $\xi (k)= k+\lfloor \log (k)\rfloor$, for each $k\in \mathbb {N}$ and $\xi (0)=0$.

Table 3 Simulated data, $\widetilde{{\mathcal {Z}}}_{30}^{obs}$

Summary statistics

In order to highlight the importance of the choice of the summary statistic in the output of ABC algorithms, we also implemented them removing one of the proposed algorithm’s coordinates, i.e., we considered the summary statistics:

$$\begin{aligned} {\mathcal {S}}_1(\widetilde{{\mathcal {Z}}}_n)&=\left( \frac{\sum _{i=1}^n Z_i}{\sum _{i=0}^{n-1} Z_i},\frac{\phi _{n-1}(Z_{n-1})}{Z_{n-1}}\right) ,\\ {\mathcal {S}}_2(\widetilde{{\mathcal {Z}}}_n)&=\left( \sum _{i=1}^n Z_i,\frac{\phi _{n-1}(Z_{n-1})}{Z_{n-1}}\right) ,\\ {\mathcal {S}}_3(\widetilde{{\mathcal {Z}}}_n)&=\left( \sum _{i=1}^n Z_i,\frac{\sum _{i=1}^n Z_i}{\sum _{i=0}^{n-1} Z_i}\right) . \end{aligned}$$

The results for the parameters $\theta$, $\gamma$, and $\tau _m$ are presented in Tables 4, 5, and 6, respectively. One observes from these tables that removing the number of progenitors in the final generation from the observed sample, and consequently from the summary statistic, i.e., ${\mathcal {S}}_3$, provides similar accuracy measures for $\tau _m$ – a consequence of the first convergence in (3). However, the estimates of its factors separately, i.e., m (dependent on $\theta$) and $\gamma$, get worse. There particularly stands out the less accurate estimate of the control parameter, $\gamma$, with a very wide 95% HPD interval. If the total progeny component is removed, i.e., ${\mathcal {S}}_1$, the integrated squared error (ISE) and the Kullback-Leibler divergence between the posterior densities (KL) for three parameters, $\theta$, $\gamma$, and $\tau _m$, increase significantly. For ISE and KL, the posterior density given by the Gibbs sampler algorithm was taken as the referent. For the latter parameter, the 95% HPD intervals do not contain the true value. Finally, removing the second coordinate in ${\mathcal {S}}$, i.e., ${\mathcal {S}}_2$, results in a worsening of the three accuracy measures for the three parameters considered.

Table 4 Summary of the estimates of the posterior density $\pi (\theta |\widetilde{{\mathcal {Z}}}_{30})$ by the different methods. RMSE: the relative mean square error; ISE: the integrated squared error; KL: the Kullback-Leibler divergence between the posterior densities. For ISE and KL, the posterior density given by the Gibbs sampler algorithm was taken as the referent

Table 5 Summary of the estimates of the posterior density $\pi (\gamma |\widetilde{{\mathcal {Z}}}_{30})$ by the different methods. RMSE: the relative mean square error; ISE: the integrated squared error; KL: the Kullback-Leibler divergence between the posterior densities. For ISE and KL, the posterior density given by the Gibbs sampler algorithm was taken as the referent

Table 6 Summary of the estimates of the posterior density $\pi (\tau _m|\widetilde{{\mathcal {Z}}}_{30})$ by the different methods. RMSE: the relative mean square error; ISE: the integrated squared error; KL: the Kullback-Leibler divergence between the posterior densities. For ISE and KL, the posterior density given by the Gibbs sampler algorithm was taken as the referent

Table 7 Summary of the sensitivity analysis regarding the choice of the prior distribution for the posterior density $\pi (\theta |\widetilde{{\mathcal {Z}}}_{30})$

Table 8 Summary of the sensitivity analysis regarding the choice of the prior distribution for the posterior density $\pi (\gamma |\widetilde{{\mathcal {Z}}}_{30})$

Table 9 Summary of the sensitivity analysis regarding the choice of the offspring and control distributions in the SMC ABC algorithm with the summary statistic and post-processing method. A geometric distribution and a Poisson distribution were used for the reproduction law, and a binomial distribution, a Poisson distribution, and a negative binomial distribution for the control laws

Table 10 Summary of the sensitivity analysis regarding the choice of the offspring and control distributions in the SMC ABC algorithm with the summary statistic and post-processing method without knowledge of the function $\xi (\cdot )$. A geometric distribution and a Poisson distribution were used for the reproduction law, and a binomial distribution, a Poisson distribution, and a negative binomial distribution for the control laws

Sensitivity analysis

We performed a sensitivity analysis regarding the choices of the prior distributions. Recall that both the offspring and the control parameters are probabilities, so that beta distributions seemed to be reasonable options as their prior distributions. In this analysis therefore, we considered different values for the shape parameters of beta distributions. The results for the Gibbs sampler algorithm and the SMC ABC algorithm with the local lineal regression adjustment using the metric $\rho _1$ are summarized in Tables 7 and 8 for the posterior densities $\pi (\theta |\widetilde{{\mathcal {Z}}}_{30})$ and $\pi (\gamma |\widetilde{{\mathcal {Z}}}_{30})$, respectively. They indicate that the estimation of the posterior densities is not sensitive to the choices of the prior distributions.

It is worth mentioning again that the ABC method relies on having ease of sampling from the model. In our implementation, this implies knowledge of the parametric families. An interesting issue is to study the performance of the algorithms when one knows that those distributions can be parametrized with a one-dimensional parameter, but the kind of parametric distribution one should use to that end is unknown. For this purpose, we shall here present a summary of the results of the sensitivity analysis regarding the choice of the parametric families for the offspring and control distributions for three different CBPs. Apart from the previous example, we simulated 30 generations of two additional CBPs starting with one individual: the first of them has a geometric distribution with parameter $\theta =0.918$ as the offspring distribution and the control variables $\phi _n(k)$ follow binomial distributions with parameters $\xi (k)$ and $\gamma =0.1$, and, in the second, the offspring distribution is a geometric distribution with parameter $\theta =0.556$ and it has control variables $\phi _n(k)$ following a binomial distribution with parameters $\xi (k)$ and $\gamma =0.9$. We proposed two offspring laws (geometric and Poisson distributions) and three control laws (binomial, Poisson, and negative binomial distributions), and ran the SMC ABC algorithm for the posterior densities $\pi (m|\widetilde{{\mathcal {Z}}}_{30})$, $\pi (\tau |\widetilde{{\mathcal {Z}}}_{30})$, and $\pi (\tau _m|\widetilde{{\mathcal {Z}}}_{30})$ in each of these six situations for the three examples considered. We opted for these three parameters since it is more reasonable to compare the results based on different distributions in terms of the offspring mean, the parameter $\tau$, and the asymptotic mean growth rate of the process, which are the stable parameters of the process.

The results are given in Table 9. They show that, for the different models used for simulating, the method usually identifies relatively well the offspring mean, the parameter $\tau$, and the asymptotic mean growth rate, which is also indicative of the goodness of our summary statistic. Particularly noteworthy is the case of the second simulated model when a Poisson distribution is taken for the control laws. The algorithm gives a good estimate for the posterior density of the parameter $\tau$, but for that of m (and consequently, of $\tau _m$) not as good as when a binomial or negative binomial control law distribution is assumed. This is related to the magnitude of the parameter $\theta =0.918$. Indeed, when we used the geometric distribution in simulating the model we also estimated the posterior density of $\theta$, resulting in a mean value of 0.9079 and variance of $2.0629\cdot 10^{-5}$, which seems to be a good estimate of that posterior density. However, the small bias observed in this estimate was amplified when we estimated m due to the fact that $m=\theta (1-\theta )^{-1}$.

Similar results were obtained when one is unaware of the existence of the function $\xi (\cdot )$. To examine that problem in this example, we considered the three models described above, and repeated the study by applying the same ABC method for the posterior distributions $\pi (m|\widetilde{{\mathcal {Z}}}_{30})$, $\pi (\tau |\widetilde{{\mathcal {Z}}}_{30})$, and $\pi (\tau _m|\widetilde{{\mathcal {Z}}}_{30})$. Again, we used geometric and Poisson offspring distributions and binomial, Poisson, and negative binomial control law distributions. The results are presented in Table 10, where one observes that there is no significant difference with respect to those in Table 9 except for a slight increase in the RMSE in most cases. This similarity is due to the fact $\varepsilon (k)=(k+\lfloor \log (k)\rfloor )\gamma \sim k\gamma$, as $k\rightarrow \infty$.