Probabilistic sensitivity analysis on Markov models with uncertain transition probabilities: an application in evaluating treatment decisions for type 2 diabetes

Abstract

Markov models are commonly used for decision-making studies in many application domains; however, there are no widely adopted methods for performing sensitivity analysis on such models with uncertain transition probability matrices (TPMs). This article describes two simulation-based approaches for conducting probabilistic sensitivity analysis on a given discrete-time, finite-horizon, finite-state Markov model using TPMs that are sampled over a specified uncertainty set according to a relevant probability distribution. The first approach assumes no prior knowledge of the probability distribution, and each row of a TPM is independently sampled from the uniform distribution on the row’s uncertainty set. The second approach involves random sampling from the (truncated) multivariate normal distribution of the TPM’s maximum likelihood estimators for its rows subject to the condition that each row has nonnegative elements and sums to one. The two sampling methods are easily implemented and have reasonable computation times. A case study illustrates the application of these methods to a medical decision-making problem involving the evaluation of treatment guidelines for glycemic control of patients with type 2 diabetes, where natural variation in a patient’s glycated hemoglobin (HbA1c) is modeled as a Markov chain, and the associated TPMs are subject to uncertainty.

Appendices

Appendix 1: Convergence properties of algorithm 1 when using algorithm 2 or 3

To describe the overall framework for the following analysis, we build on the setup for Section 2.3 in which (a) q_i has the uncertainty set $\mathcal { Q}_{i}$ contained in the hyperplane $\mathcal { H}\subset \mathbb {R}^{n}$, and $\mathcal { H}$ has the relative topology inherited from $\mathbb {R}^{n}$; and (b) $\mathcal { Q}_{i}$ is an open, bounded, convex subset of $\mathcal { H}$. Corresponding to the true TPMs in $\mathcal { P}^{\,{\boldsymbol {\pi }}}$, we let $\mathcal { U}^{{\boldsymbol {\pi }}}$ denote the uncertainty set in $\mathbb {R}^{\tau }$ (where Eq. 4 implies that $\tau \leq n^{2}(T-1)|\mathcal { A}|$); and we let $\mathcal { H}^{\boldsymbol {\pi }} \subset \mathbb {R}^{\tau }$ denote the subspace in which Algorithms 2 and 3 operate, where $\mathcal { U}^{{\boldsymbol {\pi }}} \subset \mathcal { H}^{\boldsymbol {\pi }}$ and $\mathcal { H}^{\boldsymbol {\pi }}$ has the relative topology inherited from $\mathbb {R}^{\tau }$. Hence $\mathcal { U}^{{\boldsymbol {\pi }}}$ is an open, bounded, convex subset of $\mathcal { H}^{\boldsymbol {\pi }}$; and in $\mathbb {R}^{\tau }$, we see that $\mathcal { H}^{\boldsymbol {\pi }}$ is the Cartesian product of at most $n(T-1)|\mathcal { A}|$ copies of $\mathcal { H}$ [15, pp. 98–100].

Each point $\boldsymbol {u}\in \mathcal { U}^{{\boldsymbol {\pi }}}$ represents an assignment of values to the elements of the TPMs used by the MDP (1) such that Eq. 5 is satisfied for each TPM; and we let $\mathbb {V}_{\boldsymbol {\pi }}(\boldsymbol {u})$ denote the resulting expected total reward (3) based on policy π. Let $\bar {\mathcal { U}}^{{\boldsymbol {\pi }}}$ denote the closure of $\mathcal { U}^{{\boldsymbol {\pi }}}$ in $\mathcal { H}^{\boldsymbol {\pi }}$ so that $\bar {\mathcal {U}}^{{\boldsymbol {\pi }}}$ includes the boundary points of $\mathcal { U}^{{\boldsymbol {\pi }}}$ [15, pp. 69–72]. Thus $\bar {\mathcal { U}}^{{\boldsymbol {\pi }}}$ is a closed, bounded, convex set in $\mathcal { H}^{\boldsymbol {\pi }}$ so that $\bar {\mathcal { U}}{~}^{{\boldsymbol {\pi }}}$ is also compact in $\mathcal { H}^{\boldsymbol {\pi }}$ [15, p. 233]. Equations 2 and 3 imply that $\mathbb {V}_{\boldsymbol {\pi }}(\cdot )$ is a polynomial function of its arguments, so it is continuous on $\mathbb {R}^{\tau }$; thus $\mathbb {V}_{\boldsymbol {\pi }}(\cdot )$ is also continuous when its domain is restricted to $\bar {\mathcal { U}}{~}^{{\boldsymbol {\pi }}}$ [15, p. 79]. Let $\underline {v}^{{\boldsymbol {\pi }}}$ ≡$\min \!\left \{\, \mathbb {V}_{\boldsymbol {\pi }}(\boldsymbol {u}) \, |\, \boldsymbol {u} \in \bar {\mathcal { U}}^{{\boldsymbol {\pi }}}\, \right \},$ the minimum value of Eq. 3 over $\bar {\mathcal { U}}^{{\boldsymbol {\pi }}}$ [15, p. 227]. Similarly let $\overset {\sim }{\smash {{\underline {v}}}}_{M}$ ≡$\min \! \left \{\, \mathbb {V}_{\boldsymbol {\pi }}\left (\boldsymbol {U}_{j} \right )\, |\, j = 1,\ldots , M \,\right \}$, the minimum observed value of Eq. 3 over the random vectors $\{ \boldsymbol {U}_{j} \in \mathcal { U}^{{\boldsymbol {\pi }}} \,|\, j = 1, \ldots , M\}$ corresponding to the sample $\left \{\, \overset {\sim }{\smash {\mathcal { P}}}_{j}^{{\boldsymbol {\pi }}}\, |\, j = 1,\ldots , M \,\right \}$.

Proof of Property A ₁

When Algorithm 1 uses Algorithm 2 to generate {U_j | j = 1,…,M},we show that with probability one $\,\!\overset {\sim }{\smash {{\underline {v}}}}_{M} \rightarrow \underline {v}^{{\boldsymbol {\pi }}}$as M →∞.Given any ε > 0,let $g_{M} \equiv \Pr \left \{\,\!\overset {\sim }{\smash {{\underline {v}}}}_{M} > \underline {v}^{{\boldsymbol {\pi }}} + \varepsilon \,\right \}$for M ≥ 1.Because $\mathbb {K} \equiv \mathbb {V}_{\boldsymbol {\pi }}^{-1}\left \{\,\left (\,\underline {v}^{{\boldsymbol {\pi }}},\, \underline {v}^{{\boldsymbol {\pi }}}+\varepsilon \,\right )\,\right \}$is an open subset of $\bar {\mathcal { U}}{~}^{{\boldsymbol {\pi }}}$in $\mathcal { H}^{\boldsymbol {\pi }}$and because $\bar {\mathcal { U}}{~}^{{\boldsymbol {\pi }}}$is compact in $\mathcal { H}^{\boldsymbol {\pi }}$,the sampling scheme of Algorithm 2 ensures that for each ℓ ∈{1,…,M − 1}and $\boldsymbol {u} \in \bar {\mathcal { U}}{~}^{{\boldsymbol {\pi }}}$,the function $h(\boldsymbol {u}) \equiv \Pr \left \{ \boldsymbol {U}_{\ell + 1}\;\, \not \in \mathbb {K} \, | \, \boldsymbol {U}_{\ell } =\boldsymbol {u} \, \right \}$has the following properties: (a) it depends on u but not on ℓ;(b) it is continuous at u;and (c) it attains a lower limit η > 0and an upper limit δ < 1on $\bar {\mathcal { U}}{~}^{{\boldsymbol {\pi }}}$[15, p. 227]. For $\boldsymbol {u} \in \mathbb {R}^{\tau }$and ℓ = 0,…,M,we let F_(ℓ)(u)denote thec.d.f.of U_ℓ,where $\boldsymbol {U}_{\,0} \in \mathcal { U}^{{\boldsymbol {\pi }}}$is the user-specified (fixed) initial point of Algorithm 2 so that thec.d.f. F₍₀₎(u)is degenerate [7, p. 193]. Note that $\bar {\mathcal { U}}{~}^{{\boldsymbol {\pi }}}$is the support of U_ℓfor ℓ ≥ 0[7, p. 161]. Property (c) of h(⋅)and the law of total probability [7, Eq. (33.8)] imply that

$$ \Pr\!\left\{\, \boldsymbol{U}_{\ell} \not \in \mathbb{K}\, \right\} = {\int}_{\bar{\mathcal{ U}}{~}^{{\boldsymbol{\pi}}}} h(\boldsymbol{u}) \, \text{d} F_{(\ell-1)}(\boldsymbol{u}) \geq {\eta} > 0 \text{\ \:for\ \:} \ell = 1, \ldots, M\, . $$

(A.1)

Since $g_{j} = \Pr \{\boldsymbol {U}_{\ell } \not \in \mathbb {K} \text {\ for\ } \ell = 1, \ldots , j \}$, we proceed byshowing the intermediate result g_j ≤ δ^j− 1for j = 1,…,M using inductionon j. Clearly g₁ ≤ 1.Then for M ≥ 2and ℓ = 1,…,M − 1, thelaw of total probability implies that

$$ \Pr\!\left\{ \boldsymbol{U}_{\ell+ 1}\;\, \not\in \mathbb{K}\, \text{\ and\ } \boldsymbol{U}_{\ell} \not \in \mathbb{K} \right\} = {\int}_{\bar{\mathcal{ U}}{~}^{{\boldsymbol{\pi}}}\setminus\mathbb{K}} h(\boldsymbol{u})\,\text{d} F_{(\ell)}(\boldsymbol{u}) \leq {\delta} \Pr\!\left\{ \boldsymbol{U}_{\ell} \not \in \mathbb{K} \right\} \,; $$

and by Eq. A.1, we have

$$ \Pr\!\left\{\boldsymbol{U}_{\ell+ 1} \,\,\not \in \mathbb{K} \, |\, \boldsymbol{U}_{\ell} \not \in \mathbb{K} \right\} \leq {\delta} < 1 \text{\ for\ } \ell = 1, \ldots, M-1 \, . $$

(A.2)

From Eq. A.2 for ℓ = 1 and M = 2, we see that $g_{2} = g_{1}\Pr \!\left \{ \,\boldsymbol {U}_{\,2}\not \in \mathbb {K}\, | \, \boldsymbol {U}_{1}\not \in \mathbb {K}\, \right \} \leq {\delta }$.If M ≥  3 and j =  2,…, M − 1, we assume byinduction that g_j ≤ δ^j− 1.To continue, we observe that the sampling scheme of Algorithm 2 ensures the process {U_j | j =  1,…, M} has theMarkov property. We have

$$\begin{array}{@{}rcl@{}} g_{j + 1} &=& {\int}_{\bar{\mathcal{ U}}{~}^{{\boldsymbol{\pi}}}\setminus\mathbb{K}} \Pr\left\{\, \boldsymbol{U}_{\ell} \not\in\mathbb{K} \text{\ for\ } \ell = 1, \ldots, j-1 \text{\ and\ } \boldsymbol{U}_{j + 1}\not\in\mathbb{K}\,|\, \boldsymbol{U}_{j} = \boldsymbol{u} \,\right\} \, \text{d} F_{(j)}(\boldsymbol{u}) \end{array} $$

(A.3)

$$\begin{array}{@{}rcl@{}} &=&{\int}_{\bar{\mathcal{ U}}{~}^{{\boldsymbol{\pi}}}\setminus\mathbb{K}} \Pr\left\{\, \boldsymbol{U}_{\ell} \not\in\mathbb{K} \text{\ for\ } \ell = 1, \ldots, j-1\,|\, \boldsymbol{U}_{j} = u \,\right\} h(\boldsymbol{u}) \, \text{d} F_{(j)}(\boldsymbol{u}) \end{array} $$

(A.4)

$$\begin{array}{@{}rcl@{}} &\leq& {\delta} {\int}_{\bar{\mathcal{ U}}{~}^{{\boldsymbol{\pi}}}\setminus\mathbb{K}} \Pr\left\{\, \boldsymbol{U}_{\ell} \not\in\mathbb{K} \text{\ for\ } \ell = 1, \ldots, j-1\,|\, \boldsymbol{U}_{j} = \boldsymbol{u} \,\right\}\,\text{d} F_{(j)}(\boldsymbol{u}) \\ &=&{\delta} g_{j} \leq {\delta}^{j}\,, \end{array} $$

(A.5)

where (A.3) follows by the law of total probability, Eq. A.4 follows from the Markov property[38, p. 112], and Eq. A.5 holds by the induction hypothesis. Therefore we have g_M ≤ δ^M− 1 → 0for M →∞ so we can apply the same line of reasoning as in Remark 2. Hence with probability one there exists M^∗such that$\overset {\sim }{\smash {{\underline {v}}}}_{M} \leq \underline {v}^{\boldsymbol {\pi }} + \varepsilon$for all M ≥ M^∗. Since ε is arbitrary, it follows thatwith probability one $\,\!\overset {\sim }{\smash {{\underline {v}}}}_{M} \rightarrow \underline {v}^{{\boldsymbol {\pi }}}$as M →∞. Similarly we can prove convergencewith probability one to ${\bar {v}}^{\,{\boldsymbol {\pi }}} = \max \!\left \{\, \mathbb {V}_{\boldsymbol {\pi }}(\boldsymbol {u}) \, |\, \boldsymbol {u} \in \bar {\mathcal { U}}{~}^{{\boldsymbol {\pi }}}\, \right \}$ for thecorresponding sample statistic based on {U_j | j = 1,…,M}.

Next we prove the comparable convergence result for Algorithm 1 when it uses Algorithm 3 for sampling. The notation of theprevious paragraphs is reused. In this case the sampling scheme of Algorithm 3 ensures that the random vectors {U_j | j = 1,…,M} are i.i.d. Based on Eqs. 14and 15, first we show that the {U_j | j = 1,…,M} can be expressed in terms of i.i.d. samples from a truncated u-variate normal distribution for a suitable dimension$u \leq (n-1)n(T-1)|\mathcal { A}|$. Let Z = [Z₁,…,Z_u] denote a standard normal random vector with probability density function (p.d.f.)$\phi _{u}(\boldsymbol {z}) = (2\pi )^{-u/2}\exp \left (-\frac {1}{2}{\sum }_{i = 1}^{u} {z_{i}^{2}}\right )$for$\boldsymbol {z} \in \mathbb {R}^{u}$. Let$\boldsymbol {\zeta }: \mathbb {R}^{u} \mapsto \mathcal { H}^{\boldsymbol {\pi }}$ denote the continuous function representing the u-dimensional generalization of Steps 5 and 6 in Algorithm 3. Since$\mathcal { U}^{{\boldsymbol {\pi }}}$is open in$\mathcal { H}^{\boldsymbol {\pi }}$, we seethat $\mathcal { F} \equiv \boldsymbol {\zeta }^{-1}(\mathcal { U}^{{\boldsymbol {\pi }}})$ is open in $\mathbb {R}^{u}$.

By randomly sampling the { A_j | j = 1,…,M }from the u-variate truncated normal p.d.f.

$$ f_{{\boldsymbol{A}}}(\boldsymbol{a}) \equiv \frac{\phi_{u}(\boldsymbol{a})}{{\int}_{\mathcal{ F}} \phi_{u}(\boldsymbol{z}) \, \text{d} \boldsymbol{z} } \text{\ \:for\:\ } \boldsymbol{a} \in \mathcal{ F}\, , $$

(A.6)

we obtain the i.i.d.random vectors { ζ(A_j) | j = 1,…,M }on $\mathcal { U}^{{\boldsymbol {\pi }}}$. By using an argument similar to Remark 2, we see that Steps 4–10 of Algorithm 3 generatei.i.d.random vectors {U_j | j = 1,…, M} having the samedistribution as the { ζ(A_j) | j = 1,…,M }.Since $\mathbb {K}$ is openin $\mathcal { H}^{\boldsymbol {\pi }}$, we seethat $\mathbb {L} \equiv \boldsymbol {\zeta }^{-1}(\mathbb {K})$ is open in $\mathbb {R}^{u}$. Because f_A(a) >  0 for all $\boldsymbol {a} \in \mathcal { F}$, we have

$$\begin{array}{@{}rcl@{}} \Pr\left\{\, \boldsymbol{U}_{j} \in \mathbb{K}\, \right\} &=& \Pr\left\{\, \boldsymbol{\zeta}({\boldsymbol{A}}_{j}) \in \mathbb{K}\, \right\} = \Pr\left\{\, {\boldsymbol{A}}_{j} \in \mathbb{L} \, \right\}\\ &=& {\int}_{\mathbb{L}} f_{{\boldsymbol{A}}}(\boldsymbol{a})\, \text{d} \boldsymbol{a} > 0 \end{array} $$

[37, p. 80, Eq. (10)]. Hence $\delta \equiv \Pr \left \{\, \boldsymbol {U}_{j} \not \in \mathbb {K}\, \right \} < 1$for j = 1,…,M so that g_M = δ^M → 0as M →∞. As in the previous paragraph, it follows that with probability one $\,\!\overset {\sim }{\smash {{\underline {v}}}}_{M} \rightarrow \underline {v}^{{\boldsymbol {\pi }}}$as M →∞. Similarly we can prove convergence with probability one to${\bar {v}}^{\,{\boldsymbol {\pi }}} = \max \left \{\, \mathbb {V}_{\boldsymbol {\pi }}(\boldsymbol {u}) \, |\, \boldsymbol {u} \in \bar {\mathcal { U}}{~}^{{\boldsymbol {\pi }}}\, \right \}$ for the corresponding sample statistic based on {U_j | j = 1,…,M}. □

Proof of Property A ₂

Finally we consider the use of Algorithm 1 to estimate the distribution of theexpected total reward (3) based on a given distribution on the uncertainty set$\mathcal { U}^{{\boldsymbol {\pi }}}$.When Algorithm 1 uses Algorithm 3 for sampling, the random variables $\left \{\, Y_{j} = \mathbb {V}_{\boldsymbol {\pi }}\left (\boldsymbol {U}_{j}\right ) \, |\, j = 1, \ldots , M\right \}$are independent with a common c.d.f. F(y) for $y \in \mathbb {R}$; and letting $\widehat {F}_{\!\!M}(y)$(for $y \in \mathbb {R}$) denote the empirical c.d.f. for this random sample of size M, we see that with probability one$\widehat {F}_{\!\!M}(y) \rightarrow F(y)$ uniformly in y as M →∞[7, p. 269]. Moreover the estimator $\widehat {F}_{\!\!M}(y)$ has mean F(y) and variance F(y)[1 − F(y)]/M; and for M sufficiently large, $\widehat {F}_{\!\!M}(y)$is approximately normal [7, p. 358, Ex. 27.1]. □

Proof of Property A ₃

When Algorithm 1 uses Algorithm 2 to sample the { U_j | j = 1,…, M} from the uniform distribution on $\mathcal { U}^{{\boldsymbol {\pi }}}$, the random variables $\left \{\, Y_{j} = \mathbb {V}_{\boldsymbol {\pi }}\left (\, \boldsymbol {U}_{j}\,\right )\,| \, j = 1, \ldots , M\, \right \}$are generally neither independent or identically distributed. However, Remark 2 ensures that U_jconverges in distribution to a random vector U^∗ that is uniformly distributed on $\mathcal { U}^{{\boldsymbol {\pi }}}$ as j → ∞; and because $\mathbb {V}_{\boldsymbol {\pi }}(\cdot )$is a continuous function on $\mathcal { U}^{{\boldsymbol {\pi }}}$,we see that Y_jconverges in distribution to the random variable $Y^{*} \equiv \mathbb {V}_{\boldsymbol {\pi }}(\boldsymbol {U}^{*})$having thec.d.f. F^∗(y)(for $y \in \mathbb {R}$)as j →∞(see [7, pp. 327–329] and [6, Theorem B.7.2 (c)]). Let$\widehat {F}_{\!\!M}^{*}(y)$ denote the empirical c.d.f. based on the sample { Y_j | j = 1,…,M }. Hence beyond a sufficiently long warm-up period for Algorithm 2, we may regard the process { Y_j | j = 1,…, M} as approximately covariance stationary [24, p. 445] with marginal c.d.f. F^∗(⋅). In this situation, under certain widely applicable conditions (see [9, Eq. (5.1.3)] and [24, p. 476]), for each $y\in \mathbb {R}$ separately we see that $\widehat {F}_{\!\!M}^{*}(y)$ converges in mean square and hence in probability to F^∗(y)as M →∞[33, pp. 258–260]. □

Appendix 2: MLEs and CI estimators of HbA1c transition probabilities

Table 4 The MLE, sample sizes, and lower/upper bounds of transition probabilities’ 99% overall confidence intervals for male patients

Table 5 The MLE, sample sizes, and lower/upper bounds of transition probabilities’ 99% overall confidence intervals for female patients