Nonparametric Estimation of Repeated Densities with Heterogeneous Sample Sizes

Let $\left\{X_{ij}:j=1,\dots,N_{i},i=1,\dots,n\right\}$ be the observations of a continuous variable in subpopulation $i=1,\dots,n$; e.g. the age-at-admission to the Intensive Care Unit (ICU) as described in Subsection 4.1 grouped by primary diagnosis. Here $X_{i1},\dots,X_{iN_{i}}$ are i.i.d. observations from the $i$th subpopulation with density $p_{i}$. Densities $p_{1},\dots,p_{n}$ are modeled as i.i.d. realizations of a random density $p$ taking values in family $\mathcal{P}$, where $\mathcal{P}$ consists of densities supported on a common compact interval $\mathcal{T}\subset\mathds{R}$. For example, there may exist latent characteristics that determine the susceptibility age profiles of different diseases, so the age distribution of a disease can be regarded the realization of a density-valued random element. To emphasis the grouping structure, we refer to $p_{i}$ as the density of the $i$th subpopulation, while the distribution of the random density $p$ as the population.

Our goal is to estimate the density $p_{0}$ in a subpopulation with a few i.i.d. observations $X_{0j}$, $j=1,\ldots,N$, where the number of observations $N$ can potentially be small. For example, this problem occurs when one want to estimate the age distribution for an uncommon disease. We will first obtain low-dimensional structures of $\mathcal{P}$ from subpopulations where ample observations are available, e.g. from records for more common conditions, and apply the structure to estimate the density in the subpopulation with a small sample size.

We first set the density families to be either the truncated normal or the bimodal distributions. For the truncated normal scenario, each sample was generated from $N(\mu,\sigma^{2})$ truncated on $[-3,3]$ with $\mu\sim\text{Unif}(-2,2)$ and $\sigma\sim\text{Unif}(2,4)$. For the bimodal scenario, each sample was generated according to density $p(x)\propto\exp((4+\theta)x-(26.5+\theta)x^{2}+47x^{3}-25x^{4}),x\in[0,1]$ with $\theta\sim\text{Unif}(0,10)$. In each Monte Carlo experiment, $n=50$ training samples of size $N_{i}=200$ and $n^{*}=100$ testing samples with sizes $N=10$ or 50 were created.

The mean KL divergence in Figure 3.1 shows that the proposed FPCA_MLE using either KDE as pre-smoother or with the complete training densities observed both produce estimates comparable to MLE without requiring explicit specification of a parametric family. Among all the non-parametric density estimators, the proposed methods are shown to be stable and best-performing, especially when the testing samples are small, in which case the proposed methods outperformed MLE by borrowing strength among different samples, whereas MLE does not pool information. In particular, for truncated normal family with $N^{*}=10$, on average the information loss of the proposed MAP and BLUP were around 20% less compared to that of MLE within ground truth family, and more than 40% less compared to that of KDE.

Contrary to the previous scenario where the underlying families were exponential families, we next estimate densities from a family consisting of mixtures of Gaussian distributions, which forms a multi-modal non-exponential family. This is to further demonstrate the flexibility of the proposed methods to a “mis-specified” case and also demonstrate the data-adaptive selection of the number of components.

The samples were generated following a mixture of three Gaussian distributions with density $p(x)\propto\sum_{l=1}^{3}\theta_{l}\phi((x-\mu_{l})/\sigma_{l})$ truncated to $x\in[-3,3]$, where $\phi$ is the standard normal density and the mixture probabilities $(\theta_{1},\theta_{2},\theta_{3})$ were generated from $\text{Dirichlet}(1/3,1/3,1/3)$. For $l=1,2,3$, we generated the mixture mean $\mu_{l}$ from $\text{Unif}(-5,5)$ and standard deviation $\sigma_{l}$ from $\text{Unif}(0.5,5)$. Here, instead of performing a direct MLE which is known to be unstable for mixture distributions, the expectation–maximization (EM) algorithm in the variant of Lee and Scott (2012) is used as the baseline parametric density estimator.

The left panel of Figure 3.2 compares the relative error of various density estimates to that of the EM algorithm as a baseline reference. The proposed methods enjoyed the least information loss for all but the largest sample size $N^{*}=500$ when EM method became the best. When the testing sample size $N^{*}=25$, on average, the information loss of proposed MAP and BLUP were 20% less than that of RPP, and they were all at least five times better compared to the EM algorithm. Remarkably, this is achieved by only utilizing the first few components on average, as demonstrated in the right panel of Figure 3.2. This shows that the proposed methods are effective alternatives to the EM algorithm, which overcome local maximum issues by performing convex optimization within approximating exponential families. In addition, by using the proposed AIC method for component selection, the proposed estimates utilized more components as the size of fitting sample increased, suggesting that the AIC selection adapts to the increased complexity of the underlying family as unveiled by a larger sample. Indeed, the larger the sample size, the more representative the sample will be for the underlying distribution, which has 8 free parameters in our setup. By using more components, the proposed methods can accommodate the increased complexity shown by the data such as multi-modality, which may otherwise be hidden when fitting with fewer components.

We consider the age-at-admission distributions of patients with critical conditions in the Medical Information Mart for Intensive Care (MIMIC) data (Johnson et al. 2016), collected from patients admitted to critical care units at a large hospital in Boston between 2001 and 2012. The number of patients under different primary diagnosis are highly different. Ample data points are available to establish a precise non-parametric density estimate for common diseases, as shown in the left panels of Figure 4.1, whereas observations for uncommon conditions are scarce and cannot afford flexible non-parametric approaches such as KDE, as displayed in right panels of Figure 4.1. Estimating the densities in the hard-to-observe subpopulations of patients with uncommon conditions is the emphasis in our analysis.

Gender and diagnosis code (ICD9) combinations were used to define the subpopulations, where we excluded birth-related, unclassified, or unspecified diagnosis. A total of $940$ subpopulations with more than 5 observations each were analyzed, among which $75$ subpopulations with at least $75$ observations were used as the training while the rest as testing. For example, the upper left panel in Figure 4.1 shows a training sample formed by $109$ female patients with juvenile type diabetes with ketoacidosis (DMI ketoacd uncontrold), while the bottom right panel shows a testing sample consisting of 27 male patients with abdominal aortic aneurysm rupture (Rupt abd aortic aneurysm). We focused on estimating the densities supported in the age range of $[14,85]$. KDE was utilized as the pre-smoother and the bandwidth was set to the median of the bandwidths in the training samples selected by the method of Sheather and Jones (1991).

As illustrated in right panels of Figure 4.1, the proposed density estimates produced sensible results given small test samples, which were smoother than the estimates given by KDE and avoided multiple spurious modes.

The modes of variation in the density functions obtained by varying each of the first three natural parameters of $\hat{\mathcal{P}}_{K}$ are shown in the left panels of Figure 4.2. The first three modes of variations correspond to a shift in the age distribution towards the old, the elderly, and the middle age, respectively. This suggests the existence of latent disease characteristics that affect the age profiles of different conditions.

The effects of shrinkage in the proposed MAP and BLUP estimators are demonstrated in the right panels of Figure 4.2 for the four rare conditions displayed in the right panels of Figure 4.1. Subpopulations with smaller sample sizes were typically shrunk more towards the population mean using the MAP and BLUP methods that pool population information from the training samples, in which case the estimate differed more notably from the MLE within the approximating family.

To numerically assess the density estimates, we consider cross-entropy which also targets the information loss instead of the KL divergence, since the underlying densities are unknown. Given any density estimator $\hat{p}$, we approach the cross-entropy $H(p,\hat{p})\overset{\triangle}{=}-E_{p}\log{\hat{p}}$, which is the quantity in $D(p\;\|\;\hat{p})$ that depends on the estimate $\hat{p}$, by an unbiased leave-one-out estimate

where $\hat{p}_{-j}$ is a density estimate constructed with all but the $j$th observations. A summary of cross-entropy applied on 844 testing samples using the proposed methods and KDE are reported in Table 1, where a smaller cross-entropy indicates a better estimate. Occasionally KDE resulted in non-finite cross entropy, leading to a smaller number of samples reported.

Table 1 shows that the proposed methods produced density estimates that better reflect sample information compared to KDE by having smaller mean cross-entropy, and also smaller standard deviations. To quantify the sensitivity w.r.t. the training subpopulations, we partitioned the training samples into 5 folds, held out one fold at a time in the training step, fitted the densities for the testing samples, and reported the range of the five mean cross-entropies for a sensitivity measure. The proposed methods are shown to be insensitive to varying training subpopulations, which suggests that the selected training samples are sufficient for our methods to capture the major structure of the underlying family.

We then analyze the precipitation data in the Haihe region to illustrate the versatility of the proposed methods in providing good estimates for both the body part of the density function and the tail probability. We demonstrate that our proposed methods can borrow information from stations with longer records to construct approximate distribution families for annual maxima of hourly rainfall, which lead to better estimates of densities and return levels for those stations having shorter records.

Predicting the likelihood of extreme precipitation is of great importance because of the serious damage such events can cause to economy and social life. We are in particular interested in heavy rainfall during short time period. The likelihood of extreme precipitation event in a certain region can be quantified using the $T$-year return level, which is defined as the level expected to be exceeded once in $T$ years, or equivalently the $1-1/T$ quantile of the marginal distribution under stationarity and independence assumptions. Sample quantile is a straight forward estimate for small $T$, which is referred to as the empirical return level, but this approach requires more than $T$ data points being available. Otherwise, distribution models are necessary for estimating the return level for longer return periods.

We focus on estimating the 5- to 30-year return levels at various locations in the Haihe River basin, a vital region in northern China that includes mega-cities such as Beijing and Tianjin. Precipitation data from 232 weather stations around Haihe River basin from 1961 to 2012 are available from the National Meteorological Information Center (NMIC) and the China Meteorological Administration (CMA) and are used for this analysis. All the 232 stations investigated have at least 30 years of records.

Our target is the distribution of the annual maximum rainfall per hour, where each station defines a subpopulation. As the distribution of precipitation is heavily right-skewed, the pre-smoothing density estimates could be zero within the working domain, creating an issue when we calculate the log-density. To address this issue and the positive constraint, log-precipitation is used as the working variable. See LABEL:supp-appendix:log.inf of the Supplementary Materials for detail.

We compared return levels estimation and the information loss of the proposed methods with standard approaches. To imitate the situation of newer weather stations having shorter records, we randomly selected $n=50$ stations as training sites and retained the full records in these stations for training, and another $n^{*}=100$ testing sites where records from randomly chosen 10 consecutive years were used for fitting. Estimates of density and 5- to 30-year return levels were obtained from non-parametric approaches including our proposed methods and standard KDE, as well as parametric approaches of maximum likelihood estimates within the families of Gamma and GEV distribution. The performance metric used to evaluate density estimation was cross-entropy (4.1), which measures the overall fit. For return level estimates, we used the relative difference of the estimates to the empirical return levels that are obtained from full records in the testing stations, which measures the fit in the right tail. The entire process was repeated $10$ times, and the overall average cross-entropies and the relative differences of the return levels averaged over the stations are reported in Figure 4.3.

The proposed methods in the MLE, MAP, and BLUP variants are all significantly better than the KDE and GEV approaches under the information loss. The estimates using Gamma family of distributions were comparable to our results but underperformed the MAP and BLUP methods which borrow information from the population. For return level estimation, the proposed MAP and BLUP methods had the least variation compared to MLE, Gamma and GEV methods, where the last three methods do not pool information. For shorter (5- and 10- year) return levels, the bias of all methods were comparable, but for longer (20- and 30- year) return levels, the bias of our methods were significantly smaller than Gamma and GEV methods. In summary, our proposed methods have the best performance in estimating the return levels as well. Hence, the proposed methods provide better fits compared to classical approaches in both the bulk and the tail parts of the distributions and lead to better quantification of interesting characteristics of the distribution such as return levels from a small number of discrete observations.

These numerical results come as a positive surprise, as we only used 10 years of records at the testing sites and yet were able to obtain good estimates for 20- and even 30-years return levels. The proposed approaches are able to estimate the tail structure from the data rather than relying on assumptions, and borrow information from the population enforced through parameter shrinkage.

The proposed MLE estimates within low-dimensional approximating family $\hat{\mathcal{P}}_{K}$ are shown to be consistent. The MLE is computed with an additional i.i.d. sample $X_{1},\dots,X_{N}$ of size $N$, of which the density $p$ is an independent copy of the random density. The approximating family $\hat{\mathcal{P}}_{K}$ is obtained from $n$ discretely observed training subpopulations each with $N_{i}\geq m$ observations. Theoretical results for this practical scenario are presented here, and results for MLE obtained within $\mathcal{P}_{K}$ and $\tilde{\mathcal{P}}_{K}$ are included in the Supplementary Materials. For brevity, we use $\hat{p}=\hat{p}_{N,K}\in\hat{\mathcal{P}}_{K}$ to denote the MLE density estimate with parameter $\hat{\theta}=\hat{\theta}_{N,K}$ defined in (2.14).

The consistency results rely on the following conditions imposed on the covariance of the log-transformed random density.

1. (C1)

   For process $f=\operatorname{\psi}p\in L^{2}(\mathcal{T})$, the eigenvalues and eigenfunctions of the integral operator $\mathcal{G}$ associated with the covariance function $G$ satisfy, as $J\rightarrow\infty$,

   	$\displaystyle\left(\sum_{k=1}^{J}\|\varphi_{k}^{\prime}\|_{\infty}\right)\left(\sum_{k>J}\lambda_{k}\right)$	$\displaystyle=O(1),$
   	$\displaystyle\left(\sum_{k=1}^{J}\|\varphi_{k}^{\prime}\|_{\infty}\right)^{2/3}\left(\sum_{k>J}\lambda_{k}\right)$	$\displaystyle=o(1).$

2. (C2)

   For all $k=1,\dots,\operatorname{rank}\mathcal{G}$, the non-zero eigenvalues are distinct and satisfy $\delta_{k}\overset{\triangle}{=}\inf_{l\not=k}\lvert\lambda_{k}-\lambda_{l}\rvert/2>0$; the estimated eigenfunctions are properly aligned so that $\left\langle\varphi_{k},\tilde{\varphi}_{k}\right\rangle\geq 0$ and $\left\langle\varphi_{k},\hat{\varphi}_{k}\right\rangle\geq 0$.

Condition (C1) ensures that the high frequency eigenfunctions are not overly bumpy relative to the amount of information they are carrying as quantified by the eigenvalues, so that the trajectories can be reasonably represented with the first few components. Condition (C2) is standard for simplifying presentation. We assume $\operatorname{rank}\mathcal{G}=\infty$ in discussion hereafter unless otherwise stated.

Conditions (A1), (C1), and (C2) in combine are sufficient for the consistency of MLE within $\mathcal{P}_{K}$ and $\tilde{\mathcal{P}}_{K}$ as shown in LABEL:supp-thm:kldiv.all of the Supplementary Materials. When discrete training samples and pre-smoothing are involved, in parallel to those in Petersen and Müller (2016), additional conditions (D1), (D2), and (S1) are also required as listed in the Appendix. Those conditions basically requires that the pilot density estimator performs well uniformly, provided the sample sizes are sufficiently large. Here (D1) is a uniform MISE requirement, (D2) requires the $L^{\infty}$ error to be $O_{p}(a_{m})$, and (S1) is a requirement for the rate of sample sizes as the number of samples increases, so that one can derive uniform rates of convergence for the pre-smoothed training densities and log-densities under the $L^{\infty}$ norm. The conditions (D1), (D2), and (S1) can be satisfied by the weighted kernel density estimate (KDE) in (2.11) with proper bandwidths and kernels, e.g., $h\asymp m^{-1/3}$ and the Gaussian kernel. See LABEL:supp-sec:presmooth of the Supplementary Materials for detail.

Another interesting question is the convergence of the component scores of pre-smoothed training trajectories $\check{f}_{i}=\operatorname{\psi}\check{p}_{i}$, $i=1,\ldots,n$, which are used to capture distribution of random densities in MAP (2.15) and BLUP (2.16). Notation-wise, let $\check{f}_{i,K}=\hat{\mu}+\sum_{k\leq K}\check{\eta}_{ik}\hat{\varphi}_{k}$ be the truncated pre-smoothed log-density with $\check{\eta}_{ik}=\left\langle\check{f}_{i}-\hat{\mu},\hat{\varphi}_{k}\right\rangle$ for $k=1,\ldots,K$, $i=1,\dots,n$. The next proposition shows that the estimated component scores $\check{\eta}_{ik}$ converge uniformly to the truth $\eta_{ik}=\left\langle f_{i}-\mu,\varphi_{k}\right\rangle$.

With the previous proposition, it is easy to see $\sup_{i\leq n}\sup_{k\leq K}\lvert\check{\eta}_{ik}-\eta_{ik}\rvert=O_{p}(a_{m}+n^{-1/2})$ as $n\rightarrow\infty$ if $K$ is fixed. To have a similar result with $K\rightarrow\infty$ as well, we need $n,K\rightarrow\infty$ in such a way that $(n^{-1/2}+b_{m})\sup_{k\leq K}\delta_{k}^{-1}\rightarrow 0$. Under the conditions of Proposition 5.1, we can set $K\asymp n^{\beta}$ for any $\beta<\max(1/2,r_{m}/3)/(1+r)$, then the sample component scores are consistent uniformly. In particular, $K=\hat{K}^{*}$ satisfies the requirement.

The following result further characterizes the asymptotic behavior of the proposed MLE for a given new fitting sample with a small and fixed sample size $N$, where $K$ is finite and fixed and there is a large amount of observations from training subpopulations.