# Robust explicit estimation of the log-logistic distribution with
applications

Zhuanzhuan Maa and Min Wangb,111Corresponding author. Email:
<EMAIL_ADDRESS>and Chanseok Parkc
a School of Mathematical and Statistical Sciences, The University of Texas Rio
Grande Valley, Brownsville, TX
bDepartment of Management Science and Statistics, The University of Texas at
San Antonio, TX
cDepartment of Industrial Engineering, Pusan National University, Busan, Korea

###### Abstract

The parameters of the log-logistic distribution are generally estimated based
on classical methods such as maximum likelihood estimation, whereas these
methods usually result in severe biased estimates when the data contain
outliers. In this paper, we consider several alternative estimators, which not
only have closed-form expressions, but also are quite robust to a certain
level of data contamination. We investigate the robustness property of each
estimator in terms of the breakdown point. The finite sample performance and
effectiveness of these estimators are evaluated through Monte Carlo
simulations and a real-data application. Numerical results demonstrate that
the proposed estimators perform favorably in a manner that they are comparable
with the maximum likelihood estimator for the data without contamination and
that they provide superior performance in the presence of data contamination.

Key words: Breakdown point, robustness, outlier, percentiles, repeated median
estimator, Hodges-Lehmann estimator

2000 MSC: 62F10, 62F35

##  1 Introduction

The log-logistic (LL) distribution, also called the Fisk distribution in
economics, is related to the logistic distribution in an identical fashion to
how the log-normal and normal distributions are related to each other. The LL
distribution can be generated by taking an exponential transformation on the
logistic distribution. A random variable $T$ is said to follow the LL
distribution with the scale parameter $\alpha>0$ and the shape parameter
$\beta>0$, denoted by $T\sim\mathrm{LL}(\alpha,~{}\beta)$, if its cumulative
distribution function (cdf) can be written as

$F(t\mid\alpha,\beta)=\frac{t^{\beta}}{\alpha^{\beta}+t^{\beta}},$ (1)

where $t>0$, and its probability density function (pdf) takes the form

$f(t\mid\alpha,\beta)=\frac{\beta\alpha^{\beta}t^{\beta-1}}{(t^{\beta}+\alpha^{\beta})^{2}}.$

It is worth noting that the LL distribution can be viewed as one of the
commonly used parametric distributions for survival analysis. Unlike the well-
known Weibull distribution, whose hazard rate function is $h(x\mid
k,\lambda)=k\lambda^{-k}x^{k-1}$ with scale parameter $\lambda>0$ and shape
parameter $k>0$, the LL distribution has a non-monotonic hazard rate function
given by

$h(t\mid\alpha,\beta)=\frac{(\beta/\alpha)(t/\alpha)^{\beta-1}}{1+(t/\alpha)^{\beta}},$

which can increase initially and decrease later and at times can be hump-
shaped. In addition, by allowing $\alpha$ to differ among groups, it can be
adopted as the basis of an accelerated failure time (AFT) model. These
statistical properties make the LL distribution be a favorable distribution in
a variety of fields, such as economics (Fisk 1961), hydrology (Shoukri et al.
1988), medical and biological sciences (Geskus 2001), and engineering (Ashkar
and Mahdi 2003).

The parameters of the LL distribution are generally estimated based on
classical methods such as maximum likelihood and least squares estimation
methods. For instance, Balakrishnan and Malik (1987) studied the method of
moments for the parameters of the truncated LL distribution. Kantam and
Srinivasa (2002) studied the modified maximum likelihood estimation of
$\alpha$ while assuming a known value of $\beta$. Reath et al. (2018) proposed
unbiased or nearly unbiased estimators for the parameters of the LL
distribution. He et al. (2020a) reconsidered the maximum likelihood estimators
based on simple random sampling and ranked set sampling techniques. More
recently, He et al. (2020b) studied modified best linear unbiased estimator
for estimating $\beta$. These existing estimators perform well for a complete
data without any contamination, whereas they may become unreliable and can
result in severely distorted estimates when the data contain outliers. A small
amount of data contamination even a single outlying observation could induce a
large impact on these estimators and even make them breakdown.

However, in many practical applications, there is no guarantee that the
collected data are free of any contamination due to various volatile operating
conditions when collecting data from an operating system. These observations
motivate us to consider several outlier-resistant estimators for the
parameters of the LL distribution, which include the percentile estimators,
the repeated median estimators, the sample median and median absolute
deviation estimators, the Hodges-Lehmann and Shamos estimators. It is worth
noting that these estimators not only have simple closed-form expressions, but
also are quite robust to a certain level of data contamination. Furthermore,
we investigate the breakdown point (Hampel et al. 1986) of each estimator to
measure its robustness in the presence of data contamination. Here, the
breakdown point is defined as the proportion of incorrect observations (i.e.
either arbitrarily large or small observations), the estimator of a parameter
can deal with before yielding estimated values arbitrarily close to zero
(implosion) or infinity (explosion). For example, the mean has a breakdown
point of $0\%$ and the median has a breakdown point of $50\%$.

The remainder of this paper is organized as follows. In Section 2, we review
maximum likelihood estimation (MLE) for the parameters of the LL distribution
and provide a list of alternative outlier-resistant estimators. In Section 3,
we conduct extensive Monte Carlo simulations to investigate the finite sample
performance of each estimator under consideration with and without data
contamination. In Section 4, we apply the proposed estimators to a real data
example for illustrative purposes. Finally, some concluding remarks are
provided in Section 5 with a proof deferred to the appendix.

##  2 Parameter estimation methods

Let $t_{1},t_{2},\cdots,t_{n}$ be $n$ observations from the LL distribution
with the parameters $\alpha$ and $\beta$ in (1). The log-likelihood function
can be written as

$\log
L=n\log\left(\beta\right)-n\log\left(\alpha\right)+\left(\beta-1\right)\sum_{i=1}^{n}\log
t_{i}-2\sum_{i=1}^{n}\log\left[1+\left(\frac{t_{i}}{\alpha}\right)^{\beta}\right].$

By taking the first derivative of the above log-likelihood function with
respect to $\alpha$ and $\beta$, respectively, and then equating them to zero,
the MLEs of $\alpha$ and $\beta$, denoted by $\hat{\alpha}$ and $\hat{\beta}$,
can be obtained by finding the solutions of the following two equations

$\displaystyle
2\sum_{i=1}^{n}\frac{\left(t_{i}/\alpha\right)^{\beta}}{1+\left(t_{i}/\alpha\right)^{\beta}}-n=0,$
$\displaystyle
2\beta\sum_{i=1}^{n}\frac{\left(t_{i}/\alpha\right)^{\beta}\log\left(t_{i}/\alpha\right)}{1+\left(t_{i}/\alpha\right)^{\beta}}-\beta\sum_{i=1}^{n}\log\left(t_{i}/\alpha\right)-n=0.$

Given that no explicit solutions are available to the two equations above, we
can numerically obtain the MLEs $\hat{\alpha}$ and $\hat{\beta}$ by using the
mledist($\cdot$) function from the R fitdistrplus package (Delignette-Muller
and Dutang 2015).

As mentioned above, the MLEs work well for estimating the parameters of the LL
distribution when there is no data contamination, whereas they may result in
unreliable and severely distorted parameter estimates in the presence of data
contamination (even a single outlying observation); see, for eample, Park and
Cho (2003); Park et al. (2020); Huber and Ronchetti (2009), among others. This
observation motivates us to develop several alternative outlier-resistant
estimators, including percentile estimators, repeated median estimators,
sample median and median absolute deviation estimators, Hodges-Lehmann and
Shamos estimators for the parameters of the LL distribution.

###  2.1 Percentile estimators

By taking the inversion of LL cdf in (1), we obtain

$t=\alpha\left[\frac{F(t\mid\alpha,\beta)}{1-F(t\mid\alpha,\beta)}\right]^{1/\beta}.$
(2)

To derive the percentile estimators, we may use the first ($t_{0.25}$) and
third ($t_{0.75}$) quartiles of $n$ observations to solve for $\alpha$ and
$\beta$ in (2) and obtain the following two equations

$t_{0.25}=\alpha\left(\dfrac{0.25}{1-0.25}\right)^{1/\beta}\ \mathrm{and}\
t_{0.75}=\alpha\left(\dfrac{0.75}{1-0.75}\right)^{1/\beta},$

which provide the percentile estimators of $\alpha$ and $\beta$ given by

$\hat{\alpha}_{p}=t_{0.25}\left(\dfrac{0.25}{1-0.25}\right)^{-1/\hat{\beta}_{p}}\
\mathrm{and}\
\hat{\beta}_{p}=\dfrac{2\left[\log\left(0.75\right)-\log\left(0.25\right)\right]}{\log\left(t_{0.75}/t_{0.25}\right)},$
(3)

respectively. For any given pair of $n$ observations: low ($100l^{th}$) and
high ($100h^{th}$), we can obtain the percentile estimators of $\alpha$ and
$\beta$ given by

$\displaystyle\hat{\alpha}_{p}$
$\displaystyle=t_{l}\left(\dfrac{l}{1-l}\right)^{-{1}/{\hat{\beta}_{p}}}\
\mathrm{and}\ \hat{\beta}_{p}$
$\displaystyle=\dfrac{2\left[\log\left(h\right)-\log\left(l\right)\right]}{\log\left(t_{h}/t_{l}\right)},$
(4)

respectively. To study the robustness property of the percentile estimators in
(4) based on contaminated data, we assume that all observations from the LL
distribution are distinct, which however happens with probability one, since
the LL distribution is continuous. We obtain the breakdown points of the
estimators in (4) summarized in the following proposition with its proof
provided in the appendix.

###### Proposition 1

The asymptotic breakdown point of the symmetrical percentile estimators
$\hat{\alpha}_{p}$ of $\alpha$ is given by

$\kappa=\min(1-h,2h-1).$

The highest breakdown point of $\hat{\alpha}_{p}$ is $1/3$ when $l=1/3$ and
$h=2/3$, and the corresponding highest breakdown point of the percentile
estimator $\hat{\beta}_{p}$ of $\beta$ with $\hat{\alpha}_{p}$ is $2/3$.

###  2.2 Repeated median estimators

It is worth noting that we rewrite the cdf of the LL distribution in (1) as a
linear regression form

$y_{i}=\beta_{0}+\beta_{1}x_{i}\quad\mathrm{for}~{}i=1,\cdots,n,$ (5)

where
$y_{i}=\log\left[\left(1-F(t_{(i)}\mid\alpha,\beta)\right)^{-1}-1\right]$,
$x_{i}=\log\left(t_{(i)}\right)$, $\beta_{0}=-\beta\log\alpha$,
$\beta_{1}=\beta$, and $t_{(1)},\cdots,t_{(n)}$ are the order statistics (from
the smallest to the largest) of $n$ observations from the LL distribution. To
approximate $F(t_{(i)}\mid\alpha,\beta)$ in $y_{i}$, we adopt the technique of
Ross (2000) with $i/(n+1)$, which can be easily implemented by using the
ppoints($\cdot$) function in R language (R Core Team 2018). The least squares
estimates of $\beta_{1}$ and $\beta_{0}$ in (5) are given by

$\hat{\beta}_{1}=\dfrac{\sum_{i=1}^{n}(x_{i}-\bar{x_{i}})(y_{i}-\bar{y_{i}})}{\sum_{i=1}^{n}(x_{i}-\bar{x_{i}})^{2}}\
\mathrm{and}\ \hat{\beta}_{0}=\bar{y}-\hat{\beta}_{1}\bar{x},$

respectively. Analogous to the MLEs, they are not robust to data contamination
with the breakdown points of $0\%$. For this reason, we may consider the RM
estimators (Siegel 1982) for estimating $\beta_{0}$ and $\beta_{1}$ in (5),
which are given by

$\hat{\beta}_{1}=\mathop{\mathrm{med}}_{1\leq i\leq
n}\mathop{\mathrm{med}}_{j\neq
i}\frac{y_{i}-y_{j}}{x_{i}-x_{j}}\quad\mathrm{and}\quad\hat{\beta}_{0}=\mathop{\mathrm{med}}_{1\leq
i\leq n}\big{(}y_{i}-\hat{\beta}_{1}x_{i}\big{)},$ (6)

where med${}_{i}(x_{i})$ is the sample median for the $n$ observations
$x_{1},\cdots,x_{n}$. Then we obtain the RM estimators of the parameters of
the LL distribution given by

$\hat{\alpha}_{\mathrm{RM}}=\exp\left(-\dfrac{\hat{\beta}_{0}}{\hat{\beta}_{1}}\right)\quad\mathrm{and}\quad\hat{\beta}_{\mathrm{RM}}=\hat{\beta}_{1},$
(7)

respectively. It deserves mentioning that the breakdown points of the RM
estimator in (7) are equal to $50\%$, since the breakdown points of
$\hat{\beta}_{1}$ and $\hat{\beta}_{0}$ in (6) are $50\%$.

###  2.3 Sample median and median absolute deviation estimators

As mentioned in the introduction, if $Z=\log T$, then $Z$ is said to follow a
logistic regression distribution with the location parameter $\mu=\log\alpha$
and the scale parameter $s=1/\beta$. Therefore, the estimation of the LL
parameters can be viewed as an estimation problem of $\mu$ and $s$ of the
logistic distribution based on $n$ observations $\\{z_{i}:~{}z_{i}=\log
t_{i},~{}t_{i}>0,~{}i=1,2,\ldots,n\\}$.

Since the parameter $\mu$ is the median of the logistic distribution, the
sample median (SM) appears to be a natural choice for estimating $\mu$,
denoted by $\hat{\mu}=\mathop{\mathrm{med}}_{1\leq i\leq
n}\\{z_{1},\cdots,z_{n}\\}$. A simple robust estimator for the scale parameter
$s$ is the median absolute deviation (MAD) estimator given by

$\hat{s}=\dfrac{\underset{1\leq i\leq
n}{\mathrm{med}}\\{|z_{i}-\tilde{z}|\\}}{\Phi^{-1}(3/4)},$ (8)

where $\tilde{z}=\mathrm{med}\\{z_{1},\cdots,z_{n}\\}$, and $\Phi^{-1}(\cdot)$
is the inverse of the standard normal cdf. Here, $\Phi^{-1}(3/4)$ is needed to
make the estimator Fisher-consistent (Fisher 1922) for the standard deviation
under the normal distribution (Hampel et al. 1986, Park and Cho 2003). It
deserves mentioning that the MAD estimator in (8) can be easily obtained using
the mad($\cdot$) function from R language (Bellio and Ventura 2005). Then the
sample median and MAD estimators of the parameters of the LL distribution are
given by

$\hat{\alpha}_{\mathrm{SM}}=\exp\left(\hat{\mu}\right)\ \mathrm{and}\
\hat{\beta}_{\mathrm{MAD}}=\frac{1}{\hat{s}},$ (9)

respectively. Park et al. (2022) showed that that the asymptotic breakdown
points of the estimators in (9) are equal to $50\%$ .

###  2.4 Hodges-Lehmann and Shamos estimators

Within a location framework, Hodges and Lehmann (1963) showed that the Hodges-
Lehmann (HL) estimator performs better than the sample median for an
asymmetric distribution in terms of efficiency. We may thus employ the HL
estimator for estimating $\mu$ and it is defined as the median of $n(n-1)/2$
pairwise averages of $n$ observations given by

$\tilde{\mu}=\underset{i<j}{\mathrm{med}}\left\\{\dfrac{z_{i}+z_{j}}{2}\right\\},$
(10)

which can be easily obtained with the hl.loc($\cdot$) function in the R ICSNP
package (Nordhausen et al. 2018). Note that the HL estimator has a reasonable
breakdown point of $29.3\%$ (Hettmansperger and McKean 1998). We then adopt
the Fisher-consistent Shamos estimator (Shamos 1976, Lévy-Leduc et al. 2011)
for estimating the scale parameter $s$ and it is defined as

$\tilde{s}=\dfrac{\underset{i<j}{\mathrm{med}}\\{|z_{i}-z_{j}|\\}}{\sqrt{2}\Phi^{-1}(3/4)}.$
(11)

To compute $\tilde{s}$, we let $\mathbf{Z}$ be the column vector such that
$\mathbf{Z}=[z_{1},\cdots,z_{n}]^{T}$. Similar to the outer product, we define
the outer difference as below

$\mathbf{Z}\ominus\mathbf{Z}=\begin{bmatrix}z_{1}-z_{1}&z_{1}-z_{2}&\cdots&z_{1}-z_{n}\\\
z_{2}-z_{1}&z_{2}-z_{2}&\cdots&z_{2}-z_{n}\\\ \vdots&\vdots&\ddots&\vdots\\\
z_{n}-z_{1}&z_{n}-z_{2}&\cdots&z_{n}-z_{n}\\\ \end{bmatrix},$

which can be easily obtained using R language function such as outer(Z, Z,
"-"). After obtaining the above outer difference, we take only lower or upper
triangle matrix which is also obtained by using lower.tri() or upper.tri().
Then we have a set of $\mathcal{Z}=\\{z_{i}-z_{j}:i<j\\}$ and take the
absolute value of each element so that we have
$\mathcal{W}=\\{|z_{i}-z_{j}|:i<j\\}$. Then we take the median using all the
elements in $\mathcal{W}$. We may also use the shamos($\cdot$) function in
rQCC package to achieve the computation; see Park and Wang (2020).

It is worth noting that the Shamos estimator has the same breakdown point of
$29.3\%$ as the HL estimator (Rousseeuw and Croux 1993). Then the Hodges-
Lehmann and Shamos estimators of the parameters of the LL distribution are
given by

$\hat{\alpha}_{\mathrm{HL}}=\exp\left(\tilde{\mu}\right)\ \mathrm{and}\
\hat{\beta}_{\mathrm{Shamos}}=\frac{1}{\tilde{s}},$ (12)

respectively. It can be shown that the asymptotic breakdown points of the
estimators in (12) are also equal to $29.3\%$. We here refer the interested
reader to (Park et al., 2022) for a detailed derivation.

##  3 Simulation study

In this section, we carry out Monte Carlo simulations to investigate the
finite sample performance of each considered estimator for the case without
data contamination (Section 3.1) and with data contamination (Section 3.2).
For the simplicity of notation, we let “PE1”, “PE2” and “PE3” stand for the
symmetrical percentile estimators with the $5$th and $95$th, the $10$th and
$90$th, the $33$th and $67$th, respectively.

###  3.1 Numerical results without contamination

We generate random samples of size $n=10,~{}25,~{}50,~{}75,~{}100$ for the LL
distribution with $\alpha=1$ and $\beta=1.5,~{}2.5,~{}5,~{}10$ by using
rllogis($\cdot$) from the R package actuar (Dutang et al. 2008). Without loss
of generality, we may fix $\alpha$ at $1.0$, since it is just a scale
parameter. We replicate each simulation setting $M=10,000$ times and compute
the average bias and root mean square error (RMSE) of each estimator given by

$\displaystyle\mathrm{Bias}=\frac{1}{M}\sum_{i=1}^{M}\big{(}\hat{w}_{i}-w_{T}\big{)}$
$\displaystyle\mathrm{and}\quad\mathrm{RMSE}=\sqrt{\frac{1}{M}\sum_{i=1}^{M}\big{(}\hat{w}_{i}-w_{T}\big{)}^{2}},$

where $\hat{w}_{i}$ and $w_{T}$ stand for the estimate and true parameter for
$w_{T}=\\{\alpha,\beta\\}$, respectively.

Table 1: Average bias and RMSE (in parentheses) of estimates for the scale parameter $\alpha$ based on 10,000 simulations |  | Estimation methods
---|---|---
$n$ | $\beta$ |  MLE |  PE1 |  PE2 |  PE3 |  RM |  SM |  HL
10 | 1.5 |  $0.074$ |  $0.299$ |  $0.179$ |  $0.085$ |  $0.065$ |  $0.088$ |  $0.076$
(0.424) |  (0.951) |  (0.645) |  (0.453) |  (0.425) |  (0.472) |  (0.430)
2.5 |  $0.026$ |  $0.099$ |  $0.058$ |  $0.030$ |  $0.021$ |  $0.031$ |  $0.027$
(0.233) |  (0.381) |  (0.303) |  (0.247) |  (0.235) |  (0.257) |  (0.237)
5.0 |  $0.007$ |  $0.024$ |  $0.014$ |  $0.008$ |  $0.004$ |  $0.008$ |  $0.007$
(0.113) |  (0.160) |  (0.139) |  (0.119) |  (0.114) |  (0.123) |  (0.114)
| 10.0 |  $0.002$ |  $0.006$ |  $0.004$ |  $0.002$ |  $0.000$ |  $0.002$ |  $0.002$
|  (0.056) |  (0.077) |  (0.068) |  (0.059) |  (0.057) |  (0.061) |  (0.057)
25 | 1.5 |  $0.027$ |  $0.097$ |  $0.064$ |  $0.030$ |  $0.023$ |  $0.035$ |  $0.027$
(0.241) |  (0.450) |  (0.339) |  (0.256) |  (0.242) |  (0.280) |  (0.241)
2.5 |  $0.009$ |  $0.034$ |  $0.023$ |  $0.010$ |  $0.007$ |  $0.013$ |  $0.010$
(0.140) |  (0.241) |  (0.190) |  (0.149) |  (0.141) |  (0.162) |  (0.141)
5.0 |  $0.002$ |  $0.009$ |  $0.006$ |  $0.002$ |  $0.001$ |  $0.003$ |  $0.002$
(0.069) |  (0.115) |  (0.092) |  (0.073) |  (0.070) |  (0.080) |  (0.070)
| 10.0 |  $0.001$ |  $0.002$ |  $0.001$ |  $0.001$ |  $0.000$ |  $0.001$ |  $0.001$
|  (0.035) |  (0.057) |  (0.046) |  (0.037) |  (0.035) |  (0.040) |  (0.035)
50 | 1.5 |  $0.014$ |  $0.058$ |  $0.029$ |  $0.017$ |  $0.013$ |  $0.019$ |  $0.015$
(0.167) |  (0.314) |  (0.236) |  (0.178) |  (0.168) |  (0.192) |  (0.168)
2.5 |  $0.006$ |  $0.021$ |  $0.011$ |  $0.007$ |  $0.004$ |  $0.007$ |  $0.006$
(0.099) |  (0.177) |  (0.137) |  (0.105) |  (0.099) |  (0.113) |  (0.099)
5.0 |  $0.002$ |  $0.006$ |  $0.003$ |  $0.002$ |  $0.001$ |  $0.002$ |  $0.002$
(0.049) |  (0.086) |  (0.068) |  (0.052) |  (0.049) |  (0.056) |  (0.049)
| 10.0 |  $0.000$ |  $0.002$ |  $0.001$ |  $0.001$ |  $0.000$ |  $0.001$ |  $0.000$
|  (0.024) |  (0.043) |  (0.034) |  (0.026) |  (0.025) |  (0.028) |  (0.025)
75 | 1.5 |  $0.009$ |  $0.038$ |  $0.021$ |  $0.011$ |  $0.009$ |  $0.013$ |  $0.010$
(0.136) |  (0.253) |  (0.190) |  (0.144) |  (0.137) |  (0.158) |  (0.136)
2.5 |  $0.004$ |  $0.014$ |  $0.008$ |  $0.004$ |  $0.003$ |  $0.005$ |  $0.004$
(0.080) |  (0.146) |  (0.112) |  (0.085) |  (0.081) |  (0.093) |  (0.081)
5.0 |  $0.001$ |  $0.004$ |  $0.002$ |  $0.001$ |  $0.001$ |  $0.001$ |  $0.001$
(0.040) |  (0.072) |  (0.055) |  (0.042) |  (0.040) |  (0.046) |  (0.040)
| 10.0 |  $0.000$ |  $0.001$ |  $0.001$ |  $0.000$ |  $0.000$ |  $0.000$ |  $0.000$
|  (0.020) |  (0.036) |  (0.028) |  (0.021) |  (0.020) |  (0.023) |  (0.020)
100 | 1.5 |  $0.007$ |  $0.028$ |  $0.015$ |  $0.008$ |  $0.006$ |  $0.009$ |  $0.006$
(0.118) |  (0.225) |  (0.167) |  (0.124) |  (0.119) |  (0.136) |  (0.118)
2.5 |  $0.003$ |  $0.011$ |  $0.006$ |  $0.003$ |  $0.002$ |  $0.003$ |  $0.003$
(0.070) |  (0.131) |  (0.099) |  (0.074) |  (0.071) |  (0.081) |  (0.070)
5.0 |  $0.001$ |  $0.003$ |  $0.002$ |  $0.001$ |  $0.000$ |  $0.001$ |  $0.001$
(0.035) |  (0.065) |  (0.049) |  (0.037) |  (0.035) |  (0.040) |  (0.035)
| 10.0 |  $0.000$ |  $0.001$ |  $0.001$ |  $0.000$ |  $0.000$ |  $0.000$ |  $0.000$
|  (0.017) |  (0.032) |  (0.024) |  (0.018) |  (0.018) |  (0.020) |  (0.017)
Note: Bold value indicates the smallest values of RMSE for all the considered
estimators
Table 2: Average bias and RMSE (in parentheses) of estimates for the shape parameter $\beta$ based on 10,000 simulations |  | Estimation methods
---|---|---
$n$ | $\beta$ |  MLE |  PE1 |  PE2 |  PE3 |  RM |  MAD |  Shamos
10 | 1.5 |  $0.221$ |  $0.570$ |  $0.450$ |  $0.743$ |  $-0.004$ |  $-0.334$ |  $-0.596$
(0.578) |  (0.874) |  (0.808) |  (1.836) |  (0.555) |  (0.681) |  (0.669)
2.5 |  $0.368$ |  $0.936$ |  $0.779$ |  $1.245$ |  $-0.006$ |  $-0.556$ |  $-0.993$
(0.963) |  (1.440) |  (1.362) |  (3.065) |  (0.926) |  (1.136) |  (1.115)
5.0 |  $0.735$ |  $1.857$ |  $1.589$ |  $2.496$ |  $-0.012$ |  $-1.113$ |  $-1.986$
(1.926) |  (2.862) |  (2.739) |  (6.131) |  (1.851) |  (2.272) |  (2.231)
| 10.0 |  $1.470$ |  $3.705$ |  $3.196$ |  $4.996$ |  $-0.025$ |  $-2.226$ |  $-3.972$
|  (3.852) |  (5.716) |  (5.487) |  (12.258) |  (3.702) |  (4.543) |  (4.462)
25 | 1.5 |  $0.082$ |  $0.254$ |  $0.179$ |  $0.231$ |  $-0.015$ |  $-0.497$ |  $-0.614$
(0.291) |  (0.440) |  (0.382) |  (0.640) |  (0.295) |  (0.567) |  (0.636)
2.5 |  $0.136$ |  $0.421$ |  $0.296$ |  $0.386$ |  $-0.025$ |  $-0.828$ |  $-1.023$
(0.485) |  (0.733) |  (0.636) |  (1.067) |  (0.492) |  (0.945) |  (1.059)
5.0 |  $0.272$ |  $0.839$ |  $0.592$ |  $0.772$ |  $-0.050$ |  $-1.656$ |  $-2.046$
(0.970) |  (1.465) |  (1.271) |  (2.134) |  (0.984) |  (1.891) |  (2.119)
| 10.0 |  $0.544$ |  $1.676$ |  $1.182$ |  $1.544$ |  $-0.101$ |  $-3.313$ |  $-4.092$
|  (1.941) |  (2.930) |  (2.541) |  (4.268) |  (1.968) |  (3.782) |  (4.238)
50 | 1.5 |  $0.040$ |  $0.123$ |  $0.092$ |  $0.109$ |  $-0.009$ |  $-0.542$ |  $-0.619$
(0.191) |  (0.265) |  (0.248) |  (0.384) |  (0.204) |  (0.568) |  (0.629)
2.5 |  $0.066$ |  $0.204$ |  $0.154$ |  $0.181$ |  $-0.016$ |  $-0.903$ |  $-1.031$
(0.318) |  (0.441) |  (0.413) |  (0.640) |  (0.339) |  (0.946) |  (1.048)
5.0 |  $0.132$ |  $0.408$ |  $0.308$ |  $0.363$ |  $-0.031$ |  $-1.805$ |  $-2.063$
(0.636) |  (0.882) |  (0.826) |  (1.281) |  (0.679) |  (1.893) |  (2.096)
| 10.0 |  $0.264$ |  $0.815$ |  $0.616$ |  $0.725$ |  $-0.063$ |  $-3.611$ |  $-4.125$
|  (1.272) |  (1.765) |  (1.653) |  (2.561) |  (1.357) |  (3.785) |  (4.192)
75 | 1.5 |  $0.027$ |  $0.083$ |  $0.062$ |  $0.071$ |  $-0.007$ |  $-0.554$ |  $-0.620$
(0.151) |  (0.206) |  (0.192) |  (0.291) |  (0.163) |  (0.570) |  (0.626)
2.5 |  $0.045$ |  $0.139$ |  $0.103$ |  $0.119$ |  $-0.011$ |  $-0.923$ |  $-1.033$
(0.252) |  (0.344) |  (0.319) |  (0.484) |  (0.272) |  (0.951) |  (1.044)
5.0 |  $0.090$ |  $0.279$ |  $0.205$ |  $0.238$ |  $-0.022$ |  $-1.847$ |  $-2.067$
(0.503) |  (0.688) |  (0.639) |  (0.968) |  (0.543) |  (1.901) |  (2.088)
| 10.0 |  $0.179$ |  $0.5572$ |  $0.410$ |  $0.476$ |  $-0.045$ |  $-3.693$ |  $-4.134$
|  (1.007) |  (1.376) |  (1.277) |  (1.936) |  (1.086) |  (3.803) |  (4.177)
100 | 1.5 |  $0.020$ |  $0.065$ |  $0.047$ |  $0.053$ |  $-0.005$ |  $-0.561$ |  $-0.621$
(0.129) |  (0.176) |  (0.164) |  (0.246) |  (0.141) |  (0.573) |  (0.626)
2.5 |  $0.034$ |  $0.108$ |  $0.079$ |  $0.088$ |  $-0.008$ |  $-0.935$ |  $-1.035$
(0.215) |  (0.294) |  (0.273) |  (0.410) |  (0.235) |  (0.955) |  (1.043)
5.0 |  $0.068$ |  $0.217$ |  $0.157$ |  $0.176$ |  $-0.016$ |  $-1.871$ |  $-2.070$
(0.430) |  (0.587) |  (0.547) |  (0.821) |  (0.470) |  (1.910) |  (2.086)
| 10.0 |  $0.136$ |  $0.434$ |  $0.314$ |  $0.351$ |  $-0.033$ |  $-3.741$ |  $-4.140$
|  (0.860) |  (1.174) |  (1.093) |  (1.641) |  (0.939) |  (3.820) |  (4.171)
Note: Bold value indicates the smallest values of RMSE for all the considered
estimators

Numerical results are provided in Tables 1 and 2. Several findings for the
case without data contamination can be summarized as follows. (1) We observe
that in terms of RMSE, the MLE of $\alpha$ outperforms other estimators most
often and that in terms of bias, the RM estimator of $\alpha$ is superior than
other estimators in most cases. (2) When the sample size is small (e.g.,
$n=10$), the RM estimator of $\beta$ provides better performance than other
estimators in terms of bias and RMSE, and the performance of the MLE for
$\beta$ becomes superior in terms of RMSE, as $n$ increases. (3) The
considered percentile estimators PEs of $\alpha$ usually perform worse than
other estimators in terms of bias and RMSE, whereas they perform well for
estimating of $\beta$. (4) As the sample size $n$ becomes larger, the bias and
RMSE of all the considered estimators decrease significantly, and all the
estimators behave similarly.

###  3.2 Numerical results with contamination

To study the robustness of the proposed estimators, we first generate $n=25$
LL random variables with $\alpha=1.0$, $\beta=10$ as the reference
distribution (no contamination case) and then consider the case with data
contamination under the following four scenarios: (i) $10\%$ replacement
outliers from the LL distribution with $\alpha=1.0$ and $\beta=0.1$; (ii)
$10\%$ replacement outliers from the LL distribution with $\alpha=4.0$ and
$\beta=10$; (iii) $10\%$ replacement outliers from a uniform distribution on
the interval $(0,~{}20)$, and (iv) $10\%$ extreme contaminated values from a
point mass distribution at the value of $50$.

Table 3: Average bias and RMSE (in the parentheses) of the scale ($\alpha$) and shape ($\beta$) estimators with $n=25$ under the five different simulation settings based on 10,000 simulations | MLE | PE1 | PE2 | PE3 | RM | SM/MAD | HL/Shamos
---|---|---|---|---|---|---|---
No contamination from LL$(\alpha=1.0$, $\beta=10)$ |
$\alpha$ | $0.001$ | $0.002$ | $0.001$ | $0.001$ | $0.000$ | $0.001$ | $0.001$
(0.035) | (0.057) | (0.046) | (0.037) |  (0.035) | (0.040) |  (0.035)
$\beta$ | $0.544$ | $1.676$ | $1.182$ | $1.544$ | $-0.101$ | $-3.313$ | $-4.092$
(1.941) | (2.930) | (2.541) | (4.268) | (1.968) | (3.782) | (4.238)
$10\%$ contamination from LL$(\alpha=1.0$, $\beta=0.1)$ |
$\alpha$ | $-0.168$ | $-0.788$ | $-0.309$ | $0.001$ | $-0.006$ | $0.001$ | $-0.005$
(0.170) | (0.788) | (0.310) | (0.037) |  (0.035) | (0.040) |  (0.035)
$\beta$ | $-8.836$ | $-8.377$ | $-6.160$ | $1.544$ | $-1.602$ | $-3.316$ | $-4.821$
(8.836) | (8.377) | (6.164) | (4.268) |  (2.365) | (3.784) | (4.932)
$10\%$ contamination from LL$(\alpha=4.0$, $\beta=10)$ |
$\alpha$ | $0.101$ | $0.804$ | $0.535$ | $0.058$ | $0.005$ | $0.051$ | $0.073$
(0.108) | (0.805) | (0.536) | (0.071) |  (0.036) | (0.067) | (0.083)
$\beta$ | $-4.824$ | $-5.995$ | $-5.963$ | $0.685$ | $-1.536$ | $-3.668$ | $-4.682$
(4.834) | (5.996) | (5.966) | (3.746) |  (2.330) | (4.061) | (4.800)
$10\%$ contamination from U$(0,20)$ |
$\alpha$ | $0.142$ | $2.191$ | $1.548$ | $0.058$ | $0.005$ | $0.051$ | $0.073$
(0.148) | (2.192) | (1.549) | (0.071) |  (0.036) | (0.067) | (0.083)
$\beta$ | $-6.959$ | $-7.745$ | $-7.912$ | $0.685$ | $-1.602$ | $-3.668$ | $-4.774$
(6.959) | (7.745) | (7.912) | (3.746) |  (2.365) | (4.061) | (4.888)
$10\%$ contamination from a point mass distribution at $50$ |
$\alpha$ | $0.104$ | $5.208$ | $3.988$ | $0.001$ | $0.005$ | $0.001$ | $0.006$
(0.111) | (5.213) | (3.991) | (0.037) |  (0.036) | (0.040) |  (0.036)
$\beta$ | $-7.965$ | $-8.589$ | $-8.786$ | $1.544$ | $-0.527$ | $-3.316$ | $-4.678$
(7.965) | (8.589) | (8.786) | (4.268) |  (2.081) | (3.784) | (4.797)
Note: Bold value indicates the smallest values of RMSE for all the considered
estimators

We conduct $M=10,000$ replications for each scenario above and calculate the
average bias and RMSE of all the estimators, which are summarized in Table 3.
Some conclusions for the case with data contamination can be drawn as follows.
(1) We observe that data contamination by outliers induces a large influence
on the MLE in terms of larger values of bias and RMSE. (2) The PE3 is superior
than PE1 and PE2, and this occurs because the observations on PE1 and PE2 are
outliers based on our simulation setups above. (3) The SM and HL estimators
seem performing well for estimating $\alpha$, whereas they are badly affected
by outliers when estimating $\beta$ for all the considered cases. (4) The RM
estimator provides the best performance among all the considered estimators in
most cases in terms of bias and RMSE. Finally, it deserves mentioning that we
have conducted more simulation studies with respect to other values of the
parameters and different sample sizes and obtained similar conclusions, which
are thus not presented here for simplicity.

As a result, we observe that when the data contamination is present, the MLEs
of the parameters are adversely affected by outliers. On the other hand, the
proposed estimators are quite robust to a certain level of data contamination
in terms of bias and RMSE, even when the sample size is small. Among the
proposed robust estimators, numerical results showed that the RM estimators
consistently provide more reliable results than others in most cases under
consideration. Therefore, we have a preference for the RM estimators as
alternative estimators to the MLEs for the parameters of the LL distribution,
especially when the data contain outliers.

##  4 A real-data application

In this section, we employ a real-data example to investigate the
effectiveness of the proposed robust estimators. The data are originated from
Nelson (1982) and are recently adopted by Abbas and Tang (2016) and Reath et
al. (2018) to illustrate the practical application of the LL distribution for
this modeling lifetime data. The data record the breakdown time in minutes of
an insulating fluid between electrodes at a voltage of $34$ kV. The data are
provided in Table 4. It can be seen from the boxplot of the data in Figure
1(a) that there exists an outlier in this data.

Table 4: The breakdown time in minutes of an insulating fluid data in Nelson (1982) 0.19 | 0.78 | 0.96 | 1.31 | 2.78 | 3.16 | 4.15 | 4.67 | 4.85 | 6.50
---|---|---|---|---|---|---|---|---|---
7.35 | 8.01 | 8.27 | 12.06 | 31.75 | 32.52 | 33.91 | 36.71 | 72.89 |

The parameter estimates of $\alpha$ and $\beta$ and the Kolmogorov-Smirnov
(KS) goodness of fit test (Smirnov 1948) statistics are provided in Table 5.
Here, the KS statistic is a nonparametric test that is often used to compare
different estimation methods in terms of the quality of fit: the smaller the
value of this statistic, the better the fit to the data. We observe from Table
5 that according to the KS test statistics, the best fit is obtained based on
the RM estimators. Furthermore, we depict the histogram of the data with the
fitted densities based on the MLE, PE3, RM, SM/MAD, and HL/Shamos estimation
methods in Figure 1(b). We can visualize from this figure that the best fit of
the data is also achieved by the RM method. Finally, the quantile-quantile
(Q-Q) plots of the fitted LL distributions based on the considered estimators
are presented in Figure 2. We observe that the quantile lines based on the
MLE, SM/MAD, and HL/Shamos estimators show a tendency towards the outlier,
whereas the RM estimator was not impacted. This real-data example further
demonstrates the superior performance of the RM estimators over other
considered estimators when the data contain outliers.

Table 5: Estimates of the parameters of the LL distribution for the breakdown time data Estimation method |  $\hat{\alpha}$ |  $\hat{\beta}$ |  KS ($p$-value)
---|---|---|---
MLE |  $6.2573$ |  $1.1732$ |  $0.1337~{}({0.8428})$
PE3 |  $5.8957$ |  $1.9374$ |  $0.2263~{}({0.2453})$
RM |  $6.0318$ |  $1.0153$ |  $\mathbf{0.1069}~{}({0.9654})$
SM/MAD |  $6.5000$ |  $0.7941$ |  $0.1492~{}({0.7372})$
HL/Shamos |  $6.0429$ |  $0.6014$ |  $0.1999~{}({0.3823})$
Note: The best value is marked in bold
Figure 1: (a) Box plot of the breakdown time data in Nelson (1982); (b)
Histogram with the fitted densities based on the considered estimators Figure
2: The Q-Q plots for the breakdown time data (Nelson 1982) estimated by the
considered estimation methods

##  5 Concluding remarks

In this paper, we have developed several alternative estimators to the MLEs
for the parameters of the LL distribution, which not only have simple closed-
form expressions, but also are quite robust to a certain level of data
contamination by outliers. In addition, we have investigated the breakdown
point of each estimator under consideration. Numerical results from simulation
studies showed that in terms of bias and RMSE, the proposed estimators perform
favorably in a manner that they are comparable to the MLEs in the absence of
data contamination and that they provide superior performance in the presence
of data contamination.

A natural question appears: which of the proposed robust estimators, the
percentile estimator, the repeated median, the sample median and median
absolute deviation estimators, or the Hodges-Lehmann and Shamos estimators
should be recommended for estimating the parameters of the LL distribution
when the data are contaminated by outliers. We have a preference for the RM
estimators, because they consistently offer a better performance than other
estimators under different simulation scenarios and also have a high breakdown
point of $50\%$. Finally, a real-data application further illustrated that the
RM estimators provide a better fit than other estimators in terms of the KS
test statistic and the Q-Q plots.

## Acknowledgments

The authors would like to acknowledge the comments and suggestions from two
anonymous reviewers, which have improved the quality of the manuscript. This
work was based on the first author’s dissertation research which was
supervised by the corresponding author. The work of Dr. Park was supported by
the National Research Foundation of Korea (NRF) grant funded by the Korean
government (NRF-2017R1A2B4004169).

## References

  * Abbas and Tang (2016) Abbas, K. and Tang, Y. (2016).  Objective bayesian analysis for log-logistic distribution.  Communications in Statistics-Simulation and Computation 45, 2782–2791.
  * Ashkar and Mahdi (2003) Ashkar, F. and Mahdi, S. (2003).  Comparison of two fitting methods for the log-logistic distribution.  Water Resources Research 39, 1217.
  * Balakrishnan and Malik (1987) Balakrishnan, N. and Malik, H. J. (1987).  Moments of order statistics from truncated log-logistic distribution.  Journal of Statistical Planning and Inference 17, 251–267.
  * Bellio and Ventura (2005) Bellio, R. and Ventura, L. (2005).  An introduction to robust estimation with r functions.  Proceedings of 1st International Work 1, 1–57.
  * Delignette-Muller and Dutang (2015) Delignette-Muller, M. L. and Dutang, C. (2015).  fitdistrplus: An R package for fitting distributions.  Journal of Statistical Software 64, 1–34.
  * Dutang et al. (2008) Dutang, C., Goulet, V., Pigeon, M. et al. (2008).  actuar: An R package for actuarial science.  Journal of Statistical software 25, 1–37.
  * Fisher (1922) Fisher, R. A. (1922).  On the mathematical foundations of theoretical statistics.  Philosophical Transactions of the Royal Society of London. Series A, Containing Papers of a Mathematical or Physical Character 222, 309–368.
  * Fisk (1961) Fisk, P. R. (1961).  The graduation of income distributions.  Econometrica: Journal of the Econometric Society 22, 171–185.
  * Geskus (2001) Geskus, R. B. (2001).  Methods for estimating the aids incubation time distribution when date of seroconversion is censored.  Statistics in Medicine 20, 795–812.
  * Hampel et al. (1986) Hampel, F. R., Ronchetti, E. M., Rousseeuw, P. and Stahel, W. A. (1986).  Robust Statistics: the Approach based on Influence Functions.  John Wiley & Sons, New York.
  * He et al. (2020a) He, X., Chen, W. and Qian, W. (2020a).  Maximum likelihood estimators of the parameters of the log-logistic distribution.  Statistical Papers 61, 1875–1892.
  * He et al. (2020b) He, X., Chen, W. and Yang, R. (2020b).  Modified best linear unbiased estimator of the shape parameter of log-logistic distribution.  Journal of Statistical Computation and Simulation In Press, 1–13.
  * Hettmansperger and McKean (1998) Hettmansperger, T. P. and McKean, J. W. (1998).  Robust Nonparametric Statistical Methods, vol. 5 of Kendall’s Library of Statistics.  Edward Arnold, London; copublished by John Wiley & Sons, Inc., New York.
  * Hodges and Lehmann (1963) Hodges, J. L. and Lehmann, E. L. (1963).  Estimates of location based on rank tests.  The Annals of Mathematical Statistics 34, 598–611.
  * Huber and Ronchetti (2009) Huber, P. J. and Ronchetti, E. M. (2009).  Robust Statistics.  John Wiley & Sons, New York, 2nd ed.
  * Kantam and Srinivasa (2002) Kantam, R. and Srinivasa, G. (2002).  Log-logistic distribution: Modified maximum likelihood estimation.  Gujarat Statistical Review 29, 25–36.
  * Lévy-Leduc et al. (2011) Lévy-Leduc, C., Boistard, H., Moulines, E., Taqqu, M. and Reisen, V. (2011).  Large sample behaviour of some well-known robust estimators under long-range dependence.  Statistics 45, 59–71.
  * Nelson (1982) Nelson, W. (1982).  Applied Life Data Analysis.  John Wiley & Sons, Inc., New York.  Wiley Series in Probability and Mathematical Statistics.
  * Nordhausen et al. (2018) Nordhausen, K., Sirkia, S., Oja, H. and Tyler, D. E. (2018).  ICSNP: Tools for Multivariate Nonparametrics.  R package version 1.1-1.
  * Park and Cho (2003) Park, C. and Cho, B. R. (2003).  Development of robust design under contaminated and non-normal data.  Quality Engineering 15, 463–469.
  * Park et al. (2022) Park, C., Kim, H. and Wang, M. (2022).  Investigation of finite-sample properties of robust location and scale estimators.  Communication in Statistics – Simulation and Computation 51, 2619–2645.
  * Park and Wang (2020) Park, C. and Wang, M. (2020).  rQCC: Robust Quality Control Chart R.  Package Version 1.21.0 (published in January 0, 2021).
  * Park et al. (2020) Park, C., Wang, M. and Hwang, W. Y. (2020).  A study on robustness of the paired sample tests.  Industrial Engineering & Management Systems, 19, 386–397.
  * R Core Team (2018) R Core Team (2018).  R: A Language and Environment for Statistical Computing.  R Foundation for Statistical Computing, Vienna, Austria.
  * Reath et al. (2018) Reath, J., Dong, J. and Wang, M. (2018).  Improved parameter estimation of the log-logistic distribution with applications.  Computational Statistics 33, 339–356.
  * Ross (2000) Ross, S. M. (2000).  Introduction to Probability and Statistics for Engineers and Scientists.  Harcourt/Academic Press, Burlington, MA, 2nd ed.
  * Rousseeuw and Croux (1993) Rousseeuw, P. J. and Croux, C. (1993).  Alternatives to the median absolute deviation.  Journal of the American Statistical association 88, 1273–1283.
  * Shamos (1976) Shamos, M. (1976).  Geometry and Statistics: Problems at the Interface.  New York: Academic Press, J.F. Traub ed.
  * Shoukri et al. (1988) Shoukri, M. M., Mian, I. U. H. and Tracy, D. S. (1988).  Sampling properties of estimators of the log-logistic distribution with application to Canadian precipitation data.  The Canadian Journal of Statistics 16, 223–236.
  * Siegel (1982) Siegel, A. F. (1982).  Robust regression using repeated medians.  Biometrika 69, 242–244.
  * Smirnov (1948) Smirnov, N. (1948).  Table for estimating the goodness of fit of empirical distributions.  The Annals of Mathematical Statistics 19, 279–281.

## Appendix A

Proof of Proposition 1: It is worth noting that the estimator of the shape
parameter $\hat{\beta_{p}}$ in (4) tends to infinity (explosion) when
$t_{h}=t_{l}$. This could happen if a proportion ($h-l$) of the observations
is placed to the same position as $t_{l}$. Also, the estimator
$\hat{\beta_{p}}$ goes to zero (implosion) if $t_{h}$ approaches infinity and
$t_{l}$ remains bounded or if $t_{l}$ goes to zero and $t_{h}$ remains
bounded. In this case, it suffices to place ($1-h$) observations to go to
infinity or $l$ observations to be zero. These observations indicate that the
breakdown point of $\hat{\beta_{p}}$ is given by

$\displaystyle\kappa=\min\\{h,~{}1-l,~{}h-l,~{}1-h\\}=\min\\{1-l,~{}h-l,~{}1-h\\}.$

If we consider the symmetric percentiles (i.e. $h=1-l$) and $h>1/2$, the
breakdown point of the symmetrical percentile estimator $\hat{\beta_{p}}$ is
then given by

$\kappa=\min\\{2h-1,~{}1-h\\}.$

The highest breakdown point for the symmetrical percentile estimator
$\hat{\beta_{p}}$ is obtained at the intersection of lines $2h-1=1-h$ for
$h=2/3$ and $l=1/3$. In a similar way as done for the shape estimator
$\hat{\beta_{p}}$, we can show that the breakdown point of the scale estimator
$\hat{\alpha_{p}}$ is equal to $\min\\{h,1-l\\}$. Therefore, we obtain the
highest asymptotic breakdown point the percentile estimator $\hat{\alpha}_{p}$
corresponding to the symmetrical percentile estimator $\hat{\beta}_{p}=2/3$
when $h=1-l=2/3$. This completes the proof.

## Appendix B

R code for the real-data application in Section 4 :

    library("fitdistrplus") # mledist
    library("rQCC")
    set.seed(1)
    # The time to breakdown of an insulating fluid
    X <- c(0.19, 0.78, 0.96, 1.31, 2.78, 3.16, 4.15, 4.67, 4.85, 6.50,
           7.35, 8.01, 8.27, 12.06, 31.75, 32.52, 33.91, 36.71, 72.89)
    n <- length(X)
    x <- log(X)             # rewrite log-logistic to location-scale form
    F <- ppoints(n, a=0)    # F = i/(n+1)
    y <- log((1-F)^(-1)-1)  # rewrite log-logistic to linear regression model

    # Maximum likelihood estimators
    dLogL <- function(x, beta, alpha)(beta/alpha)*
             (x/alpha)^(beta-1)*((1+(x/alpha)^beta)^-2)  # log-logistic pdf
    pLogL <- function(x, beta, alpha)
             (x^beta/(x^beta + alpha^beta))  # log-logistic cdf
    fit.mle <- mledist(X, "LogL", start = list(alpha = 0.2, beta = 0.01))
    fit.mle$estimate

    # Percentile estimators
    H <- 0.67
    L <- 0.33
    betahat.percentile <- function(data, H, L){
      PH <- quantile(data, H)
      PL <- quantile(data, L)
      betahat.percentile <- 2*(log(H) - log(L))*(log(PH/PL))^(-1)
      return(betahat.percentile)
    }
    betahat.pct <- betahat.percentile(X, H, L)
    alphahat.pct<- exp(log(quantile(X, H)) - 1/betahat.pct * log(H/L))

    # Repeated median estimators
    beta1hat.med <- function(x, y) {
      int_list <- final_list <- numeric(0)
      for (j in 1:n) {
        for (i in 1:n) {
          if (i != j) {
            int_list = append(int_list, (y[j] - y[i]) / (x[j] - x[i]))
          }
        }
        final_list <- append(final_list, median(int_list))
      }
      return(median(final_list))
    }
    beta1hat.RM <- beta1hat.med(x, y)
    beta0hat.RM <- median(y - beta1hat.RM*x)
    betahat.rm <- mean(beta1hat.RM)
    alphahat.rm<- exp(-beta0hat.RM/beta1hat.RM)

    # Median absolute deviation (MAD) & Shamos estimators for alpha (scale)
    # && sample median (SM) & Hodges-Lehmann (HL) estimators for beta (shape)
    shape_est <- function(data){
      alpha.mad <- 1/mad(data)
      alpha.shamos <- 1/shamos(data)
      shape_est.result <- cbind(alpha.mad, alpha.shamos)
      return(shape_est.result)
    }
    betahat.mad <- shape_est(x)[,1]
    betahat.shamos <- shape_est(x)[,2]

    scale_est <- function(data){
      beta.median <- exp(median(data))
      beta.hl <- exp(HL(data, "HL2"))
      scale_est.result <- cbind(beta.median, beta.hl)
      return(scale_est.result)
    }
    alphahat.median <- scale_est(x)[,1]
    alphahat.hl <- scale_est(x)[,2]


# The expected value of sample information calculations for external
validation of risk prediction models

Mohsen Sadatsafavi  Andrew J Vickers  Tae Yoon Lee  Paul Gustafson  Laure
Wynants

(January 04, 2024)

###### Abstract

In designing external validation studies of clinical prediction models,
contemporary sample size calculation methods are based on the frequentist
inferential paradigm. One of the widely reported metrics of model performance
is net benefit (NB), and the relevance of conventional inference around NB as
a measure of clinical utility is doubtful. Value of Information methodology
quantifies the consequences of uncertainty in terms of its impact on clinical
utility of decisions. We introduce the expected value of sample information
(EVSI) for validation as the expected gain in NB from conducting an external
validation study of a given size. We propose algorithms for EVSI computation,
and in a case study demonstrate how EVSI changes as a function of the amount
of current information and future study’s sample size. Value of Information
methodology provides a decision-theoretic lens to the process of planning a
validation study of a risk prediction model and can complement conventional
methods when designing such studies.

††footnotetext: From Faculty of Medicine and Pharmaceutical Sciences (MS,
TYL), and Department of Statistics (PG), The University of British Columbia,
Vancouver, British Columbia, Canada; Department of Epidemiology and
Biostatistics, Memorial Sloan Kettering Cancer Center, New York, New York, USA
(AV); Department of Epidemiology, CAPHRI Care and Public Health Research
Institute, Maastricht University, Maastricht, The Netherlands, and Department
of Development and Regeneration, KU Leuven, Leuven, Belgium
(LW)††footnotetext: *Correspondence to Mohsen Sadatsafavi, 2405 Wesbrook Mall,
Vancouver, BC, V6T1Z3, Canada; email<EMAIL_ADDRESS>
## Introduction

Clinical prediction models that quantify the risk of medically significant
events are a key aspect of individualized medicine. Before being trusted for
clinical use in a target population, a model needs to undergo validation in a
representative sample from that population. Model performance during such
external validation is assessed in various ways, which can broadly be
classified as prediction error (e.g., Brier score), discrimination (e.g.,
c-statistic), calibration (e.g., calibration intercept and slope), and net
benefit (NB)1. Among such metrics, NB is a decision-theoretic one, as it
enables direct assessment of the clinical utility of a model (i.e., if a
model’s expected NB is higher than alternatives, it is expected to confer
clinical utility). Due to such strong decision-theoretic underpinnings, and
despite its relatively new arrival in the field of predictive analytics, NB
has gained significant momentum and has become a standard component of modern
external validation studies2.

When interpreting the results of a validation study, the finite size of the
validation sample means that the assessment of model performance is
accompanied by uncertainty. Uncertainty in conventional metrics of model
performance is typically communicated using classical inferential methods
(e.g., 95% confidence interval [CI] around c-statistic or calibration slope,
p-values for comparing the c-statistic of competing models). However, the
relevance of such frequentist approach to uncertainty around NB as a decision-
theoretic metric is doubtful3–5.

Decision theory provides an alternative view to the consequences of
uncertainty by relating it to the outcome of decisions. Leaving aside, for the
sake of simplicity, the issue of risk aversion, the decision-maker should
adopt the strategy with the highest expected utility, which, if implementation
costs are comparable, is the one with the highest net benefit. While
uncertainty around NB should not affect the ‘adoption decision’, it is
associated with utility loss as it hinders our ability to identify the
decision with the highest true NB. The extent of uncertainty should thus
inform the ‘research decision’: whether the expected utility loss is large
enough to necessitate collecting further evidence. This approach towards
uncertainty quantification is referred to as Value of Information (VoI)
analysis6–8.

VoI methods are widely accepted in health policy-making7,9. However, their
relevance to clinical decision-making and risk prediction has only recently
been highlighted10. In particular, the expected value of perfect information
(EVPI), the expected gain in NB by completely eliminating uncertainty, has
been applied to both the development and validation phases of risk prediction
models10,11. If based on EVPI, the investigator concludes that uncertainty is
large and procuring more validation data is warranted, the next natural
question is ‘how much data should be collected?’. The expected value of sample
information (EVSI) is the expected gain in NB by conducting a study of a given
sample size (and design)12. In the context of model validation, EVSI can thus
inform the sample size determination of future validation studies.

The aim of the present work is to define EVSI for the external validation
phase of a risk prediction model and propose algorithms for its computation.
The rest of this manuscript is structured as follows. After outlining the
context, we define EVSI for a future study aimed at evaluating the clinical
utility of a model in a target population. We propose algorithms for EVSI
computation based on bootstrapping when preliminary data from the target
populations are available, a beta-binomial modeling approach for binary
outcomes when current information can be specified as independent Beta
distributions, and a more general algorithm when the current information is
expressed as an arbitrary probability distribution. A case study puts the
developments in context where we provide practical suggestions on reporting
and interpreting EVSI. We conclude by proposing areas of further inquiry.

## Methods

### Context and notations

We focus on validating a previously-developed risk prediction model in a well-
defined new target population. The developments are presented for models for
binary responses but are generally applicable to any context where NB can be
assessed (e.g., survival outcomes13, models for treatment benefit14). The
model is advertised as a function that maps patient characteristics to an
estimate of event risk. We acknowledge the relevance of other sources of
uncertainty (e.g., missing data, the plausibility that the sample is truly
representative of the target population) but focus the developments on
sampling uncertainty, and discuss some other uncertainty sources whenever
relevant.

VoI analysis is a Bayesian paradigm and requires explicit specification of
current information on the performance of the model in the target population.
By default, we assume this information is generated from a preliminary
validation exercise based on an iid sample $\mathbf{d}$ of $n$ individuals.
However, such information can be generated in other ways. For example, one can
extrapolate (optimism-corrected) performance of a model from its development
sample, or its performance in previous validation studies15, and accommodate
any added uncertainty about the compatibility of the populations (e.g. by
quantitatively discounting external information16). Formal expert elicitation
approaches can also be employed to construct distributions around model
performance metrics17. In such instances, $\mathbf{d}$ can be an abstract
entity representing, for example, the expert-elicited information.

Given our current information, we are considering the utility of procuring
data of sample of size $n^{*}$ for a future validation study in the same
target population. We are interested in expected utility gain from such a
study. We express utility in terms of net benefit (NB), as is quantified and
presented in decision curve analysis18. In a nutshell, NB is based on the
specification of a risk threshold, denoted by $z$, placed on predicted risks
such that individuals with predicted risk lower than this threshold are
considered low-risk, and those with predicted risk higher than this threshold
are classified as high-risk (and therefore are selected for treatment). For
individuals with predicted risk precisely equal to $z$, the decision-maker
will be ambivalent about the treatment decision. Because among these
individuals, treating all means a proportion $z$ is properly treated (true
positive) while the proportion $1-z$ is unnecessarily treated (false
positive), such ambivalence implies that each false positive has $-z/(1-z)$
times clinical utility of a true positive. Because the true and false positive
rates of the model can be estimated at this threshold, this implied exchange
rate can be used to evaluate the NB of a model in net true positive or false
positive units. Any given model should be compared with at least two ‘default’
strategies of treating no one and treating all. Treating no one has NB=0 by
definition, and NB calculations for treating all follow the same logic as
above, but by considering all individuals as positive. In practice, NB is
measured over a range of plausible risk thresholds.

Let $\theta$ be the set of parameters that fully specify the target population
for NB calculations. A common setup for binary outcomes is to define $\theta$
as the triplet of prevalence ($\theta_{p}$), sensitivity ($\theta_{se}$), and
specificity ($\theta_{sp}$)20. We index treating no one, model-based
treatment, and treating all by, respectively, subscripts 0,1,2, resulting in
the following equation for NBs (for brevity, we drop the notation that
indicates quantities are functions of $z$)

$NB(i,\mathbf{\theta})=\left\\{\begin{array}[]{lll}0&i=0&\mbox{(treat no
one)}\\\
\theta_{p}\theta_{se}-(1-\theta_{p})(1-\theta_{sp})\frac{z}{1-z}&i=1&(\mbox{use
model to decide})\\\ \theta_{p}-(1-\theta_{p})\frac{z}{1-z}&i=2&(\mbox{treat
all})\\\ \end{array}\right.,$ (1)

with $\theta=\\{\theta_{p},\theta_{se},\theta_{sp}\\}$. For time-to-event
outcomes, these quantities are time-dependent and should be defined at a time-
horizon of interest13. If there are alternative strategies they can also be
considered (the above indexing of $NB()$ can accommodate other strategies),
but this is omitted from this work for the ease of exposition.

Our current information about $\theta$ based on the available evidence is
expressed as $P(\theta|\mathbf{d})$. According to Bayes’ theorem,
$P(\theta|\mathbf{d})\propto P(\theta)P(\mathbf{d}|\theta)$, indicating that
this information is influenced by any prior information we might have had, as
well as by what we have learned from $\mathbf{d}$.

### The expected value of current information

With $P(\theta|\mathbf{d})$ quantifying their knowledge about $\theta$, a
risk-neutral decision-maker should pick the strategy with the highest expected
NB given the current information:

$\mathbb{E}NB_{\text{current
information}}=\max_{i}{\mathbb{E}_{\mathbb{\theta}|\mathbf{d}}NB(i,\mathbb{\theta})},$
(2)

where the expectation is with respect to $P(\theta|\mathbf{d})$. Generally,
evaluating this expectation might require a Monte Carlo approach (e.g., via
Markov Chain Monte Carlo sampling from $P(\theta|\mathbf{d})$). However, in
many instances, this expectation can be more readily obtained. For example, a
Bayesian interpretation of the bootstrap allows one to interpret an estimate
of $\theta$ obtained from a bootstrapped copy of $\mathbf{d}$ as a random draw
from the posterior distribution of $\theta$ given $\mathbf{d}$21. However, as
the number of bootstrap iterations increases, due to the law of large numbers
the average values from such bootstrap samples for each component of $\theta$
as well as $NB$s converge to their point estimate in the sample. This
justifies interpreting the point estimates of $NB$ from the sample, as is
conducted in conventional decision curve analysis21, in a Bayesian framework
as the posterior mean under a non-informative prior.

### The expected value of perfect information (EVPI)

Details of the reasoning and algorithms behind calculating EVPI is presented
previously11. In brief, if we could know the true values of $\theta$, we would
have picked the strategy with the highest true NB. We indeed do not have
access to the ground truth, but have partial information in terms of
$P(\theta|\mathbf{d})$. As such, we can model the consequence of knowing the
truth and take its expectation with respect to this distribution:

$\mathbb{E}NB_{\text{perfect
information}}=\mathbb{E}_{\mathbf{\theta}|\mathbf{d}}\max_{i}NB(i,\mathbf{\theta}).$
(3)

The EVPI is then the difference in expected NB of having perfect information
versus current information:

$EVPI=\mathbb{E}NB_{\text{perfect information}}-\mathbb{E}NB_{\text{current
information}}.$ (4)

The EVPI is a scalar quantity in NB units and quantifies the expected loss in
NB due to uncertainty, or alternatively, the expected gain in NB by completely
resolving uncertainty (i.e., conducting a validation study of infinite size).
We have proposed a generic algorithm based on bootstrapping to estimate this
expectation, as well as an asymptotic approach based on central limit theorem
and two-dimensional unit normal loss integral11.

### The expected value of sample information (EVSI)

The reasoning behind EVSI calculation is similar to that of EVPI, with the
modification that instead of knowing the truth, we will have more (but not
complete) information about $\theta$, with additional information coming from
a future iid sample $\mathcal{D^{*}}$ of size $n^{*}$. With this added
information we will have the opportunity to revise our current decision. We
will update our knowledge from $P(\theta|\mathbf{d})$ to
$P(\theta|\mathbf{d},\mathcal{D^{*}})$, and will select the strategy that has
the highest expected $NB$.

The NB of this approach once future data $\mathcal{D^{*}}$ are obtained is
$\max_{i}\mathbb{E}_{\mathbf{\theta|{\mathcal{D^{*}},d}}}NB(i,\mathbf{\theta})$.
We do not know what data we will observe in the future, but we can treat
$\mathcal{D^{*}}$ as a random variable, and take the expectation of this line
of reasoning with regard to its predictive distribution:

$\mathbb{E}NB_{\text{sample
information}}=\mathbb{E}_{\mathcal{D^{*}}|\mathbf{d}}\max_{i}\mathbb{E}_{\mathbf{\theta}|\mathcal{D^{*}},\mathbf{d}}NB(i,\mathbb{\theta}).$
(5)

The EVSI is the difference between this ENB and ENB under current information:

$EVSI=\mathbb{E}NB_{\text{sample information}}-\mathbb{E}NB_{\text{current
information}}.$ (6)

The EVSI is a non-negative scalar quantity in the same NB units as EVPI. It
quantifies the expected gain in NB by procuring a further validation sample of
a given size. The higher the EVSI, the higher the expected gain from the
planned validation study. EVPI, quantifying the expected NB gain by completely
resolving uncertainty, puts an upper limit on EVSI.

## EVSI computation algorithms

The challenge for EVSI calculation is the two nested expectations. The inner
expectation represents updated estimates of NBs after obtaining future data,
and the outer expectation is with respect to the predictive distribution of
future data. Below we detail three computaion algorithms.

### A bootstrap-based algorithm

If the current evidence $\mathbf{d}$ is in terms of a random sample of size
$n$ from the target population of interest, sampling from
$P(\theta|\mathbf{d})$ can be performed via the Bayesian bootstrap21. The
resulting algorithm for EVSI computation is similar to the previously proposed
method for EVSI calculations for data-driven economic evaluations22 and can be
seen as an extension of the bootstrap method for EVPI computation for external
validation studies11.

The Bayesian bootstrap of $\mathbf{d}$ assigns an improper non-informative
prior distribution $\text{Dirichlet}(0,0,...,0)$ across the $n$ data points. A
random draw from the posterior distribution of the generating population can
be materialized by drawing a random weight vector from the
$\text{Dirichlet}(1,1,...,1)$ distribution. Estimating $\theta$ in this
weighted sample is akin to drawing a value for $\theta$ from its posterior
distribution21. Further, once the generating population is at hand, the data
of the future study ($\mathcal{D^{*}}$) can be drawn by sampling from this
population. This results in a two-level resampling approach for EVSI
calculation explained in Table 1.

1. For $j$ = 1 to $N$ (number of Monte Carlo simulations) (a) Draw $\mathbf{F}$, a Bayesian bootstrapa of the data $\mathbf{d}$. This is a random draw from the posterior distribution of the population that has generated the sample. (b) Compute and record NBs from $\mathbf{F}$. (c) Record the maximum NB from those calculated in the previous step. (d) Draw $\mathcal{D^{*}}$, a random realization of the dataset of the future study, through sampling $n^{*}$ observations with replacement from $\mathbf{F}$. (e) Merge the future and existing data and estimate the expected NBs of each strategy in the pooled sample (this will be $\mathbb{E}_{\mathbf{\theta}|\mathcal{D^{*},d}}NB(i,\mathbf{\theta})$). Record maximum NB. Next $j$ 2. Take the average of NBs from step 1b and pick the maximum value. This is $\mathbb{E}NB_{\text{current information}}$b. 3. Average maximum NB from step 1c. This is $\mathbb{E}NB_{\text{perfect information}}$. From this subtract NB from step 2. This is EVPI. 4. Average maximum NB from step 1e. This is $\mathbb{E}NB_{\text{sample information}}$. From this subtract NB from step 2. This is EVSI. a Ordinary bootstrap can also be used - see text
b $\mathbb{E}NB_{\text{current information}}$ can also be directly evaluated
in the original sample without bootstrapping. However, we recommend averaging
over bootstrapped values as this approach is consistent with the way other
terms are evaluated and this prevents getting occasional negative VoI metrics

Table 1: Bootstrap-based computation of EVSI

The power of this approach is in the flexibility of the bootstrap method in
accommodating different types of outcomes and practical considerations in
validation studies. For example, if the dataset $\mathbf{d}$ has non-trivial
level of missing data, which affects the amount of information it contains,
the step 1b in the above algorithm can include an imputation component,
generating the population $\mathbf{F}$ without any missing values (and
similarly, step 1d can introduce missingness in the data to simulate a real-
world validation dataset).

The ordinary bootstrap, which assigns weights from a scaled
$\text{Multinomial}(1/n,1/n,...)$ distribution to observations, is
conceptually similar to the Bayesian bootstrap and has been interpreted in a
Bayesian way (e.g., in the approximate Bayesian bootstrap employed in missing
value imputation23). In data-driven VoI analyses it generates similar results
to the Bayesian bootstrap11,22. It thus can be used in lieu of the Bayesian
bootstrap in step 1a.

### An algorithm based on beta-binomial conjugacy for binary outcomes

In their proposal for a parametric Bayesian NB estimation for binary outcomes,
Cruz and Korthauer take advantage of the beta-binomial conjugacy and the fact
that the likelihood function for $NB$ can be factorized into terms involving
prevalence, sensitivity, and specificity. As such, they argue that specifying
prior information in terms of independent beta distributions on these
quantities will result in corresponding posterior distributions that are also
independent beta20. In this case, given that $P(\theta|\mathbf{d})$ is
composed of three independent beta distributions, by the same token
$P(\theta|\mathbf{d},\mathcal{D^{*}})$ after a realization of
$\mathcal{D^{*}}$ will also be the product of three independent beta
distributions. This results in a fast Monte-Carlo approach for EVSI
computation as explained in Table 2.

0. Specify the current information for NBs under current information as
$\theta_{p}\sim\text{Beta}(\alpha_{p},\beta_{p})$,
$\theta_{se}\sim\text{Beta}(\alpha_{se},\beta_{se})$,
$\theta_{sp}\sim\text{Beta}(\alpha_{sp},\beta_{sp})$. 1. Plug in the expected
values of $\theta$ from the above distributions in equation 1. Record the
maximum expected NB. This is $\mathbb{E}NB_{\text{current information}}$. 2.
For $j$ = 1 to $N$ (number of Monte Carlo simulations) (a) Randomly draw from
the distribution of $\theta$ as specified in step  ‣ 2. (b) Compute true NBs
from such draws. Record the maximum NB. (c) Generate $\mathcal{D^{*}}$,
defined by $\\{n^{*}_{tp},n^{*}_{fn},n^{*}_{tn},n^{*}_{fp}\\}$ given the draw
from $\theta$ obtained in step 2a as:
$n^{*}_{+}\sim\text{Binomial}(n^{*},\theta_{p})$ (number of positive cases in
the future sample),
$n^{*}_{tp}\sim\text{Binomial}(n^{*}_{+},\theta_{se})$,
$n^{*}_{fn}=n^{*}_{+}-n^{*}_{tp}$,
$n^{*}_{tn}\sim\text{Binomial}(n^{*}-n^{*}_{+},\theta_{sp})$,
$n^{*}_{fp}=n^{*}-n^{*}_{+}-n^{*}_{tn}$. (d) Calculate expected NBs given
generated $\mathcal{D^{*}}$ by plugging in the updated estimates of the
expected value of $\theta$ as follows:
${\bar{\theta}^{*}_{p}}=(\alpha_{p}+n^{*}_{tp}+n^{*}_{fn})/(\alpha_{p}+\beta_{p}+n^{*})$,
${\bar{\theta}^{*}_{se}}=(\alpha_{se}+n^{*}_{tp})/(\alpha_{se}+\beta_{se}+n^{*}_{tp}+n^{*}_{fn})$,
${\bar{\theta}^{*}_{sp}}=(\alpha_{sp}+n^{*}_{tn})/(\alpha_{sp}+\beta_{sp}+n^{*}_{tn}+n^{*}_{fp})$.
Record the maximum NB. This will be
$\mathbb{E}_{\mathbf{\theta|\mathbf{d},\mathcal{D^{*}}}}NB(i,\mathbf{\theta})$.
Next $j$ 3. Average NBs from step 2b. This is $\mathbb{E}NB_{\text{perfect
information}}$. From this subtract ENB from step 1. This is EVPI. 4. Average
NBs from step 2d. This is $\mathbb{E}NB_{\text{sample information}}$. From
this subtract ENB from step 1. This is EVSI.

Table 2: Computation of EVSI based on beta-binomial conjugacy

Independent beta distributions for each component of $\theta$ can emerge, from
example, if we ask an expert to express their belief about prevalence,
sensitivity, and specificity, in terms of fractions. That is, if their best
estimate of prevalence is 30%, and their uncertainty is akin to their opinion
comping from a sample of 10 individuals, the information can be expressed as
$\theta_{p}\sim\text{Beta}(3,7)$ (similar for sensitivity and specificity).

Importantly, for binary outcomes and in the absence of missing data, the
Bayesian bootstrap approach explained above reduces to this beta-binomial
model. This is because the information in a sample $\mathbf{d}$ can be fully
specified via four quantities: the number of true positives ($n_{tp}$), false
negatives ($n_{fn}$), true negatives ($n_{tn}$), and false positives
($n_{fp})$. The Bayesian bootstrapping of individual observations is akin to
assigning an improper prior $\text{Dirichlet}(0,0,0,0)$ to these quantities,
giving rise to the posterior distribution being of the form
$\text{Dirichlet}(n_{tp},n_{fn},n_{tn},n_{fp})$. Due to the aggregate property
of the Dirichlet distribution24, the $\text{Dirichlet}(0,0,0,0)$ prior of the
Bayesian bootstrap is equal to assigning $\text{Beta}(0,0)$ priors to
$\theta_{p}$, $\theta_{se}$, and $\theta_{sp}$. As such, the beta-binomial
algorithms will be equal to the Bayesian bootstrap and can be used instead for
VoI computations due to its computational efficiency.

Yet another utility of this approach is that instead of the improper prior
employed in the Bayesian bootstrap, the user can specify informative priors.
For example, if an independent cross-sectional study in the same target
population has estimated outcome prevalence that can be expressed as
$\text{Beta}(\alpha,\beta)$, then combined with observed data, our current
information about prevalence can be expressed as
$\theta_{p}\sim\text{Beta}(\alpha+n_{tp}+n_{fn},\beta+n_{fp}+n_{tn})$.
Further, even in the absence of prior information, proper non-informative or
weakly informative priors can be used to rectify some of the previously-
identified problems with bootstrap-based VoI calculations at extreme
situations. In small validation samples and extreme thresholds, VoI
calculations might generate counter-intuitive results11. One instance is when
sample estimates of $\theta_{se}$ or $\theta_{sp}$ are at their extreme (0 or
1). The improper $\text{Beta}(0,0)$ prior on sensitivity and specificity can
result in improper beta posteriors with one of their two parameters being 0.
Therefore, we suggest specifying proper priors (e.g., a flat prior
$\text{Beta}(1,1)$) for the three quantities in $\theta$ when current
information is from an observed sample.

### A general algorithm for arbitrary specification of $P(\theta|\mathbf{d})$

In many practical situations our current information $P(\theta|\mathbf{d})$
may not have a simple mathematical form, but one can still draw a Monte Carlo
sample of arbitrary size from its distribution (e.g., via Markov Chain Monte
Carlo methods). Examples include when $P(\theta_{p})$ is constructed from a
meta-analysis of prevalence studies, or $P(\theta_{se},\theta_{sp})$ from a
bivariate meta-analysis of the diagnostic accuracy of the test, giving rise to
a joint distribution for $P(\theta_{se},\theta_{sp})$, potentially conditional
on prevalence25.

For such a general case where we have a sample from $P(\theta|\mathbf{d})$,
the empirical distribution of the sample can be used as the distribution for
the current information. The posterior distribution after observing future
data can then be approximated by the discrete distribution that puts a mass
$w^{(i)}\propto P(\mathcal{D}^{*}|\theta^{(i)})$ (e.g., a product of Binomial
probabilities for binary outcomes) on each observation from this sample. This
will result in the generic algorithm explained in Table 3.

0. Quantify current information given arbitrary distribution $P(\theta|\mathbf{d})$. Obtain $M$ draws $(\theta^{(1)},\theta^{(2)},...,\theta^{(M)})$ from $P(\theta|\mathbf{d})$. 1. Estimate expected NBs under current information by taking expectation of NB terms (in equation 1) with respect to the empirical distribution that puts a mass $w^{(k)}$ of $\propto 1$ on each $\theta^{(k)}$:
$NB(i,\mathbf{\theta})=\left\\{\begin{array}[]{ll}0&i=0\\\
\sum_{k=1}^{M}w^{(k)}\\{\theta^{(k)}_{p}\theta^{(k)}_{se}-(1-\theta^{(k)}_{p})(1-\theta^{(k)}_{sp})\frac{z}{1-z}\\}/\sum_{k=1}^{M}w^{(k)}&i=1\\\
\sum_{k=1}^{M}w^{(k)}\\{\theta^{(k)}_{p}-(1-\theta^{(k)}_{p})\frac{z}{1-z}\\}/\sum_{k=1}^{M}w^{(k)}&i=2\\\
\end{array}\right..$ Record the maximum expected $NB$. This is
$\mathbb{E}NB_{\text{current information}}$. 2. For $j$ = 1 to $N$ (number of
Monte Carlo simulations) (a) Obtain a random draw from $P(\theta|\mathbf{d})$
(or loop over the draws from step  ‣ 3). This is taken as the true $\theta$ in
this iteration. (b) Calculate true NBs given $\theta$ obtained in the previous
step. Record the maximum NB. (c) Generate $\mathcal{D^{*}}$ of size $n^{*}$,
defined by $\\{n^{*}_{tp},n^{*}_{fn},n^{*}_{tn},n^{*}_{fp}\\}$ given the
sample from $\theta$ obtained in step 2a:
$n^{*}_{+}\sim\text{Binomial}(n^{*},\theta_{p})$ (number of positive cases in
the future sample),
$n^{*}_{tp}\sim\text{Binomial}(n^{*}_{+},\theta_{se})$,
$n^{*}_{fn}=n^{*}_{+}-n^{*}_{tp}$,
$n^{*}_{tn}\sim\text{Binomial}(n^{*}-n^{*}_{+},\theta_{sp})$,
$n^{*}_{fp}=n^{*}-n^{*}_{+}-n^{*}_{tn}$. (d) Update $P(\theta|\mathbf{d})$ by
updating the weights $w^{(k)}$ in the empirical distribution of
$\theta|\mathbf{d}$ (from step 1) based on the likelihood
$P(\mathcal{D^{*}}|\theta)$. For binary outcomes the weights are:
$w^{(k)}\propto P_{Bin}(n^{*}_{+},n^{*},\theta_{p}^{(k)})\times
P_{Bin}(n^{*}_{tp},n^{*}_{+},\theta_{se}^{(k)})\times
P_{Bin}(n^{*}_{tn},n^{*}-n^{*}_{+},\theta_{sp}^{(k)})$,
with $P_{Bin}(x,m,p)$ being the probability mass function of the binomial
distribution at value $x$ for $m$ trials and success probability $p$. (e)
Calculate updated NBs given future data by taking expectation of NB terms in
step 1 using the updated weights from step 2d. Next $j$ 3. Average NBs from
from step 2b. This is $\mathbb{E}NB_{\text{perfect information}}$. Subtract
from this ENB from step 1. This is EVPI. 4. Average NBs from step 2e. This is
$\mathbb{E}NB_{\text{sample information}}$. Subtract from this ENB from step
1. This is EVSI.

Table 3: EVSI computation for binary outcomes for general, non-conjugate
distributions

## Results

We used data from GUSTO-I, a clinical trial of multiple thrombolytic
strategies for acute myocardial infarction (AMI), as our case study26. This
dataset has been widely used for methodological research in predictive
analytics27. All analyses are conducted in R28.

For the present case study, we use a similar set-up as our previous work, with
a risk prediction model for 30-day mortality developed using the non-US sample
of this dataset, and validated in the US sub-sample. The risk prediction model
is a logit model fitted to the entire non-US subset (sample size 17796), with
the resulting equation:

logit(P(Y=1))= - 2.0842 + 0.0781[age] + 0.4027[AMI location other (vs.
inferior)] + 0.5773[AMI location anterior (v. inferior)] + 0.4678[previous
AMI] + 0.7666[AMI severity (Killip score)] - 0.0775[min(blood pressure, 100)]
+ 0.0182[pulse].

Similar to our previous work, we initially assume we have access to only
$n=500$ individuals from the US sub-sample as the source of current
information. We randomly selected 500 individuals, without replacement, from
the entire 23034 observations in the US sub-sample of GUSTO-I. Later, we will
use more observations from the US sub-sample to study how EVSI behaves with
the amount of current information. The prevalence of the outcome in the
development sample (entire non-US sub-sample) was 0.0723; in the current
validation sample it was 0.0860, while in the entire US sub-sample it was
0.0679. We consider the range of relevant thresholds to be 1%-10%.

We start by drawing the decision curve at this range of thresholds, as shown
in Figure 1. Panel (a) shows the decision curve, and panel (b) shows the
difference between the NB of the model and the best alternative default
strategy ($NB_{1}-\max(NB_{0},NB_{2})$) along with its 95% confidence interval
(CI, based on the percentile methods using 10,000 bootstraps). The model had
higher expected NB than the alternative strategies at the entire range of
thresholds of interest.

---
(a) Decision curve

---
(b) Incremental NB of the model v. best alternative strategy

Figure 1: Decision curve (a) and incremental net benefit of the model v. the
best alternative strategy (b). Dashed lines are 95% confidence intervals

Next, we draw the EVPI curve. EVPI and EVSI computations are based on 10^6
Monte Carlo simulations using the aggregate beta-binomial model, with a flat
$\text{Beta}(1,1)$ prior on each component of $\theta$. The EVPI at 0.01
threshold was 0.00039. At the 0.02 threshold it was 0.00125. The EVPI was 0 at
0.05 and 0.10 thresholds (so EVSI for any $n^{*}$ will also be 0). We
therefore focus on VoI calculations at 0.01 and 0.02 thresholds.

The EVPI curve is presented in Figure 2. EVPI quantifies the expected NB loss
per decision (every time the model might be used). The overall NB loss due to
uncertainty is therefore affected by the expected number of times the medical
decision will be made29. In scaling VoI quantities to the population, we apply
a similar approach as in our previous work11: a decision-making body in charge
of making a national recommendation can consider the model to be applicable to
all AMI cases in the US (the target population of the validation sample).
There are around 800,000 AMIs in the US every year30, and a guideline
development panel in charge of recommending risk stratification for AMI
management can consider all those instances relevant. The scaled VoI values
are provided in true positive units on the second Y-axis of the EVPI figure.
At the 0.02 threshold, a future validation study can have a maximum benefit
equal to gaining 1001 net true positive (individuals who will die within 30
days and are correctly identified as high risk) per year, or avoiding 49069
extra false positive cases per year (individuals who will not die within 30
days.

Figure 2: The expected value of perfect information (EVPI) of the validation
sample

The EVSI curves are shown in Figure 3. Generally, it is expected that EVSI
will increase with a larger future study, and will asymptote to EVPI, a
pattern that is obvious for both thresholds. The ‘diminishing return’ pattern
is also obvious for both thresholds: the expected NB gain steeply rises at
small samples and plateaus as $n^{*}$ grows. Similar to EVPI, population
scaling is applied to EVSIs and is presented as the second Y-axis of the EVSI
curve. At 0.02 threshold, a future study of size $n^{*}=500$ has a per-
decision EVSI value of 0.00100, corresponding to population value of 804 in
true positive cases gained, or 39381 in false positive cases averted.

---
(a) $z=0.01$

---
(b) $z=0.02$

Figure 3: EVSI for the case study

Figure 4 provides the results of brief simulation studies that demonstrate how
EVSI changes as a function of the amount of current information, represented
by the sample size of $\mathbf{d}$, the current validation data from which
$P(\theta|\mathbf{d})$ is constructed. Starting from a value of $500$, we
doubled the sample size of the current validation sample, specified the
distributions for $\theta_{p}$, $\theta_{se}$, and $\theta_{sp}$ based on this
sample, and performed EVSI calculations as above, until reaching $n=8,000$.
Results are provided for both 0.01 and 0.02 thresholds and are the average of
100 independent simulations, in each a new sample (with replacement) from the
US subset of GUSTO-I was obtained as the current validation data (100
simulations were considered adequate based on inspecting the standard error of
the Monte Carlo simulations).

As expected, the higher the amount of current information (the larger the
$n$), the lower the gain in NB by conducting future studies. When the current
validation sample is small, the initial gain in NB by conducting future
research can be very large, as evidenced by the steep slope of the EVSI curve.
This is not the case when the current information is relatively robust (e.g.,
based on $n$=8,000, which includes on average 544 events). The comparison
between the two thresholds is also informative, as it demonstrates that the
amount of gain from future sample is highly dependent on the threshold, even
in relative terms. For example, at $z=0.01$, the EVSI curve when current
information is based on $n$=1,000 (dotted blue line) is only slightly lower
than the EVSI curve when $n$=500 (dotted black line). At $z=0.02$, EVSI when
current information is based on $n$=1,000 (solid blue line) is significantly
lower than when the current information is based on $n$=500 (solid black
line), across the entire range of future sample sizes.

Figure 4: The expected value of sample information (EVSI) as a function of the
sample sample size of the future study ($\mathcal{D^{*}}$) for various levels
of current information, represented by the sample size of the current study
($n$) from which $P(\theta)$ is constructed

## Discussion

Evaluating the performance of a model in a finite sample is fraught with
uncertainty. In this work, we applied the VoI framework to investigate the
decision-theoretic consequences of such uncertainty. We defined validation
EVSI as the expected gain in NB by obtaining an external validation sample of
a given size from the target population of interest. In a case study we showed
the feasibility of EVSI calculations, and studied how EVSI is affected by the
amount of current information and the sample size of the future study. We
suggested scaling the EVSI to the population to quantify the overall gain in
clinical utility in true (or false) positive units from conducting an external
validation study of a given sample size.

How can EVSI analysis inform study design in predictive analytics? Conducting
a validation study is an investment in resources that will generate further
information on the performance of a clinical prediction model in a target
population. The decision whether to undertake a validation study should
ultimately hinge on whether the information gained from such a study is worth
the required investments. The EVSI, when scaled to the population, determines
the expected return on investment in NB unit. Ultimately, such return should
be contrasted against the efforts and resources required for such a study. In
decision-analytic (health policy) modeling, EVSI is typically in net monetary
units, and when scaled to population, can be compared with the budget of a
planned data collection activity. The optimal sample size will be one that
maximizes the difference between population EVSI and study costs29. The NB for
risk prediction models, on the other hand, is in net true or false positive
units, and as such cannot directly be compared with the budget of a validation
study. One can always embark on full decision analysis to translate all
outcomes to net monetary units, but this will likely require sophisticated
decision modeling and context-specific assumptions on long-term outcomes, a
process that might take significant amount of time and require further data
collection (e.g., to obtain utility weights for outcomes). To us, a main
reason for the vast popularity of decision curve analysis is that it provides
an assumption-free, reproducible method for NB calculation based on the very
same data that are used for studying model calibration and discrimination. We
think VoI analysis during model development and validation should generally
keep the same spirit. A full decision analysis should be relegated to after an
impact analysis has measured the resource-use implications of implementing the
model. This, however, means the decision rules for determining the optimal
sample size of development and validation studies based on the VoI framework
would be different than those used in health policy analysis. This paper
deliberately stayed away from proposing such rules, and instead focused on
defining concepts and proposing computation methods for EVSI. Proposing such
decision rules is detached from EVSI calculations and deserves its own airing.

There are multiple areas of further inquiry. The EVSI framework should also be
applied to the development phase of prediction models. This can guide the
investigator on whether further development, or moving to validation, should
be prioritized. We mainly focused on NB loss due to sampling uncertainty.
However, there are several sources of uncertainty, such as whether our
existing information on the performance of the model is directly applicable to
the target population, or if predictors and outcome are measured with the same
quality between the study and usual practice. The comparative statistical and
computational performance of the EVSI computations algorithms, and the
adequacy of a given number of simulations for each algorithm, should be
evaluated in dedicated studies. We also do not claim the algorithms we
proposed are the only ones that can be used for EVSI computations. Other
algorithms, such as those based on central limit theorem31,32, can prove
useful. VoI analysis in decision modeling has received significant boost in
computational speed in recent years due to the arrival of algorithms based on
non-parametric regression modeling33,34. This approach can facilitate VoI
analysis in risk prediction as well. The EVSI defined in this work is for a
single, homogeneous target population, and does not consider heterogeneous
settings. Consider the hierarchical Bayesian tri-variate model by Wynants et
al for NB estimation in multi-center validation studies19 and imagine a new
center is interested in validating this model in their local setting. The
posterior predictive distribution from the hierarchical model for a new center
(potentially combined with any local information on the performance of the
model) can act as ‘current information’, enabling EVSI calculations. This will
be valid and informative insofar as the expected NB gain in that center is
concerned. However, in a hierarchical Bayesian framework, information learned
from one cluster propagates to others, and thus conducting a validation study
in a new center results in change in the distribution of NB in other centers.
Such NB gain from learning across clusters in hierarchical settings is not
reflected in the EVSI computed for this particular setting. VoI metrics and
corresponding computation algorithms for multi-center studies should be
developed separately. Further, during external validation, often a secondary
aim is to update the model if its performance turns out to be sub-optimal35.
Such model revision can take different levels of complexity (from intercept
correction to re-estimating regression coefficients36) and one might be
interested in the expected yield of a given sample size for such model
revision. This will get connected to the VoI concepts for model development.
As stated earlier, how VoI metrics for prediction models should inform
objective functions for determining the optimal sample size should be debated
by the community in the hope of generating consensus and best practice
standards.

There is an ongoing debate on the appropriateness of conventional, frequentist
metrics of uncertainty when reporting NB of a clinical prediction model3–5.
The same concerns can logically be extended to frequentist method for sample
size and power calculations around NB37–39. VoI methodology provides a
rigorous, utilitarian response to such controversies. The toolbox of VoI
methods for clinical prediction models is growing, and perhaps it is time to
formalize the role of VoI in uncertainty quantification and design of
empirical studies in predictive analytics.

## References

re 1. Steyerberg EW, Vergouwe Y. Towards better clinical prediction models:
Seven steps for development and an ABCD for validation. _European Heart
Journal_ 2014; 35: 1925–1931.

pre 2. Fitzgerald M, Saville BR, Lewis RJ. Decision curve analysis. _JAMA_
2015; 313: 409–410.

pre 3. Kerr KF, Marsh TL, Janes H. The Importance of Uncertainty and Opt-In v.
Opt-Out: Best Practices for Decision Curve Analysis. 2019; 39: 491–492.

pre 4. Capogrosso P, Vickers AJ. A Systematic Review of the Literature
Demonstrates Some Errors in the Use of Decision Curve Analysis but Generally
Correct Interpretation of Findings. _Med Decis Making_ 2019; 39: 493–498.

pre 5. Vickers AJ, Van Claster B, Wynants L, et al. Decision curve analysis:
Confidence intervals and hypothesis testing for net benefit. _Diagn Progn Res_
2023; 7: 11.

pre 6. Felli J, Hazen G. Sensitivity analysis and the expected value of
perfect information. 1998; 18: 95–109.

pre 7. Jackson CH, Baio G, Heath A, et al. Value of Information Analysis in
Models to Inform Health Policy. _Annu Rev Stat Appl_ 2022; 9: 95–118.

pre 8. Heath A, Kunst N, Jackson C. _Value of Information for Healthcare
Decision-Making_. 1st ed. Boca Raton: Chapman; Hall/CRC. Epub ahead of print
December 2023. DOI: 10.1201/9781003156109.

pre 9. Tuffaha HW, Gordon LG, Scuffham PA. Value of information analysis in
healthcare: A review of principles and applications. _J Med Econ_ 2014; 17:
377–383.

pre 10. Sadatsafavi M, Yoon Lee T, Gustafson P. Uncertainty and the value of
information in risk prediction modeling. _Medical Decision Making: An
International Journal of the Society for Medical Decision Making_ 2022;
272989X221078789.

pre 11. Sadatsafavi M, Lee TY, Wynants L, et al. Value-of-Information Analysis
for External Validation of Risk Prediction Models. _Medical Decision Making:
An International Journal of the Society for Medical Decision Making_ 2023; 43:
564–575.

pre 12. Ades A, Lu G, Claxton K. Expected value of sample information
calculations in medical decision modeling. 2004; 24: 207–227.

pre 13. Vickers AJ, Cronin AM, Elkin EB, et al. Extensions to decision curve
analysis, a novel method for evaluating diagnostic tests, prediction models
and molecular markers. _BMC Med Inform Decis Mak_ 2008; 8: 53.

pre 14. Vickers AJ, Kattan MW, Sargent DJ. Method for evaluating prediction
models that apply the results of randomized trials to individual patients.
_Trials_ ; 8. Epub ahead of print December 2007\. DOI: 10.1186/1745-6215-8-14.

pre 15. Glynn D, Nikolaidis G, Jankovic D, et al. Constructing Relative Effect
Priors for Research Prioritization and Trial Design: A Meta-epidemiological
Analysis. _Med Decis Making_ 2023; 43: 553–563.

pre 16. Ibrahim JG, Chen M-H, Gwon Y, et al. The power prior: Theory and
applications. _Stat Med_ 2015; 34: 3724–3749.

pre 17. O’Hagan A. _Uncertain judgements: Eliciting experts’ probabilities_.
London: Wiley, 2006.

pre 18. Vickers AJ, Elkin EB. Decision curve analysis: A novel method for
evaluating prediction models. 2006; 26: 565–574.

pre 19. Wynants L, Riley RD, Timmerman D, et al. Random-effects meta-analysis
of the clinical utility of tests and prediction models. _Stat Med_ 2018; 37:
2034–2052.

pre 20. Cruz GNF, Korthauer K. Bayesian Decision Curve Analysis with bayesDCA,
http://arxiv.org/abs/2308.02067 (2023, accessed 4 January 2024).

pre 21. Rubin DB. The Bayesian Bootstrap. _Ann Statist_ ; 9. Epub ahead of
print January 1981. DOI: 10.1214/aos/1176345338.

pre 22. Sadatsafavi M, Marra C, Bryan S. Two-level resampling as a novel
method for the calculation of the expected value of sample information in
economic trials. _Health Econ_ 2013; 22: 877–882.

pre 23. Schafer J. Multiple imputation: A primer. _Statistical Methods in
Medical Research_ 1999; 8: 3–15.

pre 24. Broderick T, Jordan MI, Pitman J. Beta Processes, Stick-Breaking and
Power Laws. _Bayesian Anal_ ; 7. Epub ahead of print June 2012. DOI:
10.1214/12-BA715.

pre 25. Li J, Fine JP. Assessing the dependence of sensitivity and specificity
on prevalence in meta-analysis. _Biostatistics_ 2011; 12: 710–722.

pre 26. investigators G. An international randomized trial comparing four
thrombolytic strategies for acute myocardial infarction. _The New England
Journal of Medicine_ 1993; 329: 673–682.

pre 27. Ennis M, Hinton G, Naylor D, et al. A comparison of statistical
learning methods on the gusto database. _Statistics in Medicine_ 1998; 17:
2501–2508.

pre 28. R Core Team. _R: A Language and Environment for Statistical
Computing_. Vienna, Austria: R Foundation for Statistical Computing,
https://www.R-project.org/ (2019).

pre 29. Willan AR, Goeree R, Boutis K. Value of information methods for
planning and analyzing clinical studies optimize decision making and research
planning. _J Clin Epidemiol_ 2012; 65: 870–876.

pre 30. Virani SS, Alonso A, Aparicio HJ, et al. Heart Disease and Stroke
Statistics-2021 Update: A Report From the American Heart Association.
_Circulation_ 2021; 143: e254–e743.

pre 31. Marsh TL, Janes H, Pepe MS. Statistical inference for net benefit
measures in biomarker validation studies. _Biometrics_ 2020; 76: 843–852.

pre 32. Yoon Lee T, Gustafson P, Sadatsafavi M, et al. Closed-Form Solution of
the Unit Normal Loss Integral in 2 Dimensions, with Application in Value-of-
Information Analysis. _Med Decis Making_ 2023; 43: 621–626.

pre 33. Heath A, Manolopoulou I, Baio G. A Review of Methods for Analysis of
the Expected Value of Information. _Med Decis Making_ 2017; 37: 747–758.

pre 34. Heath A, Kunst N, Jackson C, et al. Calculating the Expected Value of
Sample Information in Practice: Considerations from 3 Case Studies. _Med Decis
Making_ 2020; 40: 314–326.

pre 35. Steyerberg EW, Borsboom GJJM, Houwelingen HC van, et al. Validation
and updating of predictive logistic regression models: A study on sample size
and shrinkage. _Statistics in Medicine_ 2004; 23: 2567–2586.

pre 36. Vergouwe Y, Nieboer D, Oostenbrink R, et al. A closed testing
procedure to select an appropriate method for updating prediction models.
_Stat Med_ 2017; 36: 4529–4539.

pre 37. Riley RD, Debray TPA, Collins GS, et al. Minimum sample size for
external validation of a clinical prediction model with a binary outcome.
_Stat Med_ 2021; 40: 4230–4251.

pre 38. Riley RD, Snell KIE, Ensor J, et al. Minimum sample size for
developing a multivariable prediction model: Part Ib -b Continuous outcomes.
_Stat Med_ 2019; 38: 1262–1275.

pre 39. Riley RD, Snell KI, Ensor J, et al. Minimum sample size for developing
a multivariable prediction model: PART II - binary and time-to-event outcomes.
_Stat Med_ 2019; 38: 1276–1296.

p