# Fast Estimation of the Composite Link Model for Multidimensional Grouped Counts

María Durbán, Carlo Giovanni Camarda

► To cite this version:

María Durbán, Carlo Giovanni Camarda. Fast Estimation of the Composite Link Model for Multi-dimensional Grouped Counts. David O'Sullivan; James Sweeney; Shirin Moghaddam; Kevin Burke. Proceedings of the 39th International Workshop on Statistical Modelling, Mathematics Applications Consortium for Science & Industry (MACSI), pp.128-132, 2025. hal-05419633

HAL Id: hal-05419633

https://hal.science/hal-05419633v1

Submitted on 16 Dec 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
# Fast Estimation of the Composite Link Model for Multidimensional Grouped Counts

María Durbán <fnref n="1" />, Carlo G. Camarda <fnref n="2" />

E-mail for correspondence: marialuz.durban@uc3m.es

**Abstract:** This paper introduces a new approach for estimating the Composite Link Model within a penalized likelihood framework, aimed to address indirect observations of grouped count data. While the model is effective in these contexts, its application becomes computationally challenging in large, high-dimensional settings. To overcome this, we propose a reformulated iterative estimation procedure that leverages Generalized Linear Array Models, enabling the disaggregation and smooth estimation of latent distributions in multidimensional data. Through applications to high-dimensional mortality datasets, we demonstrate the model's ability to capture fine-grained patterns while comparing its computational performance to the conventional algorithm.

**Keywords:** Grouped counts; Composite Link Model; Generalized Linear Array Models

## 1 The Penalized Composite Link Model (PCLM) over multidimensional arrays

Suppose that we observe an array of counts $\mathbf{Y} = (y_{i,j})$ along with explanatory variables $\mathbf{x}_1 = (1, \dots, n_1)$ and $\mathbf{x}_2 = (1, \dots, n_2)$. We assume that $Y_{i,j} \sim P(\mu_{i,j})$. In our setting, while we observe $\mathbf{y}$, a vector of length $n_1 n_2$, we aim to estimate the latent distribution $\gamma = (\gamma_{11}, \gamma_{12}, \dots, \gamma_{m_1 m_2})$, $n_d < m_d$ $d = 1, 2$, which is related to $\mu$ via the composition matrix $\mathbf{C}$, and describes how the elements of the latent vector $\gamma$ are combined to yield $\mu$. The resulting Composite Link Model for Poisson aggregated counts can thus be written as follows:

$$ y \sim P(\mu) \quad \mu = \mathbf{C}\gamma \quad \gamma = \exp(\eta), \quad (1) $$

where $\mathbf{C}$ describes how the elements of the latent vector $\gamma$ are combined to yield $\mu$, and $\eta$ is the linear predictor in a Poisson framework.

This paper was published as a part of the proceedings of the 39th International Workshop on Statistical Modelling (IWSM), Limerick, Ireland, 13–18 July 2025. The copyright remains with the author(s). Permission to reproduce or extract any parts of this abstract should be requested from the author(s).

<footnote n="1">¹ Department of Statistics, Universidad Carlos III de Madrid, Spain</footnote>

<footnote n="2">² Institut National d'Études Démographiques, Aubervillers, France</footnote>
In our first application, we will use Swedish mortality data for females spanning the years 1960 to 2019, covering ages 10 to 104, i.e. $m_1 = 95$ and $m_2 = 60$. In our analysis, the death counts are aggregated into 5-year age groups and across 5-year calendar intervals, i.e. $n_1 = 19$ and $n_2 = 12$. The consistent aggregation structure results in a composition matrix of the form $C = C_2 \otimes C_1$. Our objective is to estimate mortality at a more granular level, specifically for individual ages and single calendar years, to provide finer insights into mortality patterns.

We characterize the latent distribution $\gamma$ by modeling the linear predictor $\eta$. However, estimating a distribution with $m_1 m_2$ values, significantly exceeding the number of observed counts $n_1 n_2$, necessitates the introduction of structural assumptions to constrain the latent space and ensure model identifiability. We suggest that smoothness is a reasonable assumption for $\gamma$, as proposed by Eilers (2007), since smoothness enables the extraction useful information from data that might initially seem insufficient to provide meaningful insights. Given the two-dimensional nature of the problem, we opted for a flexible, fully interactive bivariate smoothing function over $\mathbf{x}_1$ and $\mathbf{x}_2$:

$$ \eta = f(\mathbf{x}_1, \mathbf{x}_2) = \mathbf{B}\alpha = (\mathbf{B}_2 \otimes \mathbf{B}_1)\alpha, \quad (2) $$

where $\mathbf{B}_d, d = 1, 2$ are marginal B-spline basis over $\mathbf{x}_d$, of size $m_d \times c_d, d = 1, 2$. The coefficients corresponding to each basis, denoted as $\alpha$, which need to be estimated, must be penalized in order to derive the smooth latent distributions.

## 2 Model estimation and uncertainty quantification

The penalized log-likelihood for the Poisson PCLM in (1), where the linear predictor is defined as in (2), is given by:

$$ \ell_P = \mathbf{y}' \ln[\mathbf{C} \exp(\mathbf{B}\alpha)] - \mathbf{1}'_{n_1 n_2} \mathbf{C} \exp(\mathbf{B}\alpha) - \frac{1}{2} \alpha' \mathbf{P} \alpha, \quad (3) $$

where $\mathbf{P}$ represents a penalty term that enforces smoothness across both dimensions in the two-dimensional case, and is expressed as:

$$ \mathbf{P} = \lambda_1(\mathbf{I}_{c_2} \otimes \mathbf{D}'_1 \mathbf{D}_1) + \lambda_2(\mathbf{D}'_2 \mathbf{D}_2 \otimes \mathbf{I}_{c_1}). \quad (4) $$

The derivative of (3) with respect to the vector $\alpha$ is given by

$$ \frac{\partial \ell_P}{\partial \alpha} = \mathbf{B}' \mathbf{\Gamma} \mathbf{C}' \mathbf{W}^{-1} \mathbf{y} - \mathbf{B}' \exp(\mathbf{B}\alpha) - \mathbf{P}\alpha, \quad (5) $$

where: $\mathbf{W} = \text{diag}(\mu)$ and $\mathbf{\Gamma} = \text{diag}(\gamma)$. In the original work by Thompson and Baker (1981), the solution to the system in equation (5) was approached by defining a working matrix $\tilde{\mathbf{B}} = \mathbf{W}^{-1} \mathbf{\Gamma} \mathbf{C} \mathbf{B}$ and solving the associated GLM using standard techniques. However, as the dimensionality of the data array increases, the explicit construction of matrix $\tilde{\mathbf{B}}$ becomes computationally prohibitive, causing a significant escalation in memory and processing requirements.

To address this issue, we propose bypassing the direct computation of $\tilde{\mathbf{B}}$ by introducing the concept of a working latent response $\tilde{\mathbf{y}} = \mathbf{\Gamma} \mathbf{C}' \mathbf{W}^{-1} \mathbf{y}$, which represents the redistribution of the observed counts at the desired resolution. In fact, this would correspond to an E-step in the EM algorithm. However, a key distinction between the estimation approach proposed for the PCLM and the EM algorithm
lies in their underlying assumptions about the data distribution. The PCLM approach assumes that the observed aggregated (incomplete) counts, denoted as $\mathbf{y}$, follow a Poisson distribution with mean $\mu$, and estimates an unobservable vector $\gamma$, which is related to $\mu$ via the composition matrix $\mathbf{C}$ and is assumed to be smooth. In contrast, the EM algorithm assumes that the complete (latent) data follow a Poisson distribution with mean $\gamma$ and alternates between the re-distribution of the counts proportionally to the current approximation and the estimation of the distribution parameters.

Defining the working latent response allows us to apply Generalized Linear Array Models (Currie, et al., 2006), eliminating the need to compute $\tilde{\mathbf{B}}$ or any Kronecker products involved in the estimation process. Consequently, our method is not a variant of the EM algorithm. Instead, it represents a computational strategy specifically designed to enable the estimation of a PCLM in a multidimensional framework.

By employing the standard Newton-Raphson method, we obtain the following system of equations:

$$ (\mathbf{B}'\tilde{\Gamma}\mathbf{B} + \mathbf{P})\hat{\alpha} = \mathbf{B}'\tilde{\Gamma}\tilde{\mathbf{z}}, \quad (6) $$

where $\tilde{\mathbf{z}} = \tilde{\eta} + \tilde{\Gamma}^{-1}(\tilde{\mathbf{y}} - \tilde{\gamma})$ serves a role analogous to the conventional working vector found in the GLM context, albeit evaluated at the scale of the latent response.

In order to quantify the uncertainty of the estimates we derive the approximate distribution of $\hat{\alpha}$ using a Bayesian approach which yields:

$$ \alpha|\mathbf{y} \sim N(\hat{\alpha}, (\mathbf{B}'\tilde{\Gamma}\mathbf{C}'\mathbf{W}^{-1}\mathbf{C}\mathbf{\Gamma}\mathbf{B}' + \mathbf{P})^{-1}). \quad (7) $$

In practice, equation (7) is applied by substituting the estimate of $\gamma$ (and therefore $\mu$) evaluated at $\alpha = \hat{\alpha}$ once the algorithm has reached convergence (in a similar fashion as in the GLM case).

There are two main computational advantages of our algorithm:

i) All linear functions (the elements of $\mathbf{B}\hat{\alpha}$, $\tilde{\mu}$ and $\tilde{\mathbf{y}}$), and inner products (elements of $\mathbf{B}'\tilde{\Gamma}\mathbf{B}$) can be re-arrange to fit in the GLAM framework.

ii) The variance covariance matrix in (7) is computed only once, after the algorithm has converged.

In contrast, the original estimation procedure proposed in Eilers (2007) requires $\mathbf{V}$ to be recalculated at every iteration of the algorithm. These two factors greatly reduce the computational time required to fit the model and, in some cases, enable the fitting of data that would otherwise be infeasible.

## 3 Results

By applying a PCLM to the data described previously, we estimate the latent mortality, $\gamma$, of length $m_1m_2 = 5700$.

Figure 1 illustrates the disaggregation capability of the PCLM for four selected ages, showing both the uncertainty associated with each estimated time trend (95% confidence intervals) and the disaggregation capability of the PCLM: while observed log-rates are grouped into five-by-five age-year categories, our estimates provide continuous values for each single year and age along with the uncertainty associated with each estimated time trend.
<figure><img src="image_2.png" /><figcaption>FIGURE 1. Observed death rates by 5-year age groups and 5-year calendar intervals, along with estimated ungrouped mortality by single year of age and calendar year, for two selected age groups.</figcaption></figure>

The proposed methodology produces equivalent outcomes to those obtained using the approach outlined by Eilers (2007). However, by embedding the PCLM within a GLAM framework, our approach provides a substantial advantage, particularly in large-scale applications, where its computational efficiency and scalability make it a superior choice. Fitting the model with the proposed method and described settings takes approximately 0.75 seconds, in contrast to 7.7 seconds required by the original algorithm. Furthermore, the proposed approach is highly efficient in terms of data storage. The R workspace containing all data and outcomes requires only 3.1 MB with the proposed method, whereas the conventional algorithm demands a substantially larger file size of 256 MB.

Although these improvements may appear modest given current computational power, the proposed approach provides a clear solution to potential storage limitations and delivers remarkable speed enhancements, achieving order-of-magnitude improvements over direct evaluation. Furthermore, in practical applications, the search for optimal smoothing parameters often requires repeated execution of the scoring algorithm, where the significant reduction in computational time becomes increasingly impactful. A second application, based on a more complex dataset of weekly Spanish male mortality (not shown here due to space constraints), will further demonstrate the potential of the proposed method in a three-dimensional setting.
# 4 Conclusions

In demography and related fields, data are often aggregated over time or age, obscuring important latent patterns. The Penalized Composite Link Model (PCLM) helps recover fine-scale trends from such grouped data, but traditional estimation methods become infeasible as dimensionality increases. To overcome this, we re-formulate the PCLM estimation by defining a working latent response instead of a regression matrix, similar to the EM algorithm's E-step but without assuming a latent distribution. This avoids alternating between latent and observed data and frames the problem as a Generalized Linear Array Model (GLAM), drastically improving computational efficiency. The method proves effective in large-scale applications and can potentially extend to higher-dimensional or more flexible models beyond Poisson.

**Acknowledgments:** This work is part of the grant PID2022-137243OB-I00 funded by MCIN/AEI/10.13039/501100011033 and, by "ERDF A way of making Europe". This support is gratefully acknowledged.

# References

* Currie, I. D. and Durban, M. and Eilers, P. H. C (2006). Generalized linear array models with applications to multidimensional smoothing. Journal of the Royal Statistical Society: Series B, **68**, 259–280.

* Eilers, P.H.C (2007). Ill-posed problems with counts, the composite link model and penalized likelihood. Statistical Modelling, **7**, 239–254.

* Thompson, R. and Baker, R. J. (1981). Composite link functions in generalized linear models. Journal of the Royal Statistical Society, Series C, **30**, 125–131.