Smoothing Lexis diagrams using kernel functions: A contemporary approach

Abstract

Lexis diagrams are rectangular arrays of event rates indexed by age and period. Analysis of Lexis diagrams is a cornerstone of cancer surveillance research. Typically, population-based descriptive studies analyze multiple Lexis diagrams defined by sex, tumor characteristics, race/ethnicity, geographic region, etc. Inevitably the amount of information per Lexis diminishes with increasing stratification. Several methods have been proposed to smooth observed Lexis diagrams up front to clarify salient patterns and improve summary estimates of averages, gradients, and trends. In this article, we develop a novel bivariate kernel-based smoother that incorporates two key innovations. First, for any given kernel, we calculate its singular values decomposition, and select an optimal truncation point—the number of leading singular vectors to retain—based on the bias-corrected Akaike information criterion. Second, we model-average over a panel of candidate kernels with diverse shapes and bandwidths. The truncated model averaging approach is fast, automatic, has excellent performance, and provides a variance-covariance matrix that takes model selection into account. We present an in-depth case study (invasive estrogen receptor-negative breast cancer incidence among non-Hispanic white women in the United States) and simulate operating characteristics for 20 representative cancers. The truncated model averaging approach consistently outperforms any fixed kernel. Our results support the routine use of the truncated model averaging approach in descriptive studies of cancer.

1. Introduction

The Lexis diagram is a fundamental construct in epidemiology, demography, and sociology. ¹ A Lexis diagram is a rectangular grid of square cells with age along one axis and calendar period along the other. Individuals from a surveilled population contribute person-years and events (births, deaths, cancers, etc.) to each cell. Analysis of Lexis diagrams can elucidate temporal patterns and provide clues about the etiology of an event. ² Notably, the analysis of Lexis diagrams is a cornerstone of cancer surveillance research.

Population-based cancer rates usually exhibit Poisson-type variability. ³ This is the default mode for analysis, although negative binomial models are available that accommodate more flexible mean-to-variance relationships.^4–6 Within any given Lexis diagram, intrinsic variability (“noise”) can mask important signals and limit the power of comparative analyses to identify heterogeneity. In epidemiologic research, the classic approaches deal with intrinsic variability by aggregating granular data (e.g. data originally sampled at the scale of a single-years of age and/or calendar years) into broader age and period categories, ⁷ typically, 5-year age groups. Unfortunately, there is nothing optimal about such traditional groupings.

Several methods have been proposed to de-noise or smooth an observed Lexis diagram before the usual descriptive^7–9 and analytical ¹⁰ methods are applied. Keiding ¹ pioneered the first such approach using bivariate kernel functions. This was a major advance; however, this classic method has two major limitations. First, the choice of kernel is arbitrary, hence so too is the implicit bias-variance trade-off. Second, the estimated variance-covariance matrix does not reflect uncertainty arising from kernel (model) selection. In practice, analysts tend to use small bandwidths to minimize bias rather than variance. However, there is no evidence that such choices are optimal.

Other investigators have developed more sophisticated methods. Currie et al. ¹¹ produced smooth Lexis diagrams using bivariate P-splines and second-difference penalty functions. Camarda ¹² implemented the Currie approach in R and developed methods to incorporate established demographic constraints using asymmetric penalties. ¹³ Dokumentov et al. ¹⁴ developed a hybrid approach that combines smoothing methods with additional parameters that account for abrupt changes in age incidence by period and cohort. Chien et al.^15,16 developed a fully Bayesian smoothing approach using two-dimensional Bernstein polynomials, data-dependent prior distributions, and the Metropolis-Hastings reversible jump algorithm. Martinez-Hernandez and Genton ¹⁷ developed nonparametric methods for the situation where the data can be viewed as a functional trend in one temporal dimension (e.g. age) that is modulated over the other dimension (e.g. time).

In this article, we develop a novel and complementary kernel-based approach using truncated bivariate kernel functions^18–20 and information theory. ²¹ We refer to these algorithms as “filters” and the outputs as “filtrations” because they do more than simply “smooth” the data: by design, they pass through all signals that materially reduce the within-Lexis mean squared error. Our approach has several attractive features. It is fast, and automatic, provides the smoothest possible kernel-based estimates consistent with data, yields a variance–covariance matrix that takes model selection into account, and has excellent performance.

In Section 2, we review notation and background on Lexis diagrams and kernel functions and assemble a representative panel of incident cancers. In Section 3, we describe the new filters and corresponding variance calculations. In Section 4, we illustrate our proposed methods using invasive estrogen receptor (ER)-negative breast cancer incidence among non-Hispanic white women in the United States, and in Section 6, we simulate our method's operating characteristics over the cancer incidence panel. In Section 7, we give a summary and discuss avenues for future research. We provide technical details in an online supplement. Our R code is freely available upon request.

2. Background

2.1 Event rates on a Lexis diagram

We begin with a brief overview of the Lexis diagram and introduce our notation. ²² A Lexis diagram is a rectangular field with attained age a along the y-axis and calendar time p along the x-axis. For any given population and outcome, individual event times and corresponding person-years at risk are summed beginning at age $a_0$ at calendar time $p_0$ within a grid of square cells defined by A age intervals and P period intervals of width $\mathrm{\Delta }$ . We can describe the structure of Lexis diagrams using the notation $\mathcal{L}\left(a_0,p_0,\mathrm{\Delta },A,P\right)$ . For each cell, the event rate is $y_{ap}=E_{ap}/{\mathcal{O}}_{ap},a=1,\dots ,A,p=1,\dots ,P$ , where $E_{ap}$ is the accumulated number of events and ${\mathcal{O}}_{ap}$ is the corresponding person-years. Each cell can also be referenced by its central value

$$\left(a_{*}\left(a\right),p_{*}\left(p\right)\right)=\left(\left(a_0-\frac{\mathrm{\Delta }}{2}\right)+\mathrm{\Delta }a,\left(p_0-\frac{\mathrm{\Delta }}{2}\right)+\mathrm{\Delta }p\right)$$

We will denote a rate matrix ascertained over a Lexis diagram as ${\mathit{Y}}_{A\times P}$ , and the vectorized version (stacked column-on-column) as $\mathit{y}=\text{vec}\left(\mathit{Y}\right)$ . Our working model posits that the cells contain realizations of independent quasi-Poisson random variables with over-dispersion parameter ${\varphi }^2$ . Given an estimate ${\tilde{\varphi }}^2$ of ${\varphi }^2$ (Section 3.1), the estimated variance–covariance matrix of $\mathit{y}$ is

$$\widehat{\text{Var}}\left(\mathit{y}\right)=V_{\mathit{y}}={\tilde{\varphi }}^2\text{diag}\left(\text{vec}\left({\left[\frac{E_{ap}}{{\mathcal{O}}_{ap}^2}\right]}_{A\times P}\right)\right)\equiv {\tilde{\varphi }}^2{\mathrm{\Sigma }}_{\mathit{y}}$$

2.2 Cancer incidence panel

We extracted authoritative cancer incidence data for the United States from the Surveillance, Epidemiology, and End Results Program's Thirteen Registries Database (SEER-13 ²³ ) for 50 single-years of age (ages 35–84) and 27 calendar years (1992–2018). We selected a panel of representative scenarios defined by cancer site (14 cancers associated with obesity ²⁴ ), sex, and standard SEER race/ethnicity categories (Table 1 and Supplemental Part A). We redistributed breast cancer cases with missing or unknown ER status to the corresponding age- and year-specific ER positive and negative cells using a validated approach.^25,26 The panel includes a total of 1,049,633 incident cases.

2.3 Kernel functions

See Supplemental Part B for details. The standard form of a univariate kernel function ²⁰ is represented by a canonical “shape” function $K^{\left(s\right)}\left(\mu \right)$ that is non-negative and symmetric within the interval $\left[-1,+1\right]$ with ${\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}{K}^{\left(s\right)}\left(\mu \right)du={\int }_{-1}^1{K}^{\left(s\right)}\left(\mu \right)du=1$ . The centered bandwidth-specific kernel is

$$K_{\lambda }^{\left(s\right)}\left(u|u_0\right)=\{\begin{array}{ll}{\displaystyle \frac{1}{{\int }_{-\lambda }^{\lambda }{K}^{\left(s\right)}\left(\frac{v}{\lambda }\right)dv}{K}^{\left(s\right)}\left(\frac{u-u_0}{\lambda }\right)}& {\displaystyle u\in \left[u_0-\lambda ,u_0+\lambda \right]}\\ 0& \text{otherwise}\end{array}$$

where $u_0$ is the center point and $\lambda >0$ is the bandwidth. We will use four standard kernels, the boxcar (box), triangle (trg), Epanechnikov (Epan), and triweight (triwt).

We can extend a univariate kernel for use in two isotropic spatial dimensions by expressing its formula as a function of the Euclidean distance between any given point ${\mathit{u}}_1=\left(x_1,y_1\right)\in {\mathbb{R}}^2$ from a target point ${\mathit{u}}_0=\left(x_0,y_0\right),\text{i.e.,}\Vert {\mathit{u}}_1-{\mathit{u}}_0\Vert =\sqrt{{\left(x_1-x_0\right)}^2+{\left(y_1-y_0\right)}^2}$ . Bivariate kernels have support over the region ${\odot }_{{\mathit{u}}_0}^{\lambda }=\left\{\mathit{u}\in {\mathbb{R}}^2:\Vert \mathit{u}-{\mathit{u}}_0\Vert \le \lambda \right\}$ , that is, a circle of radius $\lambda$ centered at ${\mathit{u}}_0$ . Hence

$$K_{\lambda }^{\left(s\right)}\left(\mathit{u}|{\mathit{u}}_0\right)=\{\begin{array}{ll}{\displaystyle \frac{1}{{\int \int }_{{\odot }_{{\mathit{u}}_0}^{\lambda }}{K}^{\left(s\right)}\left(\frac{\Vert \left(x,y\right)-{\mathit{u}}_0\Vert }{\lambda }\right)dxdy}{K}^{\left(s\right)}\left(\frac{\Vert \mathit{u}-{\mathit{u}}_0\Vert }{\lambda }\right)}& {\displaystyle \mathit{u}\in {\odot }_{{\mathit{u}}_0}^{\lambda }}\\ 0& \text{otherwise}\end{array}$$

3. Kernel function analysis

3.1 Kernel filtration

Our kernel function analysis is motivated by the classic domain filtering approach of image processing ²⁷ which itself is motivated by classic kernel methods in statistics.^20,28,29 For images, bivariate kernels are used to construct mean or low-pass “blurring” filters for purposes of “de-noising” the observed pixel values and/or reducing the impact of high-frequency signals.^30,31 We apply bivariate kernels to Lexis diagrams for the same purposes. A standard kernel function analysis ²⁹ is described in Algorithm 1.

Algorithm 1: kernel function analysis of rates on a Lexis diagram:

Choose a kernel shape box, trg, Epan, or triwt. Choose a bandwidth ${\lambda }_{ap}$ for each cell.

Repeat steps 3–6 for each cell.

Superimpose a masking circle of radius ${\lambda }_{ap}$ over the cell's central value $\left(a_{*}\left(a\right),p_{*}\left(p\right)\right)$ .

Calculate kernel weights for all cells whose mid-points are covered by the masking circle.

Divide each weight in step 4 by their sum (normalization).

Use the normalized weights to calculate a weighted average of the corresponding $y_{ap}$ values within the masking circle.

We will refer to the vector of weighted averages as ${\mathit{y}}_F$ ( $F$ for filter, see below), and the corresponding matrix as ${\left({\mathit{Y}}_{A\times P}\right)}_F\equiv {\mathit{Y}}_F$ . These calculations can be described by a linear operator or filter that can be represented by a sparse matrix $K_F$ such that ${\mathit{y}}_F=K_F\mathit{y}$ . Given an estimate of the over-dispersion parameter ${\tilde{\varphi }}^2$ , it follows that $\widehat{\text{Var}}\left({\mathit{y}}_F\right)={\hat{V}}_F^{\mathit{y}}={\tilde{\varphi }}^2{K}_F{\mathrm{\Sigma }}_{\mathit{y}}{K_F}^{\mathrm{\prime }}$ .

Table 1.

Cancer panel. ^a

Population and tumor					Observed cell counts				Simulation parameters
No.	Sex b	Race/ethnicity c	Site a	Age groups	Min (No. 0)	Mean	Max	Sum	Δ	ϕ2
1	F	NHW	ER− breast	35–84	12.9	62.2	143	84,013	1	1
2	F	NHW	ER + breast	35–84	33.7	292	649	394,190	1	1.5
3	F	NHB	ER− breast	35–84	0 (6)	15.7	48	21,234	1	1
4	F	NHB	ER + breast	35–84	3	32.9	85	44,476	1	1
5	F	NHW	Ovary	35–84	2	34.2	73	46,112	1	1.5
6	F	HIS	Ovary	35–84	0 (30)	5.5	17	7443	1	1.5
7	F	NHW	Corpus uteri	35–84	4	72.7	198	98,112	1	2
8	F	HIS	Rectum	35–84	0 (155)	3.1	16	4147	2	1
9	F	NHW	Thyroid	35–84	1	31.1	94	42,044	1	2
10	F	All	Gallbladder	45–84	0 (6)	12	37	14,268	1	1
11	M	HIS	Colon	40–84	0 (9)	11.2	35	13,604	1	1
12	M	NHW	Colon	35–84	0 (2)	75.6	212	102,017	1	1
13	M	NHW	Esophagus	45–84	0 (4)	15.1	44	16,267	2	2
14	M	All	Kidney	35–84	2	62.6	184	84,526	1	1
15	M	NHB	Kidney	35–84	0 (53)	6.9	30	9289	1	1
16	M	API	Liver	40–84	0 (13)	9	30	10,930	1	1
17	M	NHB	Myeloma	45–84	0 (36)	4.6	17	5020	2	2
18	M	API	Pancreas	45–84	0 (61)	4.6	20	4963	2	1
19	M	NHW	Rectum	35–84	0 (8)	25.8	66	34,790	1	2
20	M	NHW	Stomach	45–84	0 (3)	11.3	32	12,188	2	2

a

See Supplemental Part A for details.

b

Female (F), male (M).

c

Non-Hispanic white (NHW), non-Hispanic black (NHB), Hispanic (HIS), Asian and Pacific Islander (API), all races combined (all).

The evaluation grid is discrete, while the possible bandwidths $\lambda$ are continuous. We will select values of $\lambda$ based on discrete odd integer point values $k=3,5,7,9,\dots$ with corresponding bandwidths $\lambda =\frac{k+1}{2}=2,3,4,5,\dots$ etc. Our approach allows us to apply a different bandwidth depending on cell coordinates. We will use kernels defined by the size of the target Lexis, the kernel shape, and a chosen point-value $k_2$ for the outer two annuli and $k_1$ for all other cells. Hence, we can extend the generic notation $K_F$ using the specific notation: $K_{k_2-k_1\text{ - }shape}$ , or $K_{k_2-k_1\text{ - }shape}^{A\times P}$ when we want to include information about the size of the target Lexis. In what follows we consider the simplest possible kernel, $K_{3-3\text{ - }box}^{A\times P}$ , as a benchmark.

Surprisingly, the bivariate kernels used in Algorithm 1 are invertible or nearly so. It follows that little or no information is discarded by these kernels. Their application allows “smooth” signals to pass through almost unchanged, while at the same time, complex variable signals are downweighted. Although the kernel operators are sparse, their inverses are dense.

We these results in hand we can present our “rule-of-thumb” over-dispersion estimator

$${\tilde{\varphi }}^2\equiv \frac{{\left(\left[I_{AP}-2K_{3-3-box}^{A\times P}+{\left(K_{3-3-box}^{A\times P}\right)}^2\right]\mathit{y}\right)}^{\mathrm{\prime }}{\mathrm{\Sigma }}_{\mathit{y}}^{-1}\left(\left[I_{AP}-2K_{3-3-box}^{A\times P}+{\left(K_{3-3-box}^{A\times P}\right)}^2\right]\mathit{y}\right)}{AP-1}$$

See Supplemental Part C for details.

3.2 Adaptive kernel filtration

We hypothesized that we could regularize the output of a filter $K_F$ by applying additional smoothing. Intuitively, the kernel smoothers blunt local peaks and fill in local valleys in the rate matrix. The extent to which this happens is greater or lesser if the peaks or valleys cover a smaller or larger domain relative to the kernel's bandwidth. This behavior of $K_F$ is entirely characterized by its singular values decomposition (SVD). As illustrated in Supplemental Part D, the singular vectors of $K_F$ represent canonical patterns of peaks and valleys that increase in complexity as the magnitude of the corresponding singular values decreases in magnitude. This led us to consider a truncated singular values decomposition approach. ³² A key question is how many singular vectors to discard. We posited that information theory would provide a good answer. ²¹ By combining these ideas, we developed the adaptive filter described below.

Algorithm 2: Adaptive kernel filtration of rates on a Lexis diagram:

Choose a kernel $K_F$

Construct the SVD of $K_F=UL{V}^{\mathrm{\prime }}=U\text{diag}\left({\sigma }_1,{\sigma }_2,\dots ,{\sigma }_j,\dots ,{\sigma }_{AP}\right)V^{\mathrm{\prime }}$

Center the data using the inverse-variance weighted average, ${\mathit{y}}_0=\mathit{y}-\overline{y}{1}_{AP},\overline{y}={\mathit{y}}^{\mathrm{\prime }}{\mathrm{\Sigma }}_{\mathit{y}}^{-1}\mathit{y}/1_{AP}^{\mathrm{\prime }}{\mathrm{\Sigma }}_{\mathit{y}}^{-1}{1}_{AP}$ . Centering the data stabilizes calculations in the subsequent steps. Calculate the rule-of-thumb over-dispersion estimator ${\tilde{\varphi }}^2$ .

1. 3.

   Consider the sequence of truncated kernels $K_{F:s}=U\text{diag}\left({\sigma }_1,{\sigma }_2,\dots ,{\sigma }_s,0,0,\dots ,0\right)V^{\mathrm{\prime }}$

   4.

   For each value of $s$ :

   Consider s to be the effective degrees of freedom or $edf$

   Calculate the residual vector ${\mathit{r}}_0\left(s\right)={\mathit{y}}_0-K_{F:s}{\mathit{y}}_0$

   Calculate the bias-corrected Akaike information criterion statistic

   $$AI{C}_c^F\left(s\right)={\tilde{\varphi }}^2{\mathit{r}}_0{\left(s\right)}^{\mathrm{\prime }}{\mathrm{\Sigma }}_{\mathit{y}}^{-1}{\mathit{r}}_0\left(s\right)+2s+\frac{2s^2+2s}{AP-s-1},s=1,2,\dots ,AP-2$$

   1. 5.

      The best-fit model uses $K_{F:adp}$ , where $adp=\text{argmi}{\text{n}}_{s}AI{C}_c^F\left(s\right)$ .

   2. 6.

      The fitted values are ${\mathit{y}}_{F:adp}=K_{F:adp}\mathit{y}$ and the estimated variance–covariance matrix conditional on the selected model is $\widehat{\text{Var}}\left({\mathit{y}}_{F:adp}\right)={\hat{V}}_{F:adp}^{\mathit{y}}={\tilde{\varphi }}^2{K}_{F:adp}{\mathrm{\Sigma }}_{\mathit{y}}^{-1}{K_{F:adp}}^{\mathrm{\prime }}$ .

   Use of the bias-corrected penalty term is essential to avoid over-fitting, ²¹ especially for “small” Lexis diagrams with a relatively small total number of cells. The bias-corrected AIC serves as a compromise between the classic uncorrected AIC, which overfits when the number of unknown terms is large relative to the amount of data, and the classic BIC, which is known to be overly conservative producing poorly fitting results. ²¹

   3.3 Model averaging

   We now have three versions of the data: the observed or “raw” data ${\mathit{y}}_0$ , the filtered data ${\mathit{y}}_F$ based on Algorithm 1, which we will also refer to as ${\mathit{y}}_{F:all\text{ - }in}$ because all of the kernel's singular vectors contribute to the filtration, and the adaptive filtration ${\mathit{y}}_{F:adp}$ based on Algorithm 2. We obtain our penultimate version by applying model averaging as reviewed by Burnham and Anderson. ²¹

   Algorithm 3: Model averaging:

2. Apply Algorithm 2 to identify optimal truncation values $ad{p}_1\text{for; }{K}_{F_1,}ad{p}_2\text{for}{K}_{F_2},\dots ,ad{p}_G\text{; for}{K}_{F_G}$ .

3. Identify the optimal kernel $g_0=\text{argmi}{\text{n}}_{g\in \left\{1,2,\dots ,G\right\}}AI{C}_c^{F_g}\left(ad{p}_g\right)$ . Discard all kernels with ${\mathrm{\Delta }}^{F_g}\equiv AI{C}_c^{F_g}\left(ad{p}_g\right)-AI{C}_c^{F_{g_0}}\left(ad{p}_{g_0}\right)>tol$ . By default, use $tol=7$ . Values of $tol$ between 7 and 10 have a theoretical justification. ²¹ Label the H remaining kernels $K_{F_1},K_{F_2},\dots ,K_{F_H}$ , where $1\le H\le G$ .

4. For each kernel identified in step 3, identify all cut-off values c whose $AI{C}_c^{F_h}\left(c\right)$ values are within $tol$ of the global minimum attained by model $g_0$ , that is, the sets ${\mathcal{C}}_h=\left\{c\in 1,\dots ,AP:{\mathrm{\Delta }}^{F_h}\left(c\right)=AI{C}_c^{F_h}\left(c\right)-AI{C}_c^{F_{g_0}}\left(ad{p}_{g_0}\right)\le tol\right\},h=1,\dots ,H$ . Let $n_h=|{\mathcal{C}}_h|$ and $M={\sum }_{h=1}^H{n}_h$ . Denote the complete set of “acceptable” models, which includes $n_1$ members based on $K_{F_1}$ , $K_{F_1:{\mathcal{C}}_1\left\{1\right\}},\dots ,K_{F_1:{\mathcal{C}}_1\left\{n_1\right\}},n_2$ members based on $K_{F_2}$ , $K_{F_2:{\mathcal{C}}_2\left\{1\right\}},\dots ,K_{F_2:{\mathcal{C}}_2\left\{n_2\right\}}$ , etc., as $K_{{\mathcal{F}}_m},m=1,\dots ,M$ . Assemble the corresponding $AI{C}_c$ differentials ${\mathrm{\Delta }}^{{\mathcal{F}}_m}$ and filtrations ${\mathit{y}}_{{\mathcal{F}}_m}=K_{{\mathcal{F}}_m}\mathit{y}$ .

5. Calculate Akaike weights

   $$w_m=\frac{\text{exp}\left(-\frac{{\mathrm{\Delta }}^{{\mathcal{F}}_m}}{2}\right)}{{\sum }_{k=1}^M\text{exp}\left(-\frac{{\mathrm{\Delta }}^{{\mathcal{F}}_k}}{2}\right)},m=1,\dots ,M$$

6. Calculate the truncated model average filtration

   $${\mathit{y}}_{tma}=\sum _{m=1}^M{w}_m{\mathit{y}}_{{\mathcal{F}}_m}$$

7. Calculate the unconditional (model-averaged) variance-covariance matrix

   $$\widehat{\text{Var}}\left({\mathit{y}}_{tma}\right)={\hat{V}}_{tma}^{\mathit{y}}={\tilde{\varphi }}^2\sum _{m=1}^M{w}_m\left(K_{{\mathcal{F}}_m}{\mathrm{\Sigma }}_{\mathit{y}}{K_{{\mathcal{F}}_m}}^{\mathrm{\prime }}\right)+\sum _{m=1}^M{w}_m\left({\mathit{y}}_{{\mathcal{F}}_m}-{\mathit{y}}_{tma}\right){\left({\mathit{y}}_{{\mathcal{F}}_m}-{\mathit{y}}_{tma}\right)}^{\mathrm{\prime }}$$

There is no penalty other than computation time for using a comprehensive panel of kernels. We used 48 kernels in our case study and simulations (Supplemental Part B).

3.4. Variance calculations

Each algorithm produces a variance-covariance matrix that allows us to set confidence limits for any given cell, or for any set of linear combinations of cells. Let $\mathcal{F}$ denote a filtration; formally, use $\mathcal{F}=0$ to denote the observed data. Then ${\mathit{y}}_{\mathcal{F}}$ is the filtered data, and ${\hat{V}}_{\mathcal{F}}^{\mathit{y}}$ is the corresponding variance-covariance matrix ( ${\hat{V}}_{\mathcal{F}}^{\mathit{y}}={\hat{V}}_F^{\mathit{y}}$ for Algorithm 1, ${\hat{V}}_{\mathcal{F}}^{\mathit{y}}={\hat{V}}_{F:adp}^{\mathit{y}}$ for Algorithm 2, and ${\hat{V}}_{\mathcal{F}}^{\mathit{y}}={\hat{V}}_{tma}^{\mathit{y}}$ for Algorithm 3). Matrices ${\hat{V}}_F^{\mathit{y}}$ and ${\hat{V}}_{F:adp}^{\mathit{y}}$ are conditional variances, that is, conditional on the selected kernel, but ${\hat{V}}_{tma}^{\mathit{y}}$ is unconditional, that is, it incorporates uncertainty arising from the model selection process.

Given C linear combinations of ${\mathit{y}}_{\mathcal{F}}$ encoded in contrast matrix ${\mathcal{O}}_{C\times AP}$ , the vector of outputs is ${\mathit{z}}_{\mathcal{F}}=\mathcal{O}{\mathit{y}}_{\mathcal{F}}$ and ${\widehat{Var}}_{{\mathit{z}}_{\mathcal{F}}}={\hat{V}}_{\mathcal{F}}^{\mathit{z}}=\mathcal{O}{\hat{V}}_{\mathcal{F}}^{\mathit{y}}{\mathcal{O}}^{\mathrm{\prime }}$ .

To analyze log-transformed rates ${\mathit{\nu }}_{\mathcal{F}}=\text{log}\left({\mathit{y}}_{\mathcal{F}}\right)$ , the delta-method variance-covariance matrix is $\widehat{\text{Var}}\left({\mathit{\nu }}_{\mathcal{F}}\right)={\hat{V}}_{\mathcal{F}}^{\mathit{\nu }}={\hat{V}}_{\mathcal{F}}^{\mathit{y}}\circ {\left({\mathit{y}}_{\mathcal{F}}{{\mathit{y}}_{\mathcal{F}}}^{\mathrm{\prime }}\right)}^{{\circ }^{-1}}$ , where ^∘ denotes the Hadamard (elementwise) product and ${(.)}^{{\circ }^p}$ denotes the elementwise $p$ -th power. Application of a contrast matrix ${\mathcal{O}}^{*}$ to the log rates ${\mathit{\nu }}_{\mathcal{F}}$ yields ${\mathit{\zeta }}_{\mathcal{F}}={\mathcal{O}}^{*}{\mathit{\nu }}_{\mathcal{F}}$ and ${\hat{V}}_{\mathcal{F}}^{\mathit{\zeta }}={\mathcal{O}}^{*}{\hat{V}}_{\mathcal{F}}^{\mathit{\nu }}{{\mathcal{O}}^{*}}^{\mathrm{\prime }}$ .

4. Results

4.1 Case study: Invasive ER− breast cancer incidence among non-Hispanic white women

Invasive female breast cancer is the most common malignancy among women, and its epidemiology varies by tumor subtype. ³³ Two major subtypes are ER-positive (ER+) and ER-negative (ER−) tumors. ³³ For reasons that remain unclear, the incidence of ER− breast cancer has been decreasing over time in many populations around the world, and the rate of decrease over time varies by age.^5,34–37 In this case study, we apply our new methods to SEER data on invasive ER− breast cancer incidence among non-Hispanic white women (Table 1, No. 1).

Figure 1(A) presents a heat map of the observed rates per 100,000 women years. There is considerable variability from cell to cell, but the data are not over-dispersed ( ${\tilde{\varphi }}^2={1.06}^2$ ). Figure 1(B) presents a canonical plot of log incidence over time within 5-year age groups, and Figure 1(C) plots the cross-sectional log incidence by age within 5-year calendar periods. Incidence appears to be decreasing over time within each age group. However, overlapping values and error bands make it hard to distinguish between many of the curves, and it is difficult to discern in a quantitative sense how fast the rates are changing within and between age groups and calendar periods.

Figure 1.

Estrogen receptor-negative breast cancer incidence among non-Hispanic white women. Raw data (panels A–C), benchmark kernel (panels D – F), and truncated model average (panels G – I). Left panels: Lexis diagram heat maps. Center panels: Rates over time within 5-year age groups. Right panels: Rates by age within 5-year calendar periods. Shaded envelopes show 95% point-wise confidence limits.

Figure 1(D) to (F) repeats these analyses using outputs from the benchmark all-in kernel. Panel 1D presents a heat map based on ${\mathit{y}}_{3-3-box:all\text{ - }in}$ , and panels 1E-F display log-scale curves obtained from ${\mathit{\nu }}_{3-3-box:all-in}$ with confidence bands set using ${\hat{V}}_{3-3-box:all\text{ - }in}^{\mathit{\nu }}$ . The patterns are clearer, the curves more easily distinguished, and one can appreciate that age incidence is decreasing over time in every age group, especially among women 45–49 and 50–54 years old. Compared to panels 1B to C, the confidence limits are 55% narrower on average.

Figure 1(G) to (I) repeats these analyses using the filtered data ${\mathit{y}}_{tma}$ and ${\mathit{\nu }}_{tma}$ , with log-scale confidence bands set using ${\hat{V}}_{tma}^{\mathit{\nu }}$ . Algorithm 3 identified $H=6$ kernels and a total of $M=31$ acceptable truncation values. The most heavily weighted kernels were $K_{3-3-triwt:51}$ and $K_{5-7-triwt:50}$ (Akaike weights of 0.20 and 0.19, respectively). The patterns are clearest, the curves readily distinguished, and the confidence bands are 71% narrower. The TMA kernel is illustrated in Supplemental Figure D4.

4.2 Case study: Averages, gradients, and trends

Graphs provide insight, but no descriptive study is complete until salient features that may be apparent in graphs are quantified using objective and reproducible statistics. There is latitude regarding the precise feature set^8,10; widely used features obtain from averages, gradients, and trends. We consider five features that together quantify the marginal effects of period and age as well as interactions between period and age.

The marginal period curve $mp{c}_{\mathcal{F}}^{\mathit{\nu }}$ and its gradient $\mathrm{\partial }mp{c}_{\mathcal{F}}^{\mathit{\nu }}$ . The marginal period curve $mp{c}_{\mathcal{F}}^{\mathit{\nu }}$ is the log rate averaged over the age group, one value for each calendar period. The gradient $\mathrm{\partial }mp{c}_{\mathcal{F}}^{\mathit{\nu }}$ is the first difference of $mp{c}_{\mathcal{F}}^{\mathit{\nu }}$ divided by the bin width $\mathrm{\Delta }$ .

The marginal age curve $ma{c}_{\mathcal{F}}^{\mathit{\nu }}$ and its gradient $\mathrm{\partial }ma{c}_{\mathcal{F}}^{\mathit{\nu }}$ . The marginal age curve $ma{c}_{\mathcal{F}}^{\mathit{\nu }}$ is the log rate averaged over period, with one value for each age group. The gradient $\mathrm{\partial }ma{c}_{\mathcal{F}}^{\mathit{\nu }}$ is the first difference of $ma{c}_{\mathcal{F}}^{\mathit{\nu }}$ divided by the bin width $\mathrm{\Delta }$ .

The slope of the age-specific log rates over time, ${\ell }_{\mathcal{F}}^{\pi |\alpha }$ , one slope for each age group.

Each feature is a linear function of the log rates; therefore, each can be extracted using a corresponding contrast matrix ${\mathcal{O}}^{*}$ whose values depend only on the structure of the corresponding Lexis diagram $\mathcal{L}\left(a_0,p_0,\mathrm{\Delta },A,P\right)$ . See Supplemental Part E for details.

Figure 2(A) to (C) presents these features calculated from the observed data ${\mathit{\nu }}_0$ , Figure 2(D) to (F) from ${\mathit{\nu }}_{3-3-box:all-in}$ , and Figure 2(G) to (I) from ${\mathit{\nu }}_{tma}$ . Looking down the columns, the estimates become increasingly regular and the confidence limits narrower. The marginal period curve based on model averaging ( $mp{c}_{tma}^{\mathit{\nu }}$ , Figure 2(G), left axis) decreases over time, by around 2% per year through 2004 ( $\mathrm{\partial }mp{c}_{tma}^{\mathit{\nu }}$ , Figure 2(G), right axis), and by around 3% per year circa 2005–2013. Subsequently, the gradient approaches 0 as the curve approaches a plateau.

Figure 2.

Breast cancer averages, gradients, and trends. Features extracted from the data are shown in Figure 1. Raw data (panels A–C), benchmark kernel (panels D – F), and truncated model average (panels G – I). Left panels: Marginal period curve (left axis) and gradient (right axis). Center panels: Marginal age curve (left axis) and gradient (right axis). Right panels: Age-specific period trends. Gradient estimates in the left and center panels are trimmed to exclude the first and last time points. Shaded envelopes show pointwise 95% confidence limits.

The marginal age curve ( $ma{c}_{tma}^{\mathit{\nu }}$ , Figure 2(H), left axis) increases continuously until circa age 72 at a decreasing rate ( $\mathrm{\partial }ma{c}_{tma}^{\mathit{\nu }}$ , Figure 2(H), right axis). Subsequently, the curve attains a plateau and then slowly decreases.

The incidence decreases over time in every age group at a rate that varies by age ( ${\ell }_{tma}^{\pi |\alpha }$ , Figure 2(I)). The greatest decreases of around 3.25% per year occurred circa ages 45–54, consistent with the curves shown in Figure 1(H).

4.3 Simulation studies

We assessed the operating characteristics of our algorithms by simulating the 20 cancers in Table 1. See Supplemental Part F for details. In five scenarios, we blocked the data into $2\times 2$ cells, and in nine, we incorporated over-dispersion ( ${\varphi }^2>1$ ). We considered three approaches: raw data; benchmark all-in kernel (Algorithm 1 with $K_{3-3-box}^{A\times P}$ ); and truncated model average (TMA; Algorithm 3 with 48 candidates). We simulated 500 Lexis diagrams for each scenario. We tracked the performance of benchmark and TMA for heat maps (Figure 1(D) and (G)) and features (Figure 2(D) to (I)). To summarize performance, we calculated the percentage reduction in average root mean squared error ( $RMSE$ ) for benchmark versus raw and for TMA versus raw. We also examined bias and variance per cell (heat maps) and per abscissa value (features).

In every scenario, TMA produced more accurate heat maps than the benchmark kernel (Figure 3(A)). Compared to analyzing the raw data, the benchmark kernel reduced the mean percentage error by 58% on average over the 20 scenarios, versus 74% for TMA.

Figure 3.

Arrow plots of simulation results. Rows correspond to 20 cancers summarized in Table 1. Panels correspond to features. Blue circles show a percent reduction for benchmark kernel versus raw, and yellow triangles show a percent reduction for truncated model average (TMA) versus raw.

Performance gains for $mpc$ were comparatively modest (Figure 3(B)). In Scenario 15 (M Kidney NHB), the benchmark kernel performed slightly better. TMA was uniformly superior to the benchmark for estimating the gradient $\mathrm{\partial }mpc$ (Figure 3(C); 85% mean reduction versus 68%), and TMA achieved similar gains for $mac$ and $\mathrm{\partial }mac$ (Figure 3(D) and (E), respectively). TMA was uniformly superior for estimating the slopes ${\ell }^{\pi |\alpha }$ (Figure 3(F); 61% vs. 46%).

How did TMA achieve such gains in performance? The rule-of-thumb over-dispersion estimator, estimated in step 2 and used in step 4 of Algorithm 2, appears adequate (Figure 4(A)). On average over the 20 scenarios, the median value of ${\tilde{\varphi }}^2$ (squares) exceeded the true value (triangles) by just 0.043 or 3.1%, and the ninety percent error bars for ${\tilde{\varphi }}^2$ varied by just $\pm 0.16$ ( $\pm 11.0\mathrm{\%}$ ) around the corresponding median (squares).

Figure 4.

Operating characteristics of truncated model average. (A) Rows correspond to 20 cancers summarized in Table 1. Median (squares) and 90% limits (bars) of the estimated over-dispersion parameter ( ${\tilde{\varphi }}^2$ ), versus true values (triangles). (B) Box plots of effective degrees of freedom ( $edf$ ). (C) Best-fit kernels, frequency of selection across 20 cancers.

In all scenarios, $edf$ values were one or two orders of magnitude smaller than the maximum possible $edf$ of $AP-2$ (Figure 4(B)). Hence, truncated kernels were always superior to any fixed all-in kernel regardless of kernel shape. Small values of $k_2$ and $k_1$ were selected by TMA a large majority of the time (Figure 4(C)), and all four shapes contributed. Interestingly, amongst the four kernel shapes, the boxcar was selected the least (6.5%) and triweight the most (73%).

5. Discussion

We developed a novel non-parametric approach to regularize (“smooth”) event rates ascertained over Lexis diagrams. Our methods borrow smoothing concepts from time series analysis, that is, k-point moving averages, and classic multivariate kernel methods in statistics ²⁹ and image processing, that is, filtering and singular values decomposition. Our approach uses statistical information theory, specifically, the bias-corrected AIC, ²¹ to handle the model selection problem and provide a variance–covariance matrix that takes model selection into account.

Our truncated model averaging approach adds to the armamentarium and complements existing non-parametric^12,13,17 and Bayesian ¹⁵ methods. The kernel-based methods developed by Duong and Hazelton ²⁹ are closest in spirit to our approach. In our approach, we use bias-corrected AIC versus cross-validation, and we use truncated kernels versus full kernels, which, as we show, substantially increases accuracy.

Our approach to selecting effective degrees of freedom to quantify the complexity of the underlying Lexis diagram is similar to approaches used to select basis functions in generalized additive models. ³⁸ Furthermore, Gaussian process-based smoothing can be viewed as special case of our new bivariate methods since the Matérn covariance used in previous work ³⁹ encompasses both the Exponential and Gaussian kernels, both of which are similar to the triweight kernel. Interestingly, in our simulations, the triweight kernel was by far the most frequently selected kernel shape.

Examination of smoothed Lexis diagrams provides a good overview. Subsequently, scientific conclusions obtain from quantifications of averages, gradients, and trends. As illustrated by our case study, not only does the truncated model averaging approach provide appealingly smoothed heat maps (Figure 1), but corresponding linear combinations of the smoothed values are also substantially more precise (Figure 2).

As shown by our simulation studies, when the cells of the observed Lexis diagrams are statistically independent quasi-Poisson variates, Algorithm 3 (truncated model averaging) is superior to any fixed kernel. Indeed, compared to our benchmark kernel, the truncated model average reduced the intrinsic root mean squared error by 60%–86% depending on the feature.

We focused here on Lexis diagrams, but our methods and software can be applied when the data follow an approximate multivariate normal distribution with a full-rank covariance matrix known up to a scale parameter. Sparse cell counts are an issue in the Poisson case. We implemented our methods using the Normal approximation to the Poisson distribution, which worked well in all scenarios considered (Table 1). If single-year data are sparse (numerous zeros) one can increase the bin width. Future research might consider incorporating a standard or zero-inflated Poisson log-likelihood function. More advanced kernels could also be investigated, for example, steering kernels. ²⁷ Indeed, our model averaging approach could be extended to include estimates obtained using complementary methods such as steering kernels or splines.

In cancer surveillance research few studies examine a single Lexis diagram in isolation. Rather, hypotheses are generated by examining related Lexis diagrams defined by sex, race/ethnicity, geographic region, tumor characteristics, etc. Invariably the amount of information per Lexis diminishes with increasing stratification. As demonstrated here, our new truncated model averaging approach can advance such descriptive studies.

In future work, truncated model averaging might be extended to incorporate as covariates the effects of sex, race/ethnicity, etc. For example, one could start with a joint parametric fit, for example, a proportional hazards age-period-cohort model, ⁴⁰ then apply truncated model averaging to the residuals to characterize any lack-of-fit.

Supplemental Material

sj-docx-1-smm-10.1177_09622802231192950 - Supplemental material for Smoothing Lexis diagrams using kernel functions: A contemporary approach

Supplemental material, sj-docx-1-smm-10.1177_09622802231192950 for Smoothing Lexis diagrams using kernel functions: A contemporary approach by Philip S Rosenberg, Adalberto Miranda Filho, Julia Elrod, Aryana Arsham, Ana F Best and Pavel Chernyavskiy in Statistical Methods in Medical Research