An Objective Measure of Distributional Estimability as Applied to the Phase-Type Aging Model

Abstract

The phase-type aging model (PTAM) is a class of Coxian-type Markovian models that can provide a quantitative description of the effects of various aging characteristics. Owing to the unique structure of the PTAM, parametric inference on the model is affected by a significant estimability issue, its profile likelihood functions being flat. While existing methods for assessing distributional non-estimability require the subjective specification of thresholds, this paper objectively quantifies estimability in the context of general statistical models. More specifically, this is achieved via a carefully designed cumulative distribution function sensitivity measure, under which the threshold is tailored to the empirical cumulative distribution function, thus becoming an experiment-based quantity. The proposed definition, which is validated to be innately sound, is then employed to determine and enhance the estimability of the PTAM.

1. Introduction

The phase-type aging model (PTAM) belongs to a class of Coxian Markovian models that were introduced in Cheng et al. (2021). The purpose of the PTAM is to provide a quantitative description of the effects of well-known aging characteristics resulting from a genetically determined, progressive, and irreversible process. It provides a means of quantifying the heterogeneity in aging among individuals and identifying the effects of anti-selection (Cheng 2021; Cheng et al. 2021).

The PTAM has a unique structure, including a constant transition rate for the aging process and a functional form for the relationship between aging and death (Cheng et al. 2021). This structure gives rise to flat profile likelihood functions, putting the reliability of parameter estimates into question even if MLEs can be obtained (Cheng 2021). This problem, which is referred to as the non-estimability issue, was studied in Raue et al. (2009). In this context, many estimated values of m, the number of states, produced nearly the same profile likelihood values, resulting in similar model fitting. Actually, for certain statistical models such as the PTAM, the quality of model fitting is not the only consideration since the model parameters convey biological meaning. If the parameter estimates do not make sense in this context, the research problem cannot be satisfactorily addressed, even if the model fitting quality is found to be adequate. Accordingly, we investigated the problem of flat likelihood functions in the literature, focusing on the field of estimability, also referred to as practical identifiability.

Primarily, estimability must be assessed after identifiability, which can be defined as follows according to Lehmann and Casella (1998):

Definition 1.

Let $\mathcal{M}=\left\{f\right(x;\mathit{\theta }):\mathit{\theta }\in \mathbf{\Theta }\}$ be a statistical model with either finite- or infinite-dimensional parameter space Θ. We say that $\mathcal{M}$ is identifiable if

$$\begin{array}{c}\hfill \left(f(x;{\mathit{\theta }}_{\mathbf{1}})=f(x;{\mathit{\theta }}_{\mathbf{2}})\right)\Rightarrow \left({\mathit{\theta }}_{\mathbf{1}}={\mathit{\theta }}_{\mathbf{2}}\right),\forall {\mathit{\theta }}_{\mathbf{1}},{\mathit{\theta }}_{\mathbf{2}}\in \mathbf{\Theta }.\end{array}$$

If a model is assessed to be non-identifiable, then it proves unnecessary to consider its estimability (Miao et al. 2011). The non-identifiability issue relates to properties of the model structure and can only be removed analytically by adding constraints to the model parameters or proposing new model representations (Hengl et al. 2007; Raue et al. 2009). On the other hand, the non-estimability issue relates to experimental protocols involving data quality (such as being insufficient or too noisy), algorithm approximation, or other noisy measurements. Thus, one must eliminate the non-identifiability issue before assessing estimability, in which case the remaining unreliability of the parameter estimates can be ascribed to experimental protocols (Miao et al. 2011). It has been established that the PTAM is identifiable when the number of states is greater or equal to six (Nie 2022).

We now highlight this paper’s main objectives and contributions with respect to the concept of estimability.

1.1. Defining the Concept of Estimability via an Objective Threshold

In fact, estimability has remained an open problem overall since existing methods for assessing it require the subjective specification of thresholds, which are left to the experimenters to specify. To the best of our knowledge, Gontier and Pfister (2020) seems to be the only paper that investigated and tackled the issue of threshold specification. They proposed a new definition of estimability that is based on a model selection perspective, where the subjective threshold is replaced with Bayes factors. The definition advocated in this paper solves this issue as it is based on a parallel but distinct perspective. Rather than eliminating the threshold, we make it an objective quantity that is related to the experiment itself, thus making it objective.

The proposed definition of estimability bears a resemblance to the profile likelihood function method introduced in Raue et al. (2009), where the non-estimability issue is defined as occurring when the profile-likelihood-based confidence region is infinite up to a subjective threshold. In order to make the threshold experiment-based and, therefore, objective, one needs to develop a methodology that somehow relates the confidence region to the experimental protocol. The proposed definition achieves this by replacing the profile-likelihood-based confidence region with an innovative confidence region that is based on a carefully designed c.d.f. sensitivity measure, whose purpose is to relate the confidence region to the experimental error by quantifying them into single numbers. This enables one to indirectly compare the confidence region and the experimental error. In this case, the threshold is tailored to the experimental protocol and then becomes objective. Since the measure is solely experiment-based, it quantifies the degree of estimability objectively so that the experimenter knows quantitatively the extent to which the experimental design should be improved in order to make the model estimable. The key to our approach is based on the empirical c.d.f. or ECDF; this is also the case for the Anderson–Darling and Kolmogorov–Smirnov tests, which compare a hypothesized distribution with the ECDF. The ECDF constitutes an objective measure as it represents the distribution of the experimental data.

1.2. Extending the Applicability of the Concept of Estimability to Statistical Inference

The concept of estimability originated in the field of biology, where an ODE model was utilized to model dynamic biological systems. It is also interchangeably referred to as “practical identifiability” in the literature. However, very little research on estimability has been conducted in connection with statistical models. Gontier and Pfister (2020) applied the concept of estimability to a statistical model. However, the scope of that paper still remained in system biology as a binomial model was introduced in connection with synaptic transmissions. We address the applicability of the concept of estimability in statistical inference from an objective perspective. Although this is principally illustrated via an application involving the PTAM, the proposed approach is actually applicable in other statistical inference contexts, which are illustrated by means of two additional numerical examples, one involving a discrete distribution and another, a continuous one. This constitutes the second main contribution of this paper.

1.3. Structure of the Paper

This paper is organized as follows. Preliminaries on the PTAM are introduced in Section 2. Section 3 provides a literature review on estimability. Section 4 introduces the proposed definition of estimability, which is validated in Section 5, where it is established to be innately sound. This definition is implemented to assess the estimability of the PTAM via a simulation study in Section 6. Section 7 contains certain remarks and a conclusion.

2. The Phase-Type Aging Model

The phase-type aging model (PTAM) stems from the mortality model introduced by Lin and Liu (2007). Phase-type mortality models enable one to link their parameters to biological and physiological mechanisms of aging, which constitutes a definite advantage. Thus, for instance, the longevity risk facing annuity products can be measured more accurately. Experimental results show that the phase-type mortality model with a four-state developmental period and a subsequent aging period achieved very satisfactory fitting results with respect to the Swedish and USA cohort mortality data set (Lin and Liu 2007). Later on, Su and Sherris (2012) applied the phase-type mortality model to an Australian cohort mortality data set.

In a further study of the PTAM, Cheng et al. (2021) introduced a parsimonious yet flexible representation that allowed for the modeling of various aging patterns. In addition, an efficient algorithm for evaluating the likelihood of the PTAM was developed in Cheng (2023b). The main objective of the PTAM is to describe the human aging process in terms of the evolution of the distribution of physiological ages while utilizing mortality rates as aging-related variables. Thus, although the PTAM can reproduce mortality patterns, it ought not to be treated as a mortality model. In this context, the PTAM is most applicable at human ages beyond the attainment of adulthood, where, relatively speaking, the aging process is the most significant factor that contributes to the variability in lifetimes (Cheng et al. 2021).

2.1. Preliminaries

Definition 2.

Let ${\left\{X\left(t\right)\right\}}_{t\ge 0}$ be a continuous time Markov chain (CTMC) defined on a finite state space $\mathcal{S}=\mathcal{E}\cup \Delta =\{1,2,\dots ,m\}\cup \Delta$, where $\Delta =\{m+1\}$ is the absorbing state and $\mathcal{E}$ is the set of transient states. Let ${\left\{X\left(t\right)\right\}}_{t\ge 0}$ have as initial distribution over the transient states ${\mathit{\pi }}^{\mathbf{\prime }}=({\pi }_1,{\pi }_2,\dots ,{\pi }_m)$ such that ${\mathit{\pi }}^{\mathbf{\prime }}\mathit{e}=1$, and let the transition intensity matrix be

$$\begin{array}{c}\hfill \mathbf{\Lambda }=\left[\begin{array}{cc}\mathit{S}& \mathit{h}\\ {\mathbf{0}}^{\prime }& 0\end{array}\right],\end{array}$$

where $\mathit{S}$ is an $m\times m$ matrix, ${\mathbf{0}}^{\mathbf{\prime }}$ is the transpose of a null vector, and $\mathit{h}=-\mathit{S}\mathit{e}$, with $\mathit{e}$ denoting a vector of ones. Define $T:=inf\{t\ge 0|X\left(t\right)=m+1\}$ as the time until absorption. Then, T is said to follow a continuous phase-type (CPH) distribution, denoted by $CPH(\mathit{\pi },\mathit{S})$ of order m, and $\mathit{h}$ is defined as the exit vector.

There is a long history of using phase-type “absorbing time” distributions as survival models; see Aalen (1995), Asmussen et al. (1996), Lin and Liu (2007), and Su and Sherris (2012), among others.

Result 1.

Given $T\sim CPH(\mathit{\pi },\mathit{S})$ of order m,

• The p.d.f. of T is $f_T\left(t\right)={\mathit{\pi }}^{\prime }{e}^{\mathit{S}t}\mathit{h}$.
• The c.d.f. of T is $F_T\left(t\right)=1-{\mathit{\pi }}^{\prime }{e}^{\mathit{S}t}\mathit{e}$.

It is well-known that, given a CPH distribution of order m, if ${\mathit{\pi }}^{\prime }=(1,0,\dots ,0)$ and $\mathit{S}$ has the structure

$$\begin{array}{c}\hfill \mathit{S}=\left[\begin{array}{ccccccc}-({\lambda }_1+h_1)& {\lambda }_1& 0& 0& \cdots & 0& 0\\ 0& -({\lambda }_2+h_2)& {\lambda }_2& 0& \cdots & 0& 0\\ ⋮& ⋮& ⋮& ⋮& \ddots & ⋮& ⋮\\ 0& 0& 0& 0& \cdots & -({\lambda }_{m-1}+h_{m-1})& {\lambda }_{m-1}\\ 0& 0& 0& 0& \cdots & 0& -h_m\end{array}\right],\end{array}$$

(1)

where ${\lambda }_i>0,h_j>0$, $i=1,2,\dots ,m-1$, and $j=1,2,\dots ,m$, then the resulting distribution constitutes a Coxian distribution with no probability mass at zero Cox (1955a, 1955b). A phase diagram such as that displayed in Figure 1 illustrates the process.

Definition 3.

Given that $T>0$, the PTAM of order m is defined as a Coxian distribution of order m with transition intensity matrix $\mathit{S}$ and exit rate vector $\mathit{h}$ such that

$$\begin{array}{c}\hfill \mathit{S}=\left[\begin{array}{ccccccc}-(\lambda +h_1)& \lambda & 0& 0& \cdots & 0& 0\\ 0& -(\lambda +h_2)& \lambda & 0& \cdots & 0& 0\\ ⋮& ⋮& ⋮& ⋮& \ddots & ⋮& ⋮\\ 0& 0& 0& 0& \cdots & -(\lambda +h_{m-1})& \lambda \\ 0& 0& 0& 0& \cdots & 0& -h_m\end{array}\right],\mathit{h}=\left[\begin{array}{c}{h}_1\\ h_2\\ ⋮\\ h_{m-1}\\ h_m\end{array}\right],\end{array}$$

where $\lambda >0$, $h_m>h_1>0$, and

$$\begin{array}{c}\hfill h_i=\left\{\begin{array}{cc}{\left(\frac{m-i}{m-1}{h}_1^s+\frac{i-1}{m-1}{h}_m^s\right)}^{\frac{1}{s}},\hfill & s\ne 0,\hfill \\ h_1^{\frac{m-i}{m-1}}{h}_m^{\frac{i-1}{m-1}},\hfill & s=0,\hfill \end{array}\right.\end{array}$$

(2)

where $i=1,2,\dots ,m$. This is denoted by $PTAM(h_1,h_m,s,\lambda ,m)$.

As can be seen from Figure 2, the PTAM has a phase diagram similar to that of the Coxian distribution shown in Figure 1 except that the transition rates are constant and the exit rates are functionally related, as specified in (2).

Note the following:

(i)

In Figure 2, each state in the Markov process represents the physiological age—a variable that reflects an individual’s health condition or frailty level. As the aging process proceeds, the frailty level increases until the last state, where the individual’s health conditions have deteriorated to the point of causing death.

(ii)

The transition rate $\lambda$ is assumed to be constant. The exiting rates denoted by $h_i$ are the dying rates or force of mortality. With this setup, an individual is randomly located in a certain state at a given calendar age. This mathematically describes the fact that the individuals involved will have different physiological ages given the same calendar age.

(iii)

The assumption that dying rates have the structure given in (2) is somewhat reminiscent of the well-known Box–Cox transformation introduced by Box and Cox (1964). The first and last dying rates $h_1$ and $h_m$ are included in the model parameters, whereas the remaining in-between rates are interpolated in terms of the parameter s, which is a model parameter related to the curvature of the exit rate pattern. To verify this, Figure 3 and Figure 4 present the effect of s on the pattern of the exit rates for $m=50$ and $m=100$, respectively. When $s=1$, the dying rates have a linear relationship. When $s>1$, the rates are concave, and when $s<1$, the rates are convex. In particular, when $s=0$, the rates behave exponentially. In practice, it is likely that s is less than one when calibrating to mortality data (Cheng et al. 2021). That is, the dying rates increase faster than linearly as an individual ages. Throughout this paper, it is assumed that $h_i$ follows the structure given in (2) for $i=1,2,\dots ,m$.

(iv)

The value of $\lambda$ needs to be commensurate to the value of m; otherwise, there would be no need to have numerous m states when $\lambda$ is small. Accordingly, we let their ratio be a constant; that is,

$$\begin{array}{c}\hfill \frac{m}{\lambda }=\psi ,\end{array}$$

which can be seen as a reparameterization of the PTAM that involves the five parameters $(h_1,h_m,s,\psi ,m)$. We note that such a reparameterization establishes a positive covariance between $\lambda$ and m, which is more in line with the biological interpretation. From now on, this parameterization of the PTAM is utilized.

2.2. Identifiability of the PTAM

The identifiability of the PTAM was established in Nie (2022) for $m\ge 6$, where m denotes the number of states. Illustrative examples of non-identifiable PTAM were also provided for $2\le m<6$ to clarify the concept. The identifiability, or mathematical uniqueness, guarantees that no other c.d.f.-equivalent representations of the model exist and that the likelihood function has a unique global maximum.

However, identifiability does not imply estimability. Although the model representation is unique, parameter estimates can still be unreliable when the profile likelihood functions are extremely flat, as a large range of different estimates can produce nearly identical likelihood values (Raue et al. 2009). According to Cheng (2021), the profile likelihood functions of the PTAM’s parameters $h_1$, $h_m$, and s turn out to be flat, which gives rise to this estimability issue.

3. Literature Review

3.1. Methods for Assessing Estimability

Estimability is also referred to as practical identifiability in the literature. The non-estimability issue arises in the case of flat likelihood functions (Raue et al. 2009) or, equivalently, the insensitivity of the model c.d.f. with respect to its parameters (McLean and McAuley 2012), which may be due to either of the following situations:

(i)

The model c.d.f. is insensitive with respect to parameter changes. Accordingly, this aspect involves model sensitivity.

(ii)

The effect of one parameter on the c.d.f. may be offset by that of one of the other parameters, this being defined as parameter correlation.

If the model is identifiable, then unreliable parameter estimates, if any, are caused by non-estimability issues that may be due to experimental errors, including data quality (insufficient or too noisy), algorithm approximation, or other noisy measurements (Miao et al. 2011). Unlike identifiability, estimability is less well-defined, and its characterization has remained an open problem. While it is straightforward to think qualitatively that parameters can be “loosely estimated” under noisy measurements, one would need to define quantitatively what this really means (Gontier and Pfister 2020; Raue et al. 2009).

Estimability has been widely studied in system biology, where an ODE model is utilized to model dynamic biological systems. For instance, it is assumed that

$$\begin{array}{cc}\hfill \frac{\mathrm{d}\mathit{x}\left(t\right)}{\mathrm{d}t}& =f(t,\mathit{x}(t),\mathit{u}(t),\mathit{\theta }),\hfill \\ \hfill \mathit{y}\left(t\right)& =h\left(\mathit{x}\right(t),\mathit{u}(t),\mathit{\theta })+\mathit{ϵ}\left(t\right),\hfill \end{array}$$

where

$$\begin{array}{cc}\hfill & \mathit{x}\left(t\right)\in {\mathbb{R}}^n\mathrm{is}\mathrm{a}\mathrm{vector}\mathrm{of}\mathrm{state}\mathrm{variables},\hfill \\ \hfill & \mathit{y}\left(t\right)\in {\mathbb{R}}^d\mathrm{is}\mathrm{the}\mathrm{output}\mathrm{vector},\hfill \\ \hfill & \mathit{u}\left(t\right)\in {\mathbb{R}}^p\mathrm{is}\mathrm{the}\mathrm{known}\mathrm{system}\mathrm{input}\mathrm{vector},\hfill \\ \hfill & \mathit{\theta }\in {\mathbb{R}}^q\mathrm{is}\mathrm{the}\mathrm{parameter}\mathrm{vector},\mathrm{and}\hfill \\ \hfill & \mathit{ϵ}\left(t\right)\sim N(\mathbf{0},{\sigma }^2\left(t\right){\mathit{I}}_{\mathit{d}})\mathrm{is}\mathrm{the}\mathrm{measurement}\mathrm{noise}.\hfill \end{array}$$

Four broad types of methodologies are employed. We briefly describe them next.

3.1.1. The Monte Carlo Method

This approach, which was utilized by Aslett (2012) in connection with certain lifetime models, involves repeated parameter estimation from a large number of data sets that were simulated by means of the Monte Carlo method. To apply this method, threshold values for parameter uncertainty levels are required to distinguish estimable from non-estimable parameters. Let ${\mathit{\theta }}_{\mathbf{0}}$ be the nominal parameter vector obtained from fitting the model to the original data or from prior knowledge (Miao et al. 2011). Let ${\widehat{\mathit{\theta }}}_{\mathit{i}}$ be the vector of the parameter estimates at the $i^{th}$ of M trials, which are based on data simulated from the model having ${\mathit{\theta }}_{\mathbf{0}}$ as its parameter vector. Then, the average relative estimation error (ARE) associated with ${\widehat{\theta }}_i^{\left(k\right)}$, the $k^{th}$ element of ${\widehat{\mathit{\theta }}}_{\mathit{i}}$, is given by

$$\begin{array}{c}\hfill AR{E}_k:=\frac{1}{M}\sum _{i=1}^M\frac{|{\theta }_0^{\left(k\right)}-{\widehat{\theta }}_i^{\left(k\right)}|}{|{\theta }_0^{\left(k\right)}|},\end{array}$$

where ${\theta }_0^{\left(k\right)}$ is the $k^{th}$ element ${\mathit{\theta }}_{\mathbf{0}}$, $k=1,2,\dots ,q$. Miao et al. (2011) defined non-estimability as occurring when the ARE of a parameter is sufficiently high or, equivalently, exceeds a pre-selected threshold ${\Delta }_1$.

3.1.2. Methods Based on the Correlation Matrix or the Fisher Information Matrix

According to Petersen et al. (2001), the Fisher information matrix (FIM) associated with the ODE model is given by

$$\begin{array}{c}\hfill FIM=\sum _{i=1}^N{\left(\frac{\partial {\widehat{\mathit{y}}}_i}{\partial \widehat{\mathit{\theta }}}\right)}^{\prime }{\mathit{V}}^{-1}\left(\frac{\partial {\widehat{\mathit{y}}}_i}{\partial \widehat{\mathit{\theta }}}\right),\end{array}$$

(3)

where

$$\begin{array}{cc}\hfill & \left(\frac{\partial {\widehat{\mathit{y}}}_i}{\partial \widehat{\mathit{\theta }}}\right)\mathrm{is}\mathrm{defined}\mathrm{as}\mathrm{the}\mathrm{sensitivity}\mathrm{matrix}\mathrm{and}\hfill \\ \hfill & \mathit{V}\mathrm{is}\mathrm{a}\mathrm{known}\mathrm{positive}\mathrm{definite}\mathrm{matrix}\mathrm{of}\mathrm{weights}\mathrm{on}\mathrm{the}\mathrm{variances}\mathrm{and}\mathrm{covariances}.\hfill \end{array}$$

Rodriguez-Fernandez et al. (2006) proposed a correlation matrix approach for analyzing the estimability of the ODE model. By the Cramér–Rao Theorem, the covariance matrix can be obtained as

$$\begin{array}{c}\hfill \mathit{C}\approx FI{M}^{-1},\end{array}$$

the correlation between ${\theta }_i$ and ${\theta }_j$ being

$$\begin{array}{c}\hfill r_{ij}=\frac{{\mathit{C}}_{ij}}{\sqrt{{\mathit{C}}_{ii}{\mathit{C}}_{jj}}},i\ne j,1\le i,j\le q.\end{array}$$

Similarly, Quaiser and Mönnigmann (2009) proposed a total correlation measure. In this instance, ${\theta }_i$ and ${\theta }_j$ are deemed to be non-estimable if their correlation is sufficiently high or, equivalently, exceeds a certain threshold ${\Delta }_2$.

There exist several methods focusing on the FIM. Dochain and Vanrolleghem (2001) proposed that the condition number, which is defined as the ratio of the largest eigenvalue to the smallest eigenvalue of the FIM, could also be used to assess estimability. The larger the condition number, the more correlated the parameters and the less estimable they are. The model is then non-estimable if the condition number is sufficiently high or, equivalently, exceeds a certain threshold ${\Delta }_3$. In addition, Brun et al. (2001) proposed a collinearity index to measure the parameter correlations. The model is deemed less estimable if the collinearity index is relatively large, which also requires the specification of a threshold.

3.1.3. Methods Based on the Model Sensitivity

Another approach is based on the model sensitivity. As can be seen from (3), the FIM is obtained in terms of the sensitivity matrix. Thus, the sensitivity matrix may be extracted from the FIM and analyzed specifically. The sensitivity matrix $\mathit{S}$ with observation times $(t_1,\dots ,t_N)$ is defined as

$$\begin{array}{c}\hfill {\mathit{S}}_{dN\times q}(t_1,\cdots ,t_N):=\left[\begin{array}{ccc}\frac{\partial y_1(t_1;\mathit{\theta })}{\partial {\theta }_1}& \cdots & \frac{\partial y_1(t_1;\mathit{\theta })}{\partial {\theta }_q}\\ ⋮& \ddots & ⋮\\ \frac{\partial y_d(t_1;\mathit{\theta })}{\partial {\theta }_1}& \cdots & \frac{\partial y_d(t_1;\mathit{\theta })}{\partial {\theta }_q}\\ ⋮& ⋮& ⋮\\ ⋮& ⋮& ⋮\\ \frac{\partial y_1(t_N;\mathit{\theta })}{\partial {\theta }_1}& \cdots & \frac{\partial y_1(t_N;\mathit{\theta })}{\partial {\theta }_q}\\ ⋮& \ddots & ⋮\\ \frac{\partial y_d(t_N;\mathit{\theta })}{\partial {\theta }_1}& \cdots & \frac{\partial y_d(t_N;\mathit{\theta })}{\partial {\theta }_q}\end{array}\right].\end{array}$$

Several methods are based on the sensitivity matrix. Jacquez and Greif (1985) calculated a sample correlation between the matrix columns. If $\rho ({\mathit{S}}_{.i},{\mathit{S}}_{.j})$ is close to one within a certain threshold, then ${\theta }_i$ and ${\theta }_j$ are considered to be non-estimable. Other methods exist, such as the principal components analysis (PCA) technique (Degenring et al. 2004), the orthogonal method (Yao et al. 2003), and the eigenvalue method (Vajda et al. 1989), all of which rely on subjective thresholds.

If the model is sufficiently simplified and involves only a few parameters, the sensitivity function $\frac{\partial \mathit{S}\left(t\right)}{\partial \mathit{\theta }}$ can be solved analytically, in which case the sensitivity matrix is not needed. Holmberg (1982) proposed a visual inspection approach to the sensitivity function. The larger the sensitivity measure of one parameter, the greater the change in the model c.d.f. with respect to the change of that parameter. If the sensitivity functions of certain parameters are linearly dependent, then those parameters are functionally related. The drawback of this approach is that a correlation cannot be quantified based on graphs. Moreover, a visual inspection entails a subjective assessment. Determining whether the graphs of the sensitivity functions are truly dependent depends on the experimenter’s assessment.

3.1.4. Methods Based on Profile Likelihood

Raue et al. (2009) proposed an explicit definition of estimability that is based on the profile likelihood function. They defined the profile likelihood confidence interval for parameter ${\theta }_i$ as

$$\begin{array}{c}\hfill C_{i,{\Delta }^{*}}:=\left\{{\theta }_i\right|PL(\mathit{x};{\widehat{\theta }}_i)-PL(\mathit{x};{\theta }_i)<{\Delta }^{*}\},\end{array}$$

(4)

where $PL(·)$ is the profile likelihood function of ${\theta }_i$ and ${\Delta }^{*}$ is a subjectively chosen threshold. Then, ${\theta }_i$ is said to be non-estimable if $C_{i,{\Delta }^{*}}$ is infinite. In other words, given a certain threshold ${\Delta }^{*}$, there exists a ${\delta }_i>0$ such that for all $|{\theta }_i|>{\delta }_i$, $PL(\mathit{x};{\widehat{\theta }}_i)-PL(\mathit{x};{\theta }_i)<{\Delta }^{*}$ holds true. This definition is mathematically clear as it relies on a binary event: whether $C_{i,{\Delta }^{*}}$ is infinite or not. However, a subjective threshold ${\Delta }^{*}$ is still required, as is the case for the other methods.

3.2. Relationships between Identifiability, Estimability, and Sensitivity

Informally, sensitivity refers to the degree to which a model is affected by its parameter values. Graphically, the more sensitive a model is with respect to one parameter, the more noticeably the model’s c.d.f. is affected by changes in that parameter. In this paper, we quantify this concept by introducing the c.d.f. sensitivity measure that is specified in Definition 4.

Identifiability and estimability are related via the concept of sensitivity. From the perspective of sensitivity, if a statistical model $f(x;\mathit{\theta })$ is non-identifiable with the non-identifiable set being

$$\begin{array}{c}\hfill \mathcal{A}:=\{{\mathit{\theta }}_{\mathbf{1}},{\mathit{\theta }}_{\mathbf{2}}\in \mathbf{\Theta }|\forall x,f(x;{\mathit{\theta }}_{\mathbf{1}})=f(x;{\mathit{\theta }}_{\mathbf{2}})\},\end{array}$$

then f has zero sensitivity if the parameter vector changes within $\mathcal{A}$ since the model c.d.f. does not change. Thus, the non-identifiability issue cannot be overcome by improving the experimental design since the model output is then the same for all x.

On the other hand, if the sensitivity is not zero, then different parameters produce different model c.d.f. values. In this case, different values of x produce different model outputs, and the experiment contributes some information toward parameter inference. The more sensitive the model with respect to one parameter, the more estimable that parameter is. Although, as previously mentioned, different measures, such as those based on correlations, condition numbers, and eigenvalues, may be employed, in each case, a threshold must be set.

It should be emphasized that identifiability and estimability are equally important. If a model is non-identifiable, then the likelihood function has multiple global maxima. In that case, although we are relying on a numerical algorithm that aims to maximize the likelihood function, we may question the reliability of the parameter estimates produced by that algorithm since other maxima may potentially exist. In the case of estimability, although we are certain that the MLEs are unique, we may as well question the reliability of the parameter estimates produced by that algorithm since a wide variety of estimates can yield nearly the same likelihood values.

4. A Methodology for Objectively Assessing Estimability

Based on a thorough review of the literature, very little research on estimability appears to have been conducted in connection with statistical models. To the best of our knowledge, there is only one paper, namely Gontier and Pfister (2020), that studied estimability with respect to a statistical model, wherein a new definition of estimability based on a model selection perspective was proposed. A significant contribution of theirs was the elimination of the subjective threshold in (4) by introducing a Bayes factor into the new definition. Our contribution also aims to address the problem of having to set a subjective threshold. Rather than eliminating it, we shall make the threshold an objective, experiment-based quantity.

In order to do so, a methodology needs to be established that relates the confidence region to the experimental protocol. To achieve this, we rely on the following considerations:

(i)

The curvature of the likelihood function reflects the sensitivity of the model (c.d.f.) with respect to the parameters (McLean and McAuley 2012).

(ii)

The non-estimability issue is defined as occurring when the likelihood-based confidence interval is infinite (Raue et al. 2009).

In regards to the first consideration, we replace the likelihood-based confidence region introduced in Raue et al. (2009) with an innovative confidence region that is based on a carefully designed c.d.f. sensitivity measure. Such a sensitivity measure aims to relate a confidence region to the experimental error by quantifying them. This is indirectly achieved by comparing these quantified numbers under the same measure. In this case, the threshold is tailored to the experimental protocol and, thus, becomes objective. Additionally, in light of the second consideration, we define non-estimability as occurring when the proposed confidence region is infinite.

Several preliminary definitions are needed before defining estimability:

Definition 4.

Consider an identifiable statistical model $\mathcal{M}=\left\{f\right(x;\mathit{\theta }):\mathit{\theta }\in \mathbf{\Theta }\}$, where Θ is the parameter space and $F(x;\mathit{\theta })$ is the associated c.d.f. Then, for ${\mathit{\theta }}_{\mathbf{1}},{\mathit{\theta }}_{\mathbf{2}}\in \mathbf{\Theta }$, the c.d.f. sensitivity between $f(x;{\mathit{\theta }}_{\mathbf{1}})$ and $f(x;{\mathit{\theta }}_{\mathbf{2}})$ with respect to the random sample $\mathit{x}=\{x_1,x_2,\dots ,x_n\}$ is defined as

$$\begin{array}{c}\hfill {\displaystyle \underset{x_i\in \mathit{x}}{max}|F(x_i;{\mathit{\theta }}_{\mathbf{1}})-F(x_i;{\mathit{\theta }}_{\mathbf{2}})|.}\end{array}$$

Definition 5.

Consider an identifiable statistical model $\mathcal{M}=\left\{f\right(x;\mathit{\theta }):\mathit{\theta }\in \mathbf{\Theta }\}$, where Θ is the parameter space and $F(x;\mathit{\theta })$ is the associated c.d.f. Then, for all real numbers $e^{*}>0$, $f(x;{\mathit{\theta }}_{\mathbf{1}})$ and $f(x;{\mathit{\theta }}_{\mathbf{2}})$ are said to be indistinguishable with respect to $e^{*}$ if their c.d.f. sensitivity with respect to the random sample $\mathit{x}=\{x_1,x_2,\dots ,x_n\}$ is no greater than $e^{*}$. That is, for ${\mathit{\theta }}_{\mathbf{1}},{\mathit{\theta }}_{\mathbf{2}}\in \mathbf{\Theta },e^{*}>0$,

$$\begin{array}{c}\hfill {\displaystyle \underset{x_i\in \mathit{x}}{max}|F(x_i;{\mathit{\theta }}_{\mathbf{1}})-F(x_i;{\mathit{\theta }}_{\mathbf{2}})|\le e^{*}.}\end{array}$$

Definition 6.

For a given statistical model $\mathcal{M}=\left\{f\right(x;\mathit{\theta }):\mathit{\theta }\in \mathbf{\Theta }\}$, consider the procedure utilized for obtaining parameter estimates for $\mathit{\theta }$ that are based on the random sample $\mathit{x}=\{x_1,x_2,\dots ,x_n\}$ and the numerical algorithm being implemented. Such a procedure is referred to as experiment Φ.

Definition 7.

Consider an identifiable statistical model $\mathcal{M}=\left\{f\right(x;\mathit{\theta }):\mathit{\theta }\in \mathbf{\Theta }\}$, where Θ is the parameter space and $F(x;\mathit{\theta })$ is the associated c.d.f. Let $f(x;\widehat{\mathit{\theta }})$ be an estimated model of f with respect to experiment Φ and ${\widehat{F}}_n\left(t\right)$ be the empirical c.d.f. (ECDF) obtained from Φ. Then, the experimental error associated with Φ is defined as the c.d.f. sensitivity between the estimated model and the ECDF. That is,

$$\begin{array}{c}\hfill {\displaystyle ϵ:=\underset{x_i\in \mathit{x}}{max}\left|F(x_i;\widehat{\mathit{\theta }})-{\widehat{F}}_n\left(x_i\right)\right|.}\end{array}$$

(5)

Definition 8.

Consider an identifiable statistical model $\mathcal{M}=\left\{f\right(x;\mathit{\theta }):\mathit{\theta }\in \mathbf{\Theta }\}$, where Θ is the parameter space and $F(x;\mathit{\theta })$ is the associated c.d.f. Let $f(x;\widehat{\mathit{\theta }})$ be an estimated model of f with respect to experiment Φ. Then, $f(x;\mathit{\theta })$ is said to be indistinguishable with respect to Φ if $f(x;\mathit{\theta })$ and $f(x;\widehat{\mathit{\theta }})$ are indistinguishable with respect to the experimental error ϵ as defined in (5). That is, for $\mathit{\theta }\in \mathbf{\Theta }$,

$$\begin{array}{c}\hfill {\displaystyle \underset{x_i\in \mathit{x}}{max}\left|F(x_i;\widehat{\mathit{\theta }})-F(x_i;\mathit{\theta })\right|\le ϵ.}\end{array}$$

Definition 9.

Consider an identifiable statistical model $\mathcal{M}=\left\{f\right(x;\mathit{\theta }):\mathit{\theta }\in \mathbf{\Theta }\}$, where Θ is the parameter space and $F(x;\mathit{\theta })$ is the associated c.d.f. Let $f(x;\widehat{\mathit{\theta }})$ be an estimated model of f with respect to experiment Φ; then the set $\mathcal{N}(f,\mathsf{\Phi })\subset \mathbf{\Theta }$ is called a c.d.f. sensitivity-based confidence region (CSCR) with respect to Φ if

$$\begin{array}{c}\hfill \mathcal{N}:=\left\{{\displaystyle \mathit{\theta }\in \mathbf{\Theta }\left|\underset{x_i\in \mathit{x}}{max}\right|F(x_i;\widehat{\mathit{\theta }})-F(x_i;\mathit{\theta })|\le ϵ}\right\},\end{array}$$

(6)

where ϵ is as defined in (5).

Definition 10.

Consider an identifiable statistical model $\mathcal{M}=\left\{f\right(x;\mathit{\theta }):\mathit{\theta }\in \mathbf{\Theta }\}$, where Θ is the parameter space with $dim\left(\mathbf{\Theta }\right)=d_1$. Given $k\in {\mathbb{Z}}^{+}$, let ${\mathcal{P}}_k$ be a sub-space of Θ with $dim\left({\mathcal{P}}_k\right)=d_2<d_1$, and define ${\mathcal{R}}_k:=\mathbf{\Theta }\setminus {\mathcal{P}}_k$ so that $dim\left({\mathcal{R}}_k\right)=d_1-d_2$. Let the parameters in ${\mathcal{R}}_k$ be ${\mathit{r}}_k=\left\{r_k^{\left(1\right)},r_k^{\left(2\right)},\dots ,r_k^{(d_1-d_2)}\right\}$, and let the parameters in ${\mathcal{P}}_k$ be ${\mathit{p}}_k=\left\{p_k^{\left(1\right)},p_k^{\left(2\right)},\dots ,p_k^{\left(d_2\right)}\right\}$. Let ${\mathit{\theta }}_B$ denote the boundary of the domain in the parameter space Θ. Then, a statistical model ${\mathcal{M}}_k:=\left\{g_k\left(x;{\mathit{p}}_k\right):{\mathit{p}}_k\in {\mathcal{P}}_k\right\}$ is said to be the $k^{th}$ sub-model of $\mathcal{M}$ if

$$\begin{array}{c}\hfill \underset{{\mathit{r}}_k\to {\mathit{r}}_{B,k}}{lim}F(x;\mathit{\theta })=G_k\left(x;{\mathit{p}}_k\right),\end{array}$$

where ${\mathit{r}}_{B,k}$ is the boundary of ${\mathcal{M}}_k$ and $G_k$ is the associated c.d.f.

Definition 11.

Consider an identifiable statistical model $\mathcal{M}=\left\{f\right(x;\mathit{\theta }):\mathit{\theta }\in \mathbf{\Theta }\}$, where Θ is the parameter space. Then, $\mathcal{F}$ is said to be the sub-model family of $\mathcal{M}$ if it comprises all the sub-models of $\mathcal{M}$. Namely,

$$\begin{array}{c}\hfill \mathcal{F}:=\bigcup _{k=1}^K{\mathcal{M}}_k,\end{array}$$

where ${\mathcal{M}}_k$ is a sub-model of $\mathcal{M}$ and $K\ge 1$ is the number of sub-models.

This paper’s principal contribution is the definition of the estimability of a statistical model that follows.

Definition 12.

Consider an identifiable statistical model $\mathcal{M}=\left\{f\right(x;\mathit{\theta }):\mathit{\theta }\in \mathbf{\Theta }\}$ with parameter space Θ. Then, $\mathcal{M}$ is said to be non-estimable with respect to experiment Φ if its CSCR with respect to experiment Φ is infinite. Accordingly, $\mathcal{M}$ is said to be estimable if its CSCR with respect to experiment Φ is bounded.

Observe that expression (5) quantifies the experimental error as $ϵ$. It can also be interpreted as the tolerance level within which the estimated model c.d.f. can vary. The CSCR essentially includes all possible parameters such that the c.d.f. sensitivity of a model having those parameters is less than the experimental error. Clearly, the smaller the experimental error, the smaller the CSCR, this being due to the fact that the experimental error is set as an upper bound in (6). The next step is to make Definition 12 applicable in practice, as this definition may not be of practical use if utilized directly. This caveat is addressed in the next two theorems.

Theorem 1.

Consider an identifiable statistical model $\mathcal{M}=\left\{f\right(x;\mathit{\theta }):\mathit{\theta }\in \mathbf{\Theta }\}$, where Θ is the parameter space and $F(x;\mathit{\theta })$ is the associated c.d.f. Assume that $\mathcal{M}$ has a sub-model family $\mathcal{F}:={\bigcup }_{k=1}^K{\mathcal{M}}_k$. Then, $\mathcal{M}$ is non-estimable if there exists at least one sub-model ${\mathcal{M}}_k=\left\{g_k\left(x;{\mathit{p}}_k\right):{\mathit{p}}_k\in {\mathcal{P}}_k\right\}$ that satisfies both of the following conditions:

(i) 

${\displaystyle \underset{x_i\in \mathit{x}}{max}\left|F(x_i;\widehat{\mathit{\theta }})-G_k\left(x_i;{\widehat{\mathit{p}}}_k\right)\right|\le ϵ}$, where ϵ is as defined in (5), and F and $G_k$ are the associated c.d.f. values.

(ii) 

${\mathit{r}}_{B,k}$ contains ∞ or $-\infty$.

Proof. 

Let a sub-model ${\mathcal{M}}_k$ satisfy both of the conditions, where $k=1,2,\dots ,K$. By Definition 10, we have

$$\begin{array}{c}\hfill \underset{{\mathit{r}}_k\to {\mathit{r}}_{B,k}}{lim}F(x;\mathit{\theta })=G_k\left(x;{\mathit{p}}_k\right).\end{array}$$

Let ${\tilde{\mathit{p}}}_k\left({\mathit{r}}_k\right)$ be the estimates of ${\mathit{p}}_k$ for given values of ${\mathit{r}}_k$. We write it this way to reflect its dependence on ${\mathit{r}}_k$. Assuming that the parameters are estimated under the same numerical algorithm, then for all x, we have

$$\begin{array}{c}\hfill \underset{{\mathit{r}}_k\to {\mathit{r}}_{B,k}}{lim}F\left(x;{\tilde{\mathit{p}}}_k\left({\mathit{r}}_k\right)\cup {\mathit{r}}_k\right)=G_k\left(x;{\widehat{\mathit{p}}}_k\right).\end{array}$$

(7)

Subsequently, based on (7), we have

$$\begin{array}{c}\hfill {\displaystyle \underset{{\mathit{r}}_k\to {\mathit{r}}_{B,k}}{lim}\underset{x_i\in \mathit{x}}{max}\left|F(x_i;\widehat{\mathit{\theta }})-F\left(x_i;{\tilde{\mathit{p}}}_k\left({\mathit{r}}_k\right)\cup {\mathit{r}}_k\right)\right|=\underset{x_i\in {\mathit{x}}_k}{max}\left|F(x_i;{\widehat{\mathit{\theta }}}_k)-G_k\left(x_i;{\widehat{\mathit{p}}}_k\right)\right|.}\end{array}$$

(8)

Denote the $j^{th}$ element in ${\mathit{r}}_k$ by $r_k^{\left(j\right)}$. In light of condition (ii), $r_k^{\left(j\right)}$ then becomes either ∞ or $-\infty$. On applying condition (i), (8) then implies that for some j,

$$\begin{array}{c}\hfill {\displaystyle \exists {\delta }_k^{\left(j\right)}>0\mathrm{such}\mathrm{that}\underset{x_i\in \mathit{x}}{max}\left|F(x_i;\widehat{\mathit{\theta }})-F\left(x_i;{\tilde{\mathit{p}}}_k\left({\mathit{r}}_k\right)\cup {\mathit{r}}_k\right)\right|\le ϵ,\forall \left|r_k^{\left(j\right)}\right|>{\delta }_k^{\left(j\right)}.}\end{array}$$

Thus, ${\tilde{\mathit{p}}}_k\left({\mathit{r}}_k\right)\cup {\mathit{r}}_k$ belongs to the CSCR by Definition 9. In that case, the CSCR is infinite given the parameters included in ${\mathcal{R}}_{\mathit{k}}$. Then, by Definition 12, $\mathcal{M}$ is non-estimable with respect to experiment $\mathsf{\Phi }$.    □

Theorem 2.

Consider an identifiable statistical model $\mathcal{M}=\left\{f\right(x;\mathit{\theta }):\mathit{\theta }\in \mathbf{\Theta }\}$, where Θ is the parameter space and $F(x;\mathit{\theta })$ is the associated c.d.f. Assume that $\mathcal{M}$ has a sub-model family $\mathcal{F}:={\bigcup }_{k=1}^K{\mathcal{M}}_k$. Then, $\mathcal{M}$ is estimable if, for all ${\mathcal{M}}_k\in \mathcal{F}$, conditions (i) and (ii) specified in Theorem 1 cannot be simultaneously satisfied.

Proof. 

If condition (ii) is not satisfied, then ${\mathit{r}}_{B,k}$ does not contain ∞ or $-\infty$. Denote the $j^{th}$ element in ${\mathit{r}}_k$ by $r_k^{\left(j\right)}$, and denote the $j^{th}$ element in ${\mathit{r}}_{B,k}$ by $r_{B,k}^{\left(j\right)}$. Then, (7) implies that for any $e>0$ and for all $j=1,2,\dots ,d_1-d_2$, $\exists {\delta }_k^{\left(j\right)}>0$ such that ${\displaystyle \underset{x_i\in \mathit{x}}{max}\left|F\left(x_i;{\tilde{\mathit{p}}}_k\left({\mathit{r}}_k\right)\cup {\mathit{r}}_k\right)-G_k\left(x;{\widehat{\mathit{p}}}_k\right)\right|\le e}$, $\forall r_k^{\left(j\right)}\in \left(r_{B,k}^{\left(j\right)}-{\delta }_k^{\left(j\right)},r_{B,k}^{\left(j\right)}+{\delta }_k^{\left(j\right)}\right)$.

Therefore, no matter whether or not condition (i) is satisfied, ${\tilde{\mathit{p}}}_k\left({\mathit{r}}_k\right)\cup {\mathit{r}}_k$ is already contained in a bounded region. The CSCR, which is a sub-region of this bounded region, is, of course, bounded as well.

It remains to consider the case when condition (i) is not satisfied while condition (ii) is. In that case, since condition (i) is not verified, we have

$$\begin{array}{c}\hfill {\displaystyle \underset{x_i\in \mathit{x}}{max}\left|F(x_i;\widehat{\mathit{\theta }})-G_k\left(x_i;{\widehat{\mathit{p}}}_k\right)\right|>ϵ.}\end{array}$$

Since condition (ii) is satisfied, (8) then implies that for all j,

$$\begin{array}{c}\hfill {\displaystyle \exists {\delta }_k^{\left(j\right)}>0\mathrm{such}\mathrm{that}\underset{x_i\in \mathit{x}}{max}\left|F(x_i;\widehat{\mathit{\theta }})-F\left(x_i;{\tilde{\mathit{p}}}_k\left({\mathit{r}}_k\right)\cup {\mathit{r}}_k\right)\right|>ϵ,\forall \left|r_k^{\left(j\right)}\right|>{\delta }_k^{\left(j\right)}.}\end{array}$$

This is equivalent to

$$\begin{array}{c}\hfill \exists {\delta }_k^{\left(j\right)}>0\mathrm{such}\mathrm{that}\left\{{\displaystyle |r_k^{\left(j\right)}||\underset{x_i\in \mathit{x}}{max}|F(x_i;\widehat{\mathit{\theta }})-F\left(x_i;{\tilde{\mathit{p}}}_k\left({\mathit{r}}_k\right)\cup {\mathit{r}}_k\right)|\le ϵ}\right\}\subset \left(0,{\delta }_k^{\left(j\right)}\right].\end{array}$$

Therefore, for each parameter of ${\tilde{\mathit{p}}}_k\left({\mathit{r}}_k\right)\cup {\mathit{r}}_k$ that is in the CSCR by Definition 9, its range in absolute value is covered by a bounded interval $\left(0,{\delta }_k^{\left(j\right)}\right]$. We may then conclude that the entire CSCR is a bounded region. Thus, by Definition 12, $\mathcal{M}$ is estimable with respect to experiment $\Phi$.    □

Assuming that the statistical model under consideration is identifiable, the steps to follow for implementing the proposed methodology are enumerated in the next algorithm.

Algorithm 1 The estimability analysis for a statistical model $\mathcal{M}=\left\{f\right(x;\mathit{\theta }):\mathit{\theta }\in \mathbf{\Theta }\}$

1: Determine the sub-models of $\mathcal{M}$—a requirement to conduct the proposed estimability analysis. 2: Obtain a parameter estimate of $\mathcal{M}$, $\widehat{\mathit{\theta }}$. The estimator should be consistent. It is recommended to proceed as explained in Section 7.1. 3: Compute the experimental error $ϵ$, as specified in Definition 7. 4: Proceeding as in Step 2, obtain the parameter estimates for all the sub-models of the statistical model. 5: Compute the c.d.f. sensitivities between the estimated statistical model and each one of its estimated sub-models using Definition 4. 6: Compare all the c.d.f. sensitivities obtained in Step 5 with the experimental error $ϵ$ determined in Step 3. If all c.d.f. sensitivities are greater than the $ϵ$, then the statistical model is estimable with respect to the current experiment $\mathsf{\Phi }$ in light of Theorem 2. Otherwise, it is non-estimable with respect to the current experiment $\mathsf{\Phi }$ by virtue of Theorem 1. 7: If non-estimable, one may consider improving the experimental design by making use of the techniques described in Section 5.3 and then reassess estimability, starting with Step 2.

5. Validation of the Proposed Definition

In this section, we substantiate the proposed definition and conclude that it is innately sound. This is achieved by establishing the validity of the associated theoretical results, which happen to be consistent with common sense.

We are considering multiple perspectives, as discussed in Section 5.1, Section 5.2 and Section 5.3. This coherence is further supported by the illustrative examples presented in Section 5.4.

5.1. Validation of the Data Noise, the Algorithm Noise, and the Experimental Error

We first validate the data noise, the algorithm noise, and the experimental error as specified in Definitions 13 and 7. Without any loss of generality, it is assumed that the experimental error originates from the data noise and the algorithm noise (Chis et al. 2011). The data noise and the algorithm noise are next defined as c.d.f. sensitivity measures:

Definition 13.

Consider an identifiable statistical model $\mathcal{M}=\left\{f\right(x;\mathit{\theta }):\mathit{\theta }\in \mathbf{\Theta }\}$, where Θ is the parameter space and $F(x;\mathit{\theta })$ is the associated c.d.f. Let $f(x;\widehat{\mathit{\theta }})$ be an estimated model with respect to experiment Φ. Denoting the true value of the parameter estimate by ${\widehat{\mathit{\theta }}}^T$, the data noise ${ϵ}_d$ and the algorithm noise ${ϵ}_{AL}$ of $\mathcal{M}$ with respect to experiment Φ can be defined as the following c.d.f. sensitivities, respectively:

$$\begin{array}{cc}\hfill {ϵ}_d& {\displaystyle :=\underset{x_i\in \mathit{x}}{max}\left|F_X({\widehat{\mathit{\theta }}}^T;x_i)-{\widehat{F}}_n\left(x_i\right)\right|,}\hfill \\ \hfill {ϵ}_{AL}& {\displaystyle :=\underset{x_i\in \mathit{x}}{max}\left|F_X({\widehat{\mathit{\theta }}}^T;x_i)-F_X(\widehat{\mathit{\theta }};x_i)\right|.}\hfill \end{array}$$

We now show that both types of noise tend to zero as the sample size increases. As the sample size goes to infinity,

$$\begin{array}{cc}\hfill {ϵ}_d& {\displaystyle :=\underset{x_i\in \mathit{x}}{max}\left|F_X({\widehat{\mathit{\theta }}}^T;x_i)-{\widehat{F}}_n\left(x_i\right)\right|}\hfill \\ \hfill & {\displaystyle \to \underset{x\in \mathbf{\Omega }}{max}|F_X({\mathit{\theta }}^{*};x)-F_X({\mathit{\theta }}^{*};x)|}\hfill \\ \hfill & =0,\hfill \end{array}$$

where ${\mathit{\theta }}^{*}$ is the true parameter of the model (unknown, of course) and $\mathbf{\Omega }$ is the support of $F_X$. Note that the limits of $F_X({\widehat{\mathit{\theta }}}^T;x_i)$ and ${\widehat{F}}_n\left(x_i\right)$ are based on asymptotic results on the consistency of the estimator ${\widehat{\mathit{\theta }}}^T$ and the convergence of the ECDF ${\widehat{F}}_n\left(t\right)$, respectively. Thus, ${ϵ}_d$ is a valid measure of the data noise that behaves as expected.

Now, consider the algorithm noise. As the accuracy of the algorithm tends to perfection,

$$\begin{array}{cc}\hfill {ϵ}_{AL}& {\displaystyle :=\underset{x_i\in \mathit{x}}{max}\left|F_X({\widehat{\mathit{\theta }}}^T;x_i)-F_X(\widehat{\mathit{\theta }};x_i)\right|}\hfill \\ \hfill & {\displaystyle \to \underset{x_i\in \mathit{x}}{max}\left|F_X({\widehat{\mathit{\theta }}}^T;x_i)-F_X({\widehat{\mathit{\theta }}}^T;x_i)\right|}\hfill \\ \hfill & =0.\hfill \end{array}$$

Thus, as the algorithm is becoming nearly exact, the approximation of ${\widehat{\mathit{\theta }}}^T$ is more accurate, which also agrees with common sense. Accordingly, ${ϵ}_{AL}$ is a valid measure of the algorithm noise as well.

After having validated the data noise and the algorithm noise, we may now validate the experimental error. Observe that

$$\begin{array}{cc}\hfill ϵ:=& {\displaystyle \underset{x_i\in \mathit{x}}{max}\left|F_X(\widehat{\mathit{\theta }};x_i)-{\widehat{F}}_n\left(x_i\right)\right|}\hfill \\ \hfill =& {\displaystyle \underset{x_i\in \mathit{x}}{max}\left|F_X(\widehat{\mathit{\theta }};x_i)-F_X({\widehat{\mathit{\theta }}}^T;x_i)+F_X({\widehat{\mathit{\theta }}}^T;x_i)-{\widehat{F}}_n\left(x_i\right)\right|}\hfill \\ \hfill \le & {\displaystyle \underset{x_i\in \mathit{x}}{max}\left|F_X({\widehat{\mathit{\theta }}}^T;x_i)-{\widehat{F}}_n\left(x_i\right)\right|+\underset{x_i\in \mathit{x}}{max}\left|F_X(\widehat{\mathit{\theta }};x_i)-F_X({\widehat{\mathit{\theta }}}^T;x_i)\right|}\hfill \\ \hfill =& {ϵ}_d+{ϵ}_{AL}.\hfill \end{array}$$

(9)

Then, as the experimental design tends to perfection, both the data noise and the algorithm noise tend to zero. Consequently, $ϵ$ tends to zero by appealing to the so-called “Squeeze Theorem”. It is worth emphasizing that ${ϵ}_d$ and ${ϵ}_{AL}$ are unknown to the experimenter since the true value of the parameter estimate ${\widehat{\mathit{\theta }}}^T$ is also unknown. Instead, one only knows $\widehat{\mathit{\theta }}$, which is the output of the numerical algorithm. However, this does not prevent us from establishing that $ϵ$ is a valid measure of the experimental error. Otherwise, one would be able to calculate ${\widehat{\mathit{\theta }}}^T$ analytically without resorting to numerical algorithms.

5.2. Validation of the c.d.f. Sensitivity-Based Confidence Region

We next validate the definition of the CSCR. Observe that the experimental error, $ϵ$, is an upper bound in (6). Thus, as the experimental error decreases, the CSCR shrinks, making the model more estimable, which is consistent with intuition. Now, as in Section 5.1, consider the argument that the experimental design tends to perfection. Then,

$$\begin{array}{cc}\hfill \mathcal{N}& :=\left\{{\displaystyle \mathit{\theta }\in \mathbf{\Theta }\left|\underset{x_i\in \mathit{x}}{max}\right|F(x_i;\widehat{\mathit{\theta }})-F(x_i;\mathit{\theta })|\le ϵ}\right\}\hfill \\ \hfill & \to \left\{{\displaystyle \mathit{\theta }\in \mathbf{\Theta }|\underset{x\in \mathbf{\Omega }}{max}\left|F(x;{\mathit{\theta }}^{*})-F(x;\mathit{\theta })\right|\le 0}\right\}\hfill \\ \hfill & =\left\{{\displaystyle \mathit{\theta }\in \mathbf{\Theta }|\underset{x\in \mathbf{\Omega }}{max}\left|F(x;{\mathit{\theta }}^{*})-F(x;\mathit{\theta })\right|=0}\right\}\hfill \\ \hfill & =\left\{\mathit{\theta }\in \mathbf{\Theta }|F(x;{\mathit{\theta }}^{*})-F(x;\mathit{\theta })=0,\forall x\in \mathbf{\Omega }\right\}\hfill \\ \hfill & =:\mathcal{A},\hfill \end{array}$$

where ${\mathit{\theta }}^{*}$ is the true (unknown) parameter of the model.

Thus, the CSCR, $\mathcal{N}$, collapses to $\mathcal{A}$, which is exactly the parameter set for which $\mathcal{M}$ is non-identifiable. Therefore, if the experimental error goes to zero, then any remaining unreliability associated with the parameter estimates must originate from the non-identifiability issue. This agrees with the statement found in the literature to the effect that it is unnecessary to analyze estimability if the model is non-identifiable since the non-identifiability issue cannot be overcome by improvements in the experimental design. Accordingly, the demonstrated consistency with respect to the CSCR further validates the proposed definition of estimability.

5.3. Validation by Known Techniques

The proposed definition can also be substantiated by making use of known methods for improving estimability. These are discussed next.

5.3.1. Increasing the Sample Size

Increasing the sample size can decrease both ${ϵ}_d$ and ${ϵ}_{AL}$. More observations not only decrease the data noise but also decrease the algorithm noise as the likelihood functions then become more and more concave. Then, in light of the arguments presented in Section 5.1 and Section 5.2, $ϵ$ decreases, and the CSCR shrinks, which improves estimability. As the sample size tends to infinity, $ϵ$ tends to zero, and the experiment tends to perfection. Thus, the proposed definition is consistent with an increase in the sample size.

5.3.2. Increasing the Convexity of the Log-Likelihood Function

Another way to increase the sample size is to clone the simulated data multiple times. This is called the data cloning method (Lele et al. 2007, 2010; Cheng 2023a). With this approach, the likelihood functions become more and more concave, which makes the algorithm approximation more accurate. Accordingly, ${ϵ}_{AL}$ decreases. However, this approach only reduces the algorithm noise. It cannot improve ${ϵ}_d$, as ${\widehat{\mathit{\theta }}}^T$ does not change. Thus, the proposed definition is consistent with the data cloning method.

5.3.3. Improving the Algorithm Design

One can also decrease ${ϵ}_{AL}$ by improving the efficiency of the algorithm approximation. However, the data noise remains unchanged, which is consistent with intuition as the improvement only pertains to the algorithm.

5.3.4. Securing More Complete Information

Estimability can also be enhanced if more information is secured. As the data convey additional information, ${ϵ}_d$ decreases. Then, in light of the arguments presented in Section 5.1 and Section 5.2, $ϵ$ decreases, and the CSCR shrinks, which improves estimability. Thus, the proposed definition is consistent with obtaining more complete information.

5.4. Illustrative Examples

In Section 5.1, Section 5.2 and Section 5.3, we have theoretically validated our definition of estimability by establishing that it is consistent with common sense. Examples illustrating its applicability to discrete and continuous statistical models follow.

5.4.1. Example 1

Consider the following constrained binomial model:

$$\begin{array}{c}\hfill \mathcal{M}=\left\{f\left(k;m,p=\frac{\lambda }{m}\right)=\left(\genfrac{}{}{0pt}{}{m}{k}\right){\left(\frac{\lambda }{m}\right)}^k{\left(1-\frac{\lambda }{m}\right)}^{m-k}:m>\lambda >0,m\in {\mathbb{Z}}^{+}\right\}.\end{array}$$

With the binomial distribution being identifiable, its estimability may be assessed. According to Definition 11, $\mathbf{\Theta }=\{m,\lambda \}$, ${\mathcal{P}}_1=\left\{\lambda \right\}$, ${\mathcal{R}}_1=\left\{m\right\}$, ${\mathit{\theta }}_B=\{0,\infty \}$, $r_{B,1}=\{\infty \}$. Then, the sub-model family of $\mathcal{M}$ based on $r_{B,1}$ is

$$\begin{array}{c}\hfill \mathcal{F}={\mathcal{M}}_1:=\left\{g_1(k;\lambda )=\frac{e^{-\lambda }{\lambda }^k}{k!}:\lambda >0\right\}\end{array}$$

since the Poisson distribution is the limiting distribution of the binomial distribution as $m\to \infty$ with the constraint $mp=\lambda$.

Let the underlying model have as its parameters $m=10$ and $\lambda =0.7$ (or $p=0.07$), and consider following three experiments:

• Experiment 1: The sample size is $\left|\mathit{k}\right|=50$.
• Experiment 2: The sample size is $\left|\mathit{k}\right|=500$.
• Experiment 3: The sample size is $\left|\mathit{k}\right|=1000$.

The estimability assessment results are presented in

Table 1. c.d.f. sensitivity and experimental error for the constrained binomial model—Experiments 1 to 3. The asterisk indicates that the binomial model is non-estimable.

Sample Size	Sub-Model	Sub-Model Boundary	${\displaystyle \underset{{\mathit{k}}_{\mathit{i}}\in \mathit{k}}{max}\left\{\left|\mathit{F}({\mathit{k}}_{\mathit{i}};\widehat{\mathit{m}},\widehat{\mathit{\lambda }})-{\mathit{G}}_1({\mathit{k}}_{\mathit{i}};\widehat{\mathit{\lambda }})\right|\right\}}$	$\mathit{ϵ}$
50	${\mathcal{M}}_1$	$r_{B,1}$	0.00426 *	$0.05133$
500	${\mathcal{M}}_1$	$r_{B,1}$	0.00449 *	$0.02269$
1000	${\mathcal{M}}_1$	$r_{B,1}$	0.01276	$0.00480$

and Figure 5 and Figure 6. Table 1 compares the experimental error with the c.d.f. sensitivity between the fitted binomial distribution and its sub-model, which is a direct application of Theorems 1 and 2.

It can be observed that the experimental error decreases as the sample size increases, which again validates the definition of the experimental error. According to Theorem 1, $\mathcal{M}$ is non-estimable with respect to Experiments 1 and 2 since the c.d.f. sensitivity measure is less than the experimental error. On the other hand, according to Theorem 2, $\mathcal{M}$ becomes estimable with respect to Experiment 3 since the c.d.f. sensitivity measure is greater.

The above conclusions are further supported by Figure 5, where the CSCRs for Experiments 1, 2, and 3 are visualized. In line with Theorem 1, the CSCR for Experiments 1 and 2 are both infinite. As the sample size increases, the CSCR shrinks until it becomes bounded in Experiment 3, where $\mathcal{M}$ becomes estimable, which is in line with Theorem 2.

A more intuitive way of interpreting the concept of estimability is displayed in Figure 6. In Experiments 1 and 2, the model is assessed to be non-estimable. This can be interpreted as follows: the inferential power from the data displayed in the histogram is not sufficient to distinguish the estimated p.m.f.s of the underlying model and its sub-model. This is exactly why estimability is also referred to as “practical identifiability” since the shape of the histogram cannot “practically identify” the estimated p.m.f. values of the underlying model and its sub-model. On the other hand, in Experiment 3, the model is assessed to be estimable. Thus, the inferential power is sufficient to distinguish the estimated p.m.f. values. In this case, the shape of the histogram is sufficient to favor the underlying model and dismiss its sub-model.

However, these conclusions cannot be reached by only inspecting Figure 6, which is why the CSCR associated with the proposed definition is crucial.

5.4.2. Example 2

Consider the following constrained Lomax Pareto Type II (henceforth Pareto) model:

$$\begin{array}{c}\hfill \mathcal{M}=\left\{f\left(x;\alpha ,\theta =\frac{\alpha }{\lambda }\right)=\frac{\alpha {\left(\frac{\alpha }{\lambda }\right)}^{\alpha }}{{\left(x+\frac{\alpha }{\lambda }\right)}^{\alpha +1}}:\lambda >0,\alpha >0\right\}.\end{array}$$

With the Pareto distribution being identifiable, its estimability may be assessed. According to Definition 11, $\mathbf{\Theta }=\{\alpha ,\lambda \}$, ${\mathcal{P}}_1=\left\{\lambda \right\}$, ${\mathcal{R}}_1=\left\{\alpha \right\}$, ${\mathit{\theta }}_B=\{0,\infty \}$, $r_{B,1}=\{\infty \}$. Then, the sub-model family of $\mathcal{M}$ based on $r_B$ is

$$\begin{array}{c}\hfill \mathcal{F}={\mathcal{M}}_1:=\left\{g_1(x;\lambda )=\lambda e^{-\lambda x}:\lambda >0\right\}\end{array}$$

since the exponential distribution is the limiting distribution of the Pareto distribution as $\alpha \to \infty$ with the constraint $\frac{\alpha }{\theta }=\lambda$.

Let the underlying model have as its parameters $\alpha =3$ and $\lambda =0.1$ (or $\theta =30$), and consider following three experiments:

• Experiment 1: The sample size is $\left|\mathit{x}\right|=100$.
• Experiment 2: The sample size is $\left|\mathit{x}\right|=200$.
• Experiment 3: The sample size is $\left|\mathit{x}\right|=500$.

The estimability assessment results are presented in

Table 2. c.d.f. sensitivity and experimental error for the constrained Pareto model—Experiments 1 to 3. The asterisk indicates that the Pareto model is non-estimable.

Sample Size	Sub-Model	Sub-Model Boundary	${\displaystyle \underset{{\mathit{x}}_{\mathit{i}}\in \mathit{x}}{max}\left\{\left|\mathit{F}({\mathit{x}}_{\mathit{i}};\widehat{\mathit{\alpha }},\widehat{\mathit{\lambda }})-{\mathit{G}}_1({\mathit{x}}_{\mathit{i}};\widehat{\mathit{\lambda }})\right|\right\}}$	$\mathit{ϵ}$
100	${\mathcal{M}}_1$	$r_{B,1}$	0.02081 *	$0.05438$
200	${\mathcal{M}}_1$	$r_{B,1}$	0.02628 *	$0.02685$
500	${\mathcal{M}}_1$	$r_{B,1}$	$0.03654$	$0.02500$

and Figure 7 and Figure 8. Table 2 compares the experimental error with the c.d.f. sensitivity between the fitted Pareto distribution and its sub-model, which is a direct application of Theorems 1 and 2. The results yield the same conclusions that were obtained in Example 1.

It should be emphasized that the visualization of the CSCR displayed in Figure 5 and Figure 7 is not achievable if the parameter space extends to more than three dimensions, such as in the case of the PTAM, which is considered next. In that case, we have to fully rely on Theorems 1 and 2 and compare the c.d.f. sensitivity with the experimental error. The two-dimensional illustrative examples that were previously presented provide supporting evidence corroborating the validity of Theorems 1 and 2. Another important aspect to point out is that the algorithm utilized to obtain the parameter estimates coincides with that recommended in Section 7.1. This aspect is discussed further in the remainder of this paper.

6. Estimability of the PTAM

In Section 5, we establish that the proposed definition of estimability is innately sound. In this section, we apply this definition to assess the estimability of the PTAM.

6.1. Identifiability of the PTAM

Nie (2022) proved that the PTAM is identifiable when the number of states is greater or equal to six. We may thus proceed to assess its estimability when that is the case.

6.2. Sub-Models of the PTAM

According to Definition 3, it can be concluded that the parameter space for the PTAM, $\mathbf{\Theta }$, is as follows:

$$\begin{array}{cc}\hfill \mathbf{\Theta }& =\left\{\mathit{\theta }=\left(h_1,h_m,s,\psi ,m\right)|h_m>h_1>0,\frac{m}{\psi }=\lambda ,s\in \mathbb{R},\psi >0,m\ge 6\right\},\hfill \end{array}$$

with boundary ${\mathit{\theta }}_B=\{0,6,h_1,h_m,\infty ,-\infty \}$. Next, one must obtain sub-models of the PTAM in order to investigate its parameter estimability. According to Theorems 1 and 2, it suffices to investigate sub-models whose boundaries contain ∞ or $-\infty$. The results are presented in the following propositions.

Proposition 1.

Consider the PTAM as specified in Definition 3. Given m, the limiting distribution as $s\to \infty$ is Coxian of order 2 with

$$\begin{array}{c}\hfill \mathit{S}=\left[\begin{array}{cc}-(\lambda +h_1)& \lambda \\ 0& -h_m\end{array}\right],\mathit{h}=\left[\begin{array}{c}{h}_1\\ h_m\end{array}\right].\end{array}$$

Denote this as sub-model ${\mathcal{M}}_1$ with corresponding boundary $r_{B,1}$. In this case, ${\mathcal{R}}_1=\left\{s\right\}$, ${\mathcal{P}}_1=\{h_1,h_m,\psi ,m\}$, and $r_{B,1}=\{\infty \}$.

Proposition 2.

Consider the PTAM as specified in Definition 3. Given m, the limiting distribution as $h_1\to 0$ and $s\to \infty$ is Coxian of order 2 with

$$\begin{array}{c}\hfill \mathit{S}=\left[\begin{array}{cc}-\lambda & \lambda \\ 0& -h_m\end{array}\right],\mathit{h}=\left[\begin{array}{c}0\\ h_m\end{array}\right].\end{array}$$

Denote this as sub-model ${\mathcal{M}}_2$ with corresponding boundary ${\mathit{r}}_{B,2}$. In this case, ${\mathcal{R}}_2=\{h_1,s\}$, ${\mathcal{P}}_2=\{h_m,\psi ,m\}$, and ${\mathit{r}}_{B,2}=\{0,\infty \}$.

Proposition 3.

Consider the PTAM as specified in Definition 3. Given m, the limiting distribution as $s\to -\infty$ is Coxian of order m with

$$\begin{array}{c}\hfill \mathit{S}=\left[\begin{array}{ccccccc}-(\lambda +h_1)& \lambda & 0& 0& \cdots & 0& 0\\ 0& -(\lambda +h_1)& \lambda & 0& \cdots & 0& 0\\ ⋮& ⋮& ⋮& ⋮& \ddots & ⋮& ⋮\\ 0& 0& 0& 0& \cdots & -(\lambda +h_1)& \lambda \\ 0& 0& 0& 0& \cdots & 0& -h_m\end{array}\right],\mathit{h}=\left[\begin{array}{c}{h}_1\\ h_1\\ ⋮\\ h_1\\ h_m\end{array}\right].\end{array}$$

Denote this as sub-model ${\mathcal{M}}_3$ with corresponding boundary $r_{B,3}$. In this case, ${\mathcal{R}}_3=\left\{s\right\}$, ${\mathcal{P}}_3=\{h_1,h_m,\psi ,m\}$, and $r_{B,3}=\{-\infty \}$.

Proposition 4.

Consider the PTAM as specified in Definition 3. Given m, the limiting distribution as $s\to -\infty$ and $m\to 6$ is Coxian of order 6 with

$$\begin{array}{c}\hfill \mathit{S}=\left[\begin{array}{ccccccc}-(\lambda +h_1)& \lambda & 0& 0& \cdots & 0& 0\\ 0& -(\lambda +h_1)& \lambda & 0& \cdots & 0& 0\\ ⋮& ⋮& ⋮& ⋮& \ddots & ⋮& ⋮\\ 0& 0& 0& 0& \cdots & -(\lambda +h_1)& \lambda \\ 0& 0& 0& 0& \cdots & 0& -h_m\end{array}\right],\mathit{h}=\left[\begin{array}{c}{h}_1\\ h_1\\ ⋮\\ h_1\\ h_m\end{array}\right].\end{array}$$

Denote this as sub-model ${\mathcal{M}}_4$ with corresponding boundary $r_{B,4}$. In that case, ${\mathcal{R}}_4=\{s,m\}$, ${\mathcal{P}}_4=\{h_1,h_m,\psi \}$ and ${\mathit{r}}_{B,4}=\{-\infty ,6\}$.

Proposition 5.

Consider the PTAM as specified in Definition 3. Given m, the limiting distribution as $h_1\to 0$ and $s\to -\infty$ is Coxian of order m with

$$\begin{array}{c}\hfill \mathit{S}=\left[\begin{array}{ccccccc}-\lambda & \lambda & 0& 0& \cdots & 0& 0\\ 0& -\lambda & \lambda & 0& \cdots & 0& 0\\ ⋮& ⋮& ⋮& ⋮& \ddots & ⋮& ⋮\\ 0& 0& 0& 0& \cdots & -\lambda & \lambda \\ 0& 0& 0& 0& \cdots & 0& -h_m\end{array}\right],\mathit{h}=\left[\begin{array}{c}0\\ 0\\ ⋮\\ 0\\ h_m\end{array}\right].\end{array}$$

Denote this as sub-model ${\mathcal{M}}_5$ with corresponding boundary ${\mathit{r}}_{B,5}$. In that case, ${\mathcal{R}}_5=\{h_1,s\}$, ${\mathcal{P}}_5=\{h_m,\psi ,m\}$ and ${\mathit{r}}_{B,5}=\{0,-\infty \}$.

Proposition 6.

Consider the PTAM as specified in Definition 3. Given m, the limiting distribution as $h_1\to 0$, $s\to -\infty$, and $m\to 6$ is Coxian of order 6 with

$$\begin{array}{c}\hfill \mathit{S}=\left[\begin{array}{ccccccc}-\lambda & \lambda & 0& 0& \cdots & 0& 0\\ 0& -\lambda & \lambda & 0& \cdots & 0& 0\\ ⋮& ⋮& ⋮& ⋮& \ddots & ⋮& ⋮\\ 0& 0& 0& 0& \cdots & -\lambda & \lambda \\ 0& 0& 0& 0& \cdots & 0& -h_m\end{array}\right],\mathit{h}=\left[\begin{array}{c}0\\ 0\\ ⋮\\ 0\\ h_m\end{array}\right].\end{array}$$

Denote this as sub-model ${\mathcal{M}}_6$ with corresponding boundary ${\mathit{r}}_{B,6}$. In this case, ${\mathcal{R}}_6=\{h_1,s,m\}$, ${\mathcal{P}}_6=\{h_m,\psi \}$, and ${\mathit{r}}_{B,6}=\{0,-\infty ,6\}$.

Proposition 7.

Consider the PTAM as specified in Definition 3. The limiting distribution under either one of the following cases is exponential with rate $h_1+\lambda$:

(i) 

$h_m\to \infty$ and $s\to \infty$, in which case ${\mathcal{R}}_7=\{h_m,s\}$, ${\mathcal{P}}_7=\{h_1,\psi ,m\}$, and ${\mathit{r}}_{B,7}=\{\infty ,\infty \}$.

(ii) 

$h_m\to \infty$, $s\to \infty$, and $m\to 6$, in which case ${\mathcal{R}}_7=\{h_m,s,m\}$, ${\mathcal{P}}_7=\{h_1,\psi \}$, and ${\mathit{r}}_{B,7}=\{\infty ,\infty ,6\}$.

Denote this as sub-model ${\mathcal{M}}_7$ with corresponding boundary ${\mathit{r}}_{B,7}$.

Proposition 8.

Consider the PTAM as specified in Definition 3. The limiting distribution under either one of the following cases is exponential with rate λ:

(i) 

$h_1\to 0$, $h_m\to \infty$, and $s\to \infty$, in which case ${\mathcal{R}}_8=\{h_1,h_m,s\}$, ${\mathcal{P}}_8=\{\psi ,m\}$, and ${\mathit{r}}_{B,8}=\{0,\infty ,\infty \}$.

(ii) 

$h_1\to 0$, $h_m\to \infty$, $s\to \infty$, and $m\to 6$, in which case ${\mathcal{R}}_8=\{h_1,h_m,s,m\}$, ${\mathcal{P}}_8=\left\{\psi \right\}$, and ${\mathit{r}}_{B,8}=\{0,\infty ,\infty ,6\}$.

Denote this as sub-model ${\mathcal{M}}_8$ with corresponding boundary ${\mathit{r}}_{B,8}$.

Proposition 9.

Consider the PTAM as specified in Definition 3. Given m, the limiting distribution as $h_m\to \infty$ and $s\to -\infty$ is Coxian of order $m-1$ with

$$\begin{array}{c}\hfill \mathit{S}=\left[\begin{array}{ccccccc}-(\lambda +h_1)& \lambda & 0& 0& \cdots & 0& 0\\ 0& -(\lambda +h_1)& \lambda & 0& \cdots & 0& 0\\ ⋮& ⋮& ⋮& ⋮& \ddots & ⋮& ⋮\\ 0& 0& 0& 0& \cdots & -(\lambda +h_1)& \lambda \\ 0& 0& 0& 0& \cdots & 0& -(\lambda +h_1)\end{array}\right],\mathit{h}=\left[\begin{array}{c}{h}_1\\ h_1\\ ⋮\\ h_1\\ h_1+\lambda \end{array}\right].\end{array}$$

Denote this as sub-model ${\mathcal{M}}_9$ with corresponding boundary ${\mathit{r}}_{B,9}$. In this case, ${\mathcal{R}}_9=\{h_m,s\}$, ${\mathcal{P}}_9=\{h_1,\psi ,m\}$, and ${\mathit{r}}_{B,9}=\{\infty ,-\infty \}$.

Proposition 10.

Consider the PTAM as specified in Definition 3. Given m, the limiting distribution as $h_m\to \infty$, $s\to -\infty$, and $m\to 6$ is Coxian of order 5 with

$$\begin{array}{c}\hfill \mathit{S}=\left[\begin{array}{ccccccc}-(\lambda +h_1)& \lambda & 0& 0& \cdots & 0& 0\\ 0& -(\lambda +h_1)& \lambda & 0& \cdots & 0& 0\\ ⋮& ⋮& ⋮& ⋮& \ddots & ⋮& ⋮\\ 0& 0& 0& 0& \cdots & -(\lambda +h_1)& \lambda \\ 0& 0& 0& 0& \cdots & 0& -(\lambda +h_1)\end{array}\right],\mathit{h}=\left[\begin{array}{c}{h}_1\\ h_1\\ ⋮\\ h_1\\ h_1+\lambda \end{array}\right].\end{array}$$

Denote this as sub-model ${\mathcal{M}}_{10}$ with corresponding boundary ${\mathit{r}}_{B,10}$. In this case, ${\mathcal{R}}_{10}=\{h_m,s,m\}$, ${\mathcal{P}}_{10}=\{h_1,\psi \}$, and ${\mathit{r}}_{B,10}=\{\infty ,-\infty ,6\}$.

Proposition 11.

Consider the PTAM as specified in Definition 3. Given m, the limiting distribution as $h_1\to 0$, $h_m\to \infty$, and $s\to -\infty$ is a gamma distribution with shape parameter $m-1$ and rate parameter λ. Denote this as sub-model ${\mathcal{M}}_{11}$ with corresponding boundary ${\mathit{r}}_{B,11}$. In this case, ${\mathcal{R}}_{11}=\{h_1,h_m,s\}$, ${\mathcal{P}}_{11}=\{\psi ,m\}$, and ${\mathit{r}}_{B,11}=\{0,\infty ,-\infty \}$.

Proposition 12.

Consider the PTAM as specified in Definition 3. Given m, the limiting distribution as $h_1\to 0$, $h_m\to \infty$, $s\to -\infty$, and $m\to 6$ is a gamma distribution with shape parameter 5 and rate parameter λ. Denote this as sub-model ${\mathcal{M}}_{12}$ with corresponding boundary ${\mathit{r}}_{B,12}$. In this case, ${\mathcal{R}}_{12}=\{h_1,h_m,s,m\}$, ${\mathcal{P}}_{12}=\left\{\psi \right\}$, and ${\mathit{r}}_{B,12}=\{0,\infty ,-\infty ,6\}$.

Proposition 13.

Consider the PTAM as specified in Definition 3, and $f(x;\mathit{\theta })$ is the associated p.d.f. According to Cheng et al. (2021), the limiting distribution as $m\to \infty$ is

$$\begin{array}{c}\hfill f(t;h_1,h_m,s,\psi )=e^{-{\int }_0^{t}h(u;h_1,h_m,s,\psi )\mathrm{d}u}h(t;h_1,h_m,s,\psi ),\end{array}$$

where

$$\begin{array}{c}\hfill h(t;h_1,h_m,s,\psi )=\left\{\begin{array}{cc}{\left((h_m^s-h_1^s)\frac{t}{\psi }+h_1^s\right)}^{\frac{1}{s}},\hfill & s\ne 0,\hfill \\ h_1^{1-\frac{t}{\psi }}{h}_m^{\frac{t}{\psi }},\hfill & s=0,\hfill \end{array}\right.\end{array}$$

(10)

is the hazard rate of the limiting distribution. Denote this as sub-model ${\mathcal{M}}_9$. In this case, ${\mathcal{R}}_{13}=\left\{m\right\}$, ${\mathcal{P}}_{13}=\{h_1,h_m,s,\psi \}$, and $r_{B,13}=r_{B,1}=\{\infty \}$.

Interestingly, the hazard rates that are given in (10) for $s\ne 0$ and $s=0$ correspond to the generalized Weibull distribution and the Gompertz law of mortality, respectively (Gompertz 1825; Pham and Lai 2007).

Proposition 14.

Consider the PTAM as specified in Definition 3. The limiting distribution under either one of the following cases is exponential with rate $h_m$:

(i) 

$h_1\to h_m$, $s\to \infty$, and $m\to \infty$, in which case ${\mathcal{R}}_{14}=\{h_1,s,\psi ,m\}$, ${\mathcal{P}}_{14}=\left\{h_m\right\}$, and ${\mathit{r}}_{B,14}=\{h_m,\infty ,\infty \}$.

(ii) 

$h_1\to h_m$, $s\to -\infty$, and $m\to \infty$, in which case ${\mathcal{R}}_{14}=\{h_1,s,\psi ,m\}$, ${\mathcal{P}}_{14}=\left\{h_m\right\}$, and ${\mathit{r}}_{B,14}=\{h_m,-\infty ,\infty \}$.

Denote this as sub-model ${\mathcal{M}}_{14}$ with corresponding boundary ${\mathit{r}}_{B,14}$.

Proposition 15.

Consider the PTAM as specified in Definition 3. The limiting distribution under either one of the following cases is exponential with rate $h_1$:

(i) 

$h_m\to h_1$, $s\to \infty$, and $m\to \infty$, in which case ${\mathcal{R}}_{15}=\{h_m,s,\psi ,m\}$, ${\mathcal{P}}_{15}=\left\{h_1\right\}$, and ${\mathit{r}}_{B,15}=\{h_1,\infty ,\infty \}$.

(ii) 

$h_m\to h_1$, $s\to -\infty$, and $m\to \infty$, in which case ${\mathcal{R}}_{15}=\{h_m,s,\psi ,m\}$, ${\mathcal{P}}_{15}=\left\{h_1\right\}$, and ${\mathit{r}}_{B,15}=\{h_1,-\infty ,\infty \}$.

Denote this as sub-model ${\mathcal{M}}_{15}$ with corresponding boundary ${\mathit{r}}_{B,15}$.

Since ${\mathcal{M}}_7$, ${\mathcal{M}}_8$, ${\mathcal{M}}_{14}$, and ${\mathcal{M}}_{15}$ are all exponential distributions, they can be grouped as one sub-model, which is denoted by $\mathcal{B}$. Combining all these sub-models, the sub-model family of the PTAM becomes

$$\begin{array}{c}\hfill \mathcal{F}:=\left(\bigcup _{k=1}^6{\mathcal{M}}_k\right)\bigcup \left(\bigcup _{k=9}^{13}{\mathcal{M}}_k\right)\bigcup \mathcal{B}.\end{array}$$

6.3. Simulation Study

In this section, the proposed definition of estimability is applied to assess the estimability of the PTAM via simulation studies.

We consider an underlying PTAM with $h_1=0.0007$, $h_m=0.1542$, $s=-2.2152,\lambda =0.3122$, and $m=10$. This underlying PTAM is utilized in Cheng (2021) and Cheng et al. (2021) to model the aging process of individuals residing in Channing House, a retirement community in Palo Alto, California (Hyde 1980). The Channing House data are chosen because all the residents in the community are subject to approximately the same circumstances, so, relatively speaking, the aging process is the most significant factor that contributes to the variability in their lifetimes, which is the process we intend to model using the PTAM. Furthermore, female data are chosen to preclude the effects of gender differences.

In

Table 3. c.d.f. sensitivity and experimental error for the PTAM—Experiments 1 to 6. The asterisk indicates that the PTAM is non-estimable.

Sample Size	Sub-Model	Sub-Model Boundary	100,000 IVs		200,000 IVs
			${\displaystyle \underset{\mathit{i}}{max}\left|\mathit{F}({\mathit{x}}_{\mathit{i}};\widehat{\mathit{\theta }})-{\mathit{G}}_{\mathit{k}}({\mathit{x}}_{\mathit{i}};\widehat{\mathit{p}})\right|}$	$\mathit{ϵ}$	${\displaystyle \underset{\mathit{i}}{max}\left|\mathit{F}({\mathit{x}}_{\mathit{i}};\widehat{\mathit{\theta }})-{\mathit{G}}_{\mathit{k}}({\mathit{x}}_{\mathit{i}};\widehat{\mathit{p}})\right|}$	$\mathit{ϵ}$
50	${\mathcal{M}}_1$	$r_{B,1}$	$0.23261$	$0.05807$	$0.22661$	$0.05310$
	${\mathcal{M}}_2$	${\mathit{r}}_{B,2}$	$0.22322$	$0.05807$	$0.21630$	$0.05310$
	${\mathcal{M}}_3$	$r_{B,3}$	0.02414 *	$0.05807$	${0.01698 }^{*}$	$0.05310$
	${\mathcal{M}}_4$	${\mathit{r}}_{B,4}$	$0.12693$	$0.05807$	$0.12693$	$0.05310$
	${\mathcal{M}}_5$	${\mathit{r}}_{B,5}$	0.05007 *	$0.05807$	0.05014 *	$0.05310$
	${\mathcal{M}}_6$	${\mathit{r}}_{B,6}$	$0.12054$	$0.05807$	$0.11965$	$0.05310$
	${\mathcal{M}}_9$	${\mathit{r}}_{B,9}$	0.04271 *	$0.05807$	0.01863 *	$0.05310$
	${\mathcal{M}}_{10}$	${\mathit{r}}_{B,10}$	$0.14590$	$0.05807$	$0.14300$	$0.05310$
	${\mathcal{M}}_{11}$	${\mathit{r}}_{B,11}$	0.04976 *	$0.05807$	0.04981 *	$0.05310$
	${\mathcal{M}}_{12}$	${\mathit{r}}_{B,12}$	$0.14360$	$0.05807$	$0.14360$	$0.05310$
	${\mathcal{M}}_{13}$	$r_{B,13}$	$0.16902$	$0.05807$	$0.12708$	$0.05310$
	$\mathcal{B}$	${\mathit{r}}_{B,7},{\mathit{r}}_{B,8},{\mathit{r}}_{B,14}$ or ${\mathit{r}}_{B,15}$	$0.36318$	$0.05807$	$0.35797$	$0.05310$
500	${\mathcal{M}}_1$	$r_{B,1}$	$0.32839$	$0.04051$	$0.31173$	$0.03845$
	${\mathcal{M}}_2$	${\mathit{r}}_{B,2}$	$0.21331$	$0.04051$	$0.21132$	$0.03845$
	${\mathcal{M}}_3$	$r_{B,3}$	0.02024 *	$0.04051$	0.01520 *	$0.03845$
	${\mathcal{M}}_4$	${\mathit{r}}_{B,4}$	$0.09340$	$0.04051$	$0.08392$	$0.03845$
	${\mathcal{M}}_5$	${\mathit{r}}_{B,5}$	0.02903 *	$0.04051$	0.03192 *	$0.03845$
	${\mathcal{M}}_6$	${\mathit{r}}_{B,6}$	$0.07648$	$0.04051$	$0.07648$	$0.03845$
	${\mathcal{M}}_9$	${\mathit{r}}_{B,9}$	0.02352 *	$0.04051$	0.02048 *	$0.03845$
	${\mathcal{M}}_{10}$	${\mathit{r}}_{B,10}$	$0.10873$	$0.04051$	$0.09800$	$0.03845$
	${\mathcal{M}}_{11}$	${\mathit{r}}_{B,11}$	0.02827 *	$0.04051$	0.02827 *	$0.03845$
	${\mathcal{M}}_{12}$	${\mathit{r}}_{B,12}$	$0.09986$	$0.04051$	$0.09987$	$0.03845$
	${\mathcal{M}}_{13}$	$r_{B,13}$	$0.30614$	$0.04051$	$0.15558$	$0.03845$
	$\mathcal{B}$	${\mathit{r}}_{B,7},{\mathit{r}}_{B,8},{\mathit{r}}_{B,14}$ or ${\mathit{r}}_{B,15}$	$0.33164$	$0.04051$	$0.34158$	$0.03845$
2000	${\mathcal{M}}_1$	$r_{B,1}$	$0.20239$	$0.02972$	$0.26365$	$0.01824$
	${\mathcal{M}}_2$	${\mathit{r}}_{B,2}$	$0.17499$	$0.02972$	$0.17335$	$0.01824$
	${\mathcal{M}}_3$	$r_{B,3}$	0.02612 *	$0.02972$	$0.02270$	$0.01824$
	${\mathcal{M}}_4$	${\mathit{r}}_{B,4}$	$0.05942$	$0.02972$	$0.05942$	$0.01824$
	${\mathcal{M}}_5$	${\mathit{r}}_{B,5}$	0.01999 *	$0.02972$	$0.02004$	$0.01824$
	${\mathcal{M}}_6$	${\mathit{r}}_{B,6}$	$0.04716$	$0.02972$	$0.04716$	$0.01824$
	${\mathcal{M}}_9$	${\mathit{r}}_{B,9}$	$0.03645$	$0.02972$	$0.02312$	$0.01824$
	${\mathcal{M}}_{10}$	${\mathit{r}}_{B,10}$	$0.07025$	$0.02972$	$0.07025$	$0.01804$
	${\mathcal{M}}_{11}$	${\mathit{r}}_{B,11}$	0.01936 *	$0.02972$	$0.01936$	$0.01824$
	${\mathcal{M}}_{12}$	${\mathit{r}}_{B,12}$	$0.06873$	$0.02972$	$0.06872$	$0.01824$
	${\mathcal{M}}_{13}$	$r_{B,13}$	$0.11637$	$0.02972$	$0.15976$	$0.01824$
	$\mathcal{B}$	${\mathit{r}}_{B,7},{\mathit{r}}_{B,8},{\mathit{r}}_{B,14}$ or ${\mathit{r}}_{B,15}$	$0.34331$	$0.02972$	$0.36079$	$0.01824$

, the proposed definition is applied to assess the model estimability with respect to six experiments:

• Experiment 1: A sample size of 50 with 100,000 initial values (IVs).
• Experiment 2: A sample size of 50 with 200,000 IVs.
• Experiment 3: A sample size of 500 with 100,000 IVs.
• Experiment 4: A sample size of 500 with 200,000 IVs.
• Experiment 5: A sample size of 2000 with 100,000 IVs.
• Experiment 6: A sample size of 2000 with 200,000 IVs.

The initial values (IVs) were randomly generated from the parameter space $\mathbf{\Theta }$. Those producing the smallest value of $g(·)$ in (11) were selected as the parameter estimates $\widehat{\mathit{\theta }}$. The following can be observed from Table 3:

(i)

According to the Theorems 1 and 2, the PTAM is non-estimable with respect to Experiments 1 to 5, whereas it is estimable with respect to Experiment 6.

(ii)

One potential barrier to the estimability of the PTAM occurs when the following two conditions are satisfied simultaneously: (i) $s\to -\infty$ and (ii) m is one of the model’s free parameters, which correspond to ${\mathcal{M}}_3$, ${\mathcal{M}}_5$, ${\mathcal{M}}_9$, and ${\mathcal{M}}_{11}$ in Figure 3.

(iii)

The algorithm noise decreases as the number of initial values increases. This again validates the proposed definition. However, this trend does not always hold unless the algorithm is selected, as is recommended in Section 7.1.

In conclusion, the non-estimability issue of the PTAM has been thoroughly investigated by applying the proposed definition. One may then arrive at an experimental design that makes the PTAM estimable, such as that utilized in Experiment 6.

7. Remarks and Conclusions

7.1. A Recommendation Regarding the Methodology

It is recommended that a parameter estimate ${\widehat{\mathit{\theta }}}^T$ be such that the following function g is minimized:

$$\begin{array}{c}\hfill {\displaystyle g\left(\mathit{\theta }\right)=\underset{x_i\in \mathit{x}}{max}\left|F_X(\mathit{\theta };x_i)-{\widehat{F}}_n\left(x_i\right)\right|,}\end{array}$$

(11)

where $F_X,{\widehat{F}}_n$, and $\mathit{x}$ are as specified in Definition 7. Otherwise, there will possibly exist implausible situations where the experimental error actually increases as the algorithm accuracy improves. To verify this, first notice that

$$\begin{array}{cc}\hfill ϵ& =g\left(\widehat{\mathit{\theta }}\right)\mathrm{and}{ϵ}_d=g\left({\widehat{\mathit{\theta }}}^T\right).\hfill \end{array}$$

Accordingly, if ${\widehat{\mathit{\theta }}}^T$ is based on optimization criteria other than (11), for instance, the maximization of the likelihood, then it is plausible that $ϵ=g\left(\widehat{\mathit{\theta }}\right)<g\left({\widehat{\mathit{\theta }}}^T\right)={ϵ}_d$ before improving the algorithm accuracy. This is due to the fact that ${ϵ}_d$ is not the minimum of g. Now, let ${ϵ}^{*}$ be the experimental error after the algorithm accuracy has improved. As explained in Section 5.3, increasing the algorithm accuracy does not change ${ϵ}_d$, but the experimental error converges to ${ϵ}_d$, making $ϵ<{ϵ}^{*}<{ϵ}_d$ possible, which is unreasonable. Such a situation can be avoided by determining ${\widehat{\mathit{\theta }}}^T$ via the minimization of g, in which case ${ϵ}_d$ is the minimum of g.

7.2. Future Research Avenues

The concept of estimability being introduced in this paper is entirely original. Two main avenues for future research are discussed in Section 7.2.1 and Section 7.2.2.

7.2.1. Refinements Regarding the c.d.f. Sensitivity

Interestingly, minimizing g is essentially equivalent to fitting a regression model to the ECDF ${\widehat{F}}_n$ while treating g as the loss function. This perspective then suggests that the concept of estimability can be redefined in terms of other loss functions. However, another valid g, if it exists, must also pass all the validation criteria presented in Section 5.1, Section 5.2 and Section 5.3. For example, the mean absolute deviation, that is,

$$\begin{array}{c}\hfill g\left(\mathit{\theta }\right)=\frac{1}{n}\sum _{i=1}^n\left|F_X(\mathit{\theta };x_i)-{\widehat{F}}_n\left(x_i\right)\right|,\end{array}$$

(12)

where $F_X,{\widehat{F}}_n$, and $\mathit{x}$ are as specified in Definition 7, is a valid option. Accordingly, the c.d.f. sensitivity in Definition 4 can then be modified as follows:

$$\begin{array}{c}\hfill \frac{1}{n}\sum _{i=1}^n|F(x_i;{\mathit{\theta }}_{\mathbf{1}})-F(x_i;{\mathit{\theta }}_{\mathbf{2}})|.\end{array}$$

Similarly, Definitions 5–9, 13 and Theorems 1 and 2 have to be modified to be consistent with (12) and then provide a valid alternative definition of estimability. Other definitions can also possibly exist for g. We are planning to investigate such research avenues. An alternative framework may yield a more precise definition of the concept of estimability, which can be readily applicable.

Such refinements of the definition of estimability can lead to improved strategies in situations where, for instance, the samples are small or one or more parameters lie on the boundary and certain linear/non-linear constraints are involved.

7.2.2. Refinements for Left-Truncated and Right-Censored Data

Another research avenue that we plan to pursue pertains to left-truncated and right-censored data. In such instances, the Kaplan–Meier survival curve replaces the ECDF. The experimental error specified in Definition 7 then becomes

$$\begin{array}{c}\hfill {\displaystyle ϵ:=\underset{x_i\in \mathit{x}}{max}\left|F(x_i;\widehat{\mathit{\theta }})-\left(1-{\widehat{S}}_n\left(x_i\right)\right)\right|,}\end{array}$$

(13)

where ${\widehat{S}}_n\left(x_i\right)$ denotes the Kaplan–Meier estimates. The remaining definitions and Theorems 1 and 2 should then be modified accordingly.

This generalizes the proposed methodology since $1-{\widehat{S}}_n\left(x_i\right)$ then boils down to the ECDF, ${\widehat{F}}_n\left(x_i\right)$, without left truncation and right censoring. Furthermore, ${\widehat{S}}_n\left(x_i\right)$ converges in distribution to $S\left(x_i\right)$ as the consistency property still holds for the Kaplan–Meier estimator, and thus, the validation process described in Section 5.1 also applies. Since real data often involve left truncation and right censoring, once related refinements are incorporated in the definition of estimability, we shall apply it to carry out further studies involving actual data sets, such as the Channing House data set, which was previously analyzed by Hyde (1980).

7.3. Concluding Remarks

A novel definition of estimability is proposed to objectively quantify this concept in the context of statistical models. The proposed definition is akin to an existing method in that the non-estimability is defined as occurring when the confidence region is infinite. However, unlike this method, which utilizes a likelihood-based confidence region, a c.d.f. sensitivity-based confidence region is proposed with a view of linking the confidence region to the experimental protocol. Under this setting, the threshold becomes objective as the experimental error becomes an experiment-based quantity under the proposed methodology. The validated definition is then applied to assess the estimability of two statistical models and the PTAM. In the latter case, this solves a potential non-estimability issue. The calculations were carried out with the R software package, the code being available upon request from the first author.