OVERVIEW

This Chapter Explains…

• Likelihood methods for fitting log-location-scale distributions (especially the Weibull and lognormal distributions)

• Likelihood confidence intervals/regions for model parameters and for functions of model parameters

• Normal-approximation confidence intervals/regions

• Estimation and confidence intervals for log-location-scale distributions with a given shape parameter

8.1 - INTRODUCTION

This chapter extends the methods presented in Chapter 7

• Distributions with two parameters

• Based on location-scale distributions

• These methods can be applied to any distribution in the location-scale family

• Particular attention is paid to the Weibull and lognormal distributions

• For exact failure times the density approximation to the likelihood (Eq 7.13) is adequate

8.2.1 - \(\mathcal{L}\) for Location-Scale Distributions

• For data containing only exact failures and right censored observations, the likelihood function for a location-scale distribution may be expressed as

\[ \begin{aligned} \mathcal{L}(\mu,\sigma|\underline{t})&=\prod_{i=1}^n\left[f(t_i|\mu,\sigma)\right]^{\delta_{i}}\left[1-F(t_i|\mu,\sigma)\right]^{1-\delta_i}\\\\ &=\prod_{i=1}^n\left[\frac{1}{\sigma}\phi\left(\frac{t_i-\mu}{\sigma}\right)\right]^{\delta_i}\times\left[1-\Phi\left(\frac{t_i-\mu}{\sigma}\right)\right]^{1-\delta_i} \end{aligned} \]

• Where

\[ \delta_i=\begin{cases}1 &\mbox{for an exact observation}\\\\ 0 &\mbox{for a right-censored observation}\end{cases} \]

8.2.2 - \(\mathcal{L}\) for Log-Location-Scale Distributions

• For data containing only exact failures and right censored observations, the likelihood function for a location-scale distribution may be expressed as

\[ \begin{aligned} \mathcal{L}(\mu,\sigma|\underline{t})&=\prod_{i=1}^n\left[f(\log(t_i)|\mu,\sigma)\right]^{\delta_{i}}\left[1-F(\log(t_i)|\mu,\sigma)\right]^{1-\delta_i}\\\\ &=\prod_{i=1}^n\left[\frac{1}{t_i\sigma}\phi\left(\frac{\log(t_i)-\mu}{\sigma}\right)\right]^{\delta_i}\times\left[1-\Phi\left(\frac{\log(t_i)-\mu}{\sigma}\right)\right]^{1-\delta_i} \end{aligned} \]

• Where

\[ \delta_i=\begin{cases}1 &\mbox{for exact observations}\\\\ 0 &\mbox{for right-censored observations}\end{cases} \]

Example 8.1 Shock Absorber

SMRD::smrd_app('shock_absorber8')

8.3.1 - Joint Confidence Regions for \(\mu\) and \(\sigma\)

Confidence Intervals & Confidence Regions

• For one unknown parameter (i.e. \(\overline{\theta}=\theta\)) an approximate \(100(1-\alpha)\%\) confidence interval includes all values of \(\theta\) such that

\[ R(\theta) \ge \exp\left[-\chi^2_{(1-\alpha;1)}/2\right] \]

• A plot of a likelihood-based confidence interval for a single parameter depicts \(R(\theta)\) for a range of \(\theta\) values

• The shiny app below shows a likelihood-based confidence interval plot for the exponential parameter \(\theta\) using the berkson200 data set

• For two unknown parameters (i.e. \(\overline{\theta}=\theta_1,\theta_2\)) an approximate \(100(1-\alpha)\%\) confidence region is all values of \(\theta_1\) and \(\theta_2\) such that

\[ R(\theta_1, \theta_2) \ge \exp\left[-\chi^2_{(1-\alpha;2)}/2\right] \]

Figures 8.3 & 8.4 - Joint Confidence & Relative Likelihood Regions

par(mfrow = c(1,2), las = 1)

library(SMRD)

zoom <- 1.5
threeD <- F

Shock.ld <- 
frame.to.ld(shockabsorber,
            response.column = 1, 
            censor.column = 3,
            time.units = 'Kilometers')

simple.contour(Shock.ld, 
               distribution = 'weibull', 
               show.confidence = T, 
               zoom.level = zoom,
               original.par = F,
               threeD = threeD)

simple.contour(Shock.ld, 
               distribution = 'weibull', 
               show.confidence = F, 
               zoom.level = zoom,
               threeD = threeD)

8.3.2 - Individual CI for \(\mu\) and \(\sigma\)

Procedure for computing likelihood-based CI’s for \(\mu\)

1. Estimate \(\hat{\mu}_{_{MLE}}\;\text{and}\;\hat{\sigma}_{_{MLE}}\)

2. Select a value \(\mu_0\)

3. For the selected value of \(\mu_0\), numerically find the value of \(\sigma\) that maximizes the relative likelihood, i.e.

\[ R(\mu_0)=\max_{\sigma}\left[\frac{\mathcal{L}(\mu_0,\sigma)}{\mathcal{L}(\hat{\mu}_{_{MLE}},\hat{\sigma}_{_{MLE}})}\right] \]

4. Repeat steps (2) and (3) for many values of \(\mu_0\). Record and plot \(\left(\mu_0, R(\mu_0)\right)\) - these values form what is known as the profile likelihood for \(\mu\)

Procedure for computing likelihood-based CI’s for \(\sigma\)

1. Estimate \(\hat{\mu}_{_{MLE}}\;\text{and}\;\hat{\sigma}_{_{MLE}}\)

2. Select a value \(\sigma_0\)

3. For the selected value of \(\sigma_0\), numerically find the value of \(\mu\) that maximizes the relative likelihood, i.e.

\[ R(\sigma_0)=\max_{\mu}\left[\frac{\mathcal{L}(\mu,\sigma_0)}{\mathcal{L}(\hat{\mu}_{_{MLE}},\hat{\sigma}_{_{MLE}})}\right] \]

4. Repeat steps (2) and (3) for many values of \(\sigma_0\). Record and plot \((\sigma_0, R(\sigma_0))\) - these values form what is known as the profile likelihood for \(\sigma\)

Figures 8.5 & 8.6 - Profile Likelihoods \(R(\mu)\) & \(R(\sigma)\)

• The code chunk below shows how to generate plots of the profile likelihood for both \(\mu\) and \(\sigma\)

par(mfrow = c(1,2), las = 1)

library(SMRD)

Shock.ld <- 
frame.to.ld(shockabsorber,
            response.column = 1, 
            censor.column = 3,
            time.units = 'Kilometers')

simple.contour(Shock.ld, 
               distribution = 'weibull', 
               profile = 'x', 
               size = 300,
               print.ci = F)

simple.contour(Shock.ld, 
               distribution = 'weibull', 
               profile = 'y', 
               size = 300,
               print.ci = F)

8.3.3 - Likelihood CI for Functions of \(\mu,\sigma\)

Once \(\hat{\mu}_{_{MLE}}\;\text{and}\;\hat{\sigma}_{_{MLE}}\) have been estimated

• Confidence intervals may be computed for functions of \(\mu\) and/or \(\sigma\)

  • Define a \(1-1\) transformation \(\mathbf{g[\mu,\sigma]}=[g_1(\mu,\sigma),g_2(\mu,\sigma)],\) where \(g_1\;\text{or}\;g_2\) is the function of interest

  • Due to the invariance properties of MLE’s, likelihood-based CI methods can be applied as was described for the parameters

Likelihood Confidence Interval for \(t_p\)

• For the quantile function \(t_p\)

\[ R(t_p)=\max_{\sigma}\left[\frac{\mathcal{L}(f(\mu),\sigma)}{\mathcal{L}(\hat{\mu}_{_{MLE}},\hat{\sigma}_{_{MLE}})}\right] \]

• Recall, \(t_p=\exp\left[\mu+\Phi^{-1}(p)\sigma\right]\), rearranging gives

\[ \mu = \log(t_p)-\Phi^{-1}(p)\sigma \]

• Substituting the right-side of the above expression for \(\mu\) into the relative likelihood gives

\[ R(t_p)=\max_{\sigma}\left[\frac{\mathcal{L}\left[\log(t_p)-\Phi^{-1}(p)\sigma,\sigma\right]}{\mathcal{L}(\hat{\mu}_{_{MLE}},\hat{\sigma}_{_{MLE}})}\right] \]

• Therefore, to find the likelihood-based confidence interval for \(t_p\)

  1. Choose a value for \(\mu_0\) and find \(R(\mu_0)\) for each value of \(\sigma\)

  2. Choose a quantile value \(p\), we know that \(R(\mu_0) = R(t_{p_{0}})\)

  3. Record and plot \((t_{p_{0}},R(t_{p_{0}}))\) to get the likelihood profile for \(t_p\)

Figure 8.7 & 8.8 - Likelihood Contour & Profile for \(t_{0.1}\)

par(mfrow = c(1,2), las = 1)

library(SMRD)

Shock.ld <- 
frame.to.ld(shockabsorber,
            response.column = 1, 
            censor.column = 3,
            time.units = 'Kilometers')

simple.contour(Shock.ld, 
               distribution = 'weibull', 
               quantile = 0.1, 
               size = 300, 
               zoom.level = 2.25, 
               show.confidence = F,
               original.par = F)

simple.contour(Shock.ld, 
               distribution = 'weibull', 
               quantile = 0.1, 
               profile = 'x', 
               size = 300)

Likelihood Confidence Interval For \(F(t_e)\)

• Similarly, \(F(t_e)=\Phi[\left(\log(t_e)-\mu\right)/\sigma]\), rearranging gives

\[ \mu = \log(t_e)-\Phi^{-1}\left(F(t_e)\right)\sigma \]

• substituting the right-side of this expression for \(\mu\) into the relative likelihood

\[ R\left(F(t_e)\right)=\max_{\sigma}\left[\frac{\mathcal{L}\left[\log(t_e)-\Phi^{-1}\left(F(t_e)\right)\sigma,\sigma\right]}{\mathcal{L}(\hat{\mu}_{_{MLE}},\hat{\sigma}_{_{MLE}})}\right] \]

• To find the likelihood-based CI for \(F(t_e)\)

  • Choose a value for \(\mu_0\) and find \(R(\mu_0)\) for each value of \(\sigma\)

  • For the time of interest \(t_e\), we know that \(R(\mu_0) = R\left(F(t_{e_{0}})\right)\)

  • Record and plot \(\left(F(t_{e_{0}}),R\left(F(t_{e_{0}})\right)\right)\) to get the likelihood profile for \(F(t_e)\)

par(family = 'serif', mfrow = c(1,2), las = 1, cex = 1.25)

library(SMRD)

ShockAbsorber.ld <- 
frame.to.ld(SMRD::shockabsorber,
            response.column = 1, 
            censor.column = 3,
            time.units = 'Kilometers')

simple.contour(ShockAbsorber.ld, 
               distribution = 'weibull', 
               quantile = 0.1, 
               size = 300, 
               zoom = 2.25, 
               rel.or.conf = 'Relative likelihood',
               original.par = F)

simple.contour(ShockAbsorber.ld, 
               distribution = 'weibull', 
               quantile = 0.1, 
               profile = 'x', 
               size = 300)

8.4.1 - Parameter Variance-Covariance Matrix

• For location-scale distributions, a local estimate of \(\Sigma_{\hat{\theta}}\) may be computed by inverting the observed information matrix

\[ \mathbf{\widehat{\Sigma}_{\hat{\mu},\hat{\sigma}}} = \left[ \begin{array}{cc} \widehat{Var}(\hat{\mu}) & \widehat{Cov}(\hat{\mu},\hat{\sigma}) \\ \widehat{Cov}(\hat{\mu},\hat{\sigma}) & \widehat{Var}(\hat{\sigma}) \\ \end{array} \right]=\left[ \begin{array}{cc} -\frac{\partial^2\mathcal{L}}{\partial\mu^2} & -\frac{\partial^2\mathcal{L}}{\partial\mu\partial\sigma} \\ -\frac{\partial^2\mathcal{L}}{\partial\sigma\partial\mu} & -\frac{\partial^2\mathcal{L}}{\partial\sigma^2} \\ \end{array} \right]^{-1} \]

• As was discussed in Chapter 7, \(\mathbf{\widehat{\Sigma}_{\hat{\mu},\hat{\sigma}}}\) describes the curvature of \(\mathcal{L}\) at \(\hat{\mu},\hat{\sigma}\)

• A higher degree of curvature implies a more concentrated likelihood near \(\hat{\mu},\hat{\sigma}\) and thereby greater precision.

Once Again, The SMRD Package Makes This Easy

• We can create a list of several important maximum likelihood estimate properties using the mlest function

• The mlest function requires that we provide

  • A life data object

  • A distribution for which we want the ML properties

shock.ld <- frame.to.ld(shockabsorber,
                        response.column = 1,
                        censor.column = 3)

shock.mlest <- mlest(shock.ld, 
                     distribution = "weibull")

• We can then return \(\widehat{\Sigma}_{\hat{\mu},\hat{\sigma}}\) for the shockabsorber data set (assuming a Weibull distribution) with the following code

shock.vcv <- print(shock.mlest)$vcv.matrix ; shock.vcv

             mu     sigma
mu    0.0120759 0.0039904
sigma 0.0039904 0.0053532

• Finally, we can convert this object into a formatted table using

knitr::kable(shock.vcv, digits = c(6,6))

• Similarly, we can return the parameter correlation \(\widehat{\rho}_{\hat{\mu},\hat{\sigma}}=\widehat{Cov}(\hat{\mu},\hat{\sigma})/\sqrt{\widehat{Var}(\hat{\mu})\widehat{Var}(\hat{\sigma})}\)

knitr::kable(print(shock.mlest)$param.corr.matrix,
             format = 'markdown',
             digits = c(1,6,6))

8.4.2 - Confidence Intervals for Model Parameters

For location-scale and log-location-scale distributions

• Confidence intervals based on the normal approximation assume that

\[ Z_{\hat{\mu}}=\frac{\hat{\mu}-\mu}{\widehat{se}_{\hat{\mu}}}\sim NOR(0,1) \]

• and

\[ Z_{\log(\hat{\sigma})}=\frac{\log(\hat{\sigma})-\log(\sigma)}{\widehat{se}_{\log(\hat{\sigma})}}\sim NOR(0,1) \]

• Where the \(\log(\hat{\sigma})\) transformation is made to improve the coverage probability of the CI since \(\sigma\) is a strictly positive quantity

• Therefore a \(100(1-\alpha)\%\) 2-sided CI for \(\mu\) is expressed as

\[ \left[\underline{\mu},\overline{\mu}\right]=\hat{\mu}\pm z_{(1-\alpha/2)}\widehat{se}_{\hat{\mu}} \]

• A two-sided \(100(1-\alpha)\%\) CI for \(\sigma\) is expressed as

\[ \left[\underline{\sigma},\overline{\sigma}\right]=[\hat{\sigma}/w, \hat{\sigma}\times w] \]

• where

\[ w = \exp[z_{(1-\alpha/2)}\widehat{se}_{\hat{\sigma}}/\hat{\sigma}] \]

8.4.3 - Normal-Approximation CI for Functions of \(\mu,\sigma\)

For location-scale and log-location-scale distributions

• Computing confidence intervals for \(g_1=f(\mu,\sigma)\) based on the normal approximation assumes that

\[ Z_{\hat{g}_1}=\frac{\hat{g}_1-g_1}{\widehat{se}_{\hat{g}_1}}\sim NOR(0,1) \]

• Where approximate \(100(1-\alpha)\%\) confidence interval for \(g_1\) is then expressed as

\[ \left[\underline{g_1}, \overline{g_1}\right]=\hat{g_1}\pm z_{(1-\alpha/2)}\widehat{se}_{\hat{g_1}} \]

• And

\[ \begin{aligned} \widehat{se}_{\hat{g}_1}&=\sqrt{\widehat{Var}(\hat{g}_1)}\\ &=\left[\left(\frac{\partial g_1}{\partial\mu}\right)^2\widehat{Var}(\hat{\mu})+2\left(\frac{\partial g_1}{\partial\mu}\right)\left(\frac{\partial g_1}{\partial\sigma}\right)\widehat{Cov}(\hat{\mu},\hat{\sigma})+\left(\frac{\partial g_1}{\partial\sigma}\right)^2\widehat{Var}(\hat{\sigma})\right]^{1/2} \end{aligned} \]

• Note any monotone function of \(\mu\) transmits directly to \(\left[\underline{\mu},\overline{\mu}\right]\) for example

\[ \left[\underline{\eta},\overline{\eta}\right]=\left[\exp(\underline{\mu}),\exp(\overline{\mu})\right] \]

Normal-Approximation CI for Quantile \(t_p\)

• An approximate \(100(1-\alpha)\%\) CI for \(t_p=\exp[\mu+\Phi^{-1}(p)\sigma]\) based on the normal approximation is

\[ \left[\underline{t}_p,\overline{t}_p\right]=\left[\hat{t}_p/w,\hat{t}_p\times w\right] \]

• where \(w=\exp[z_{(1-\alpha/2)}\widehat{se}_{\hat{t}_p}]\) and

\[ \begin{aligned} \widehat{se}_{\hat{t}_p}&=\sqrt{\widehat{Var}(\hat{t}_p)}=\sqrt{\hat{t}_p^2\widehat{Var}[\log(\hat{t}_p)]}\\ &=\hat{t}_p\left[\widehat{Var}(\hat{\mu})+2\Phi^{-1}(p)\widehat{Cov}(\hat{\mu},\hat{\sigma})+[\Phi^{-1}(p)]^2\widehat{Var}(\hat{\sigma})\right]^{1/2} \end{aligned} \]

\[ \begin{array}{rrrrr} \hline p & Quanhat & Std.Err. & 95\% Lower & 95\% Upper \\ \hline 0.001 & 3116 & 1436 & 1263 & 7689 \\ 0.005 & 5189 & 1797 & 2632 & 10230 \\ 0.010 & 6466 & 1925 & 3608 & 11590 \\ 0.050 & 10830 & 2043 & 7482 & 15680 \\ 0.100 & 13600 & 1981 & 10220 & 18090 \\ 0.200 & 17240 & 1901 & 13890 & 21400 \\ 0.300 & 20000 & 1954 & 16520 & 24220 \\ 0.400 & 22410 & 2141 & 18580 & 27030 \\ 0.500 & 24680 & 2452 & 20320 & 29990 \\ 0.600 & 26960 & 2881 & 21870 & 33240 \\ 0.700 & 29400 & 3446 & 23360 & 36990 \\ 0.800 & 32220 & 4212 & 24940 & 41630 \\ 0.900 & 36090 & 5408 & 26900 & 48410 \\ 0.990 & 44940 & 8615 & 30860 & 65430 \\ \hline \end{array} \]

Normal-Approximate CI for the CDF \(F(t_e)\)

• Define \(t_e\) as a time at which an estimate of the failure probability \(F(t)\) is desired

• An approximate \(100(1-\alpha)\%\) CI for \(\hat{F}(t_e)\) based on \(Z_{\hat{F}}=\frac{\hat{F}(t_e)-F(t_e)}{\widehat{se}_{\hat{F}}}\sim NOR(0,1)\)

\[ \left[\underline{F}(t_e),\overline{F}(t_e)\right]=\hat{F}(t_e)\pm z_{(1-\alpha/2)}\widehat{se}_{\hat{F}} \]

• where

\[ \widehat{se}_{\hat{F}}=\frac{\phi(\hat{\zeta}_e)}{\hat{\sigma}}\left[\widehat{Var}(\hat{\mu})+2\hat{\zeta}_e\widehat{Cov}(\hat{\mu},\hat{\sigma})+\hat{\zeta}_e^2\widehat{Var}(\hat{\sigma})\right]^{1/2} \]

• and

\[ \hat{\zeta}_{e}=[\log(t_e)-\hat{\mu}_{_{MLE}}]/\hat{\sigma}_{_{MLE}} \]

• As discussed in Chapter 3, assuming \(Z_{\hat{F}}=\frac{\hat{F}(t_e)-F(t_e)}{\widehat{se}_{\hat{F}}}\sim NOR(0,1)\) can result in a CI \(\notin (0,1)\) for small sample sizes

• To improve the coverage probability the following transformation was made

\[ \text{logit}[F(t_e|\hat{\mu},\hat{\sigma})]=\log\left[\frac{F(t_e|\hat{\mu},\hat{\sigma})}{1-F(t_e|\hat{\mu},\hat{\sigma})}\right] \]

• where

\[ Z_{\text{logit}[\hat{F}]}=\frac{\text{logit}[\hat{F}(t_e)]-\text{logit}[F(t_e)]}{\widehat{se}_{\text{logit}[\hat{F}]}}\sim NOR(0,1) \]

• An inverse logit transformation is then used to give an approximate \(100(1-\alpha)\%\) CI for \(F(t_e)\) as

\[ \left[\underline{F}(t_e),\overline{F}(t_e)\right]=\left[\frac{\hat{F}}{\hat{F}+(1-\hat{F})\times w}, \frac{\hat{F}}{\hat{F}+(1-\hat{F})/w}\right] \]

• where

\[ w=\exp\left[\frac{z_{(1-\alpha/2)}\widehat{se}_{\hat{F}}}{\hat{F}(1-\hat{F})}\right] \]

• The endpoints of this interval will always be \(\in (0,1)\)

\[ \begin{array}{rrrrr} \hline miles & Fhat & Std.Err. & 95\% Lower & 95\% Upper \\ \hline 6000 & 0.0079 & 0.0078 & 0.0011 & 0.0540 \\ 8000 & 0.0195 & 0.0154 & 0.0041 & 0.0893 \\ 10000 & 0.0391 & 0.0248 & 0.0112 & 0.1322 \\ 12000 & 0.0685 & 0.0352 & 0.0246 & 0.1826 \\ 14000 & 0.1091 & 0.0457 & 0.0471 & 0.2412 \\ 16000 & 0.1615 & 0.0559 & 0.0804 & 0.3092 \\ 18000 & 0.2255 & 0.0660 & 0.1244 & 0.3883 \\ 20000 & 0.2999 & 0.0770 & 0.1769 & 0.4794 \\ 22000 & 0.3823 & 0.0897 & 0.2342 & 0.5810 \\ 24000 & 0.4697 & 0.1039 & 0.2926 & 0.6871 \\ 26000 & 0.5582 & 0.1179 & 0.3499 & 0.7876 \\ 28000 & 0.6439 & 0.1291 & 0.4048 & 0.8718 \\ 30000 & 0.7231 & 0.1349 & 0.4569 & 0.9328 \\ \hline \end{array} \]

Normal-Approximation CI for the Hazard Function \(h(t_e)\)

\[ \begin{array}{rrrrr} \hline Time & Hazard & StdError & 95\% lower & 95\% upper \\ \hline 6000 & 0.0000042 & 0.0000032 & 0.0000009 & 0.0000191 \\ 8000 & 0.0000078 & 0.0000046 & 0.0000025 & 0.0000246 \\ 10000 & 0.0000126 & 0.0000057 & 0.0000052 & 0.0000308 \\ 12000 & 0.0000187 & 0.0000068 & 0.0000091 & 0.0000383 \\ 14000 & 0.0000261 & 0.0000082 & 0.0000140 & 0.0000485 \\ 16000 & 0.0000348 & 0.0000105 & 0.0000193 & 0.0000628 \\ 18000 & 0.0000449 & 0.0000141 & 0.0000242 & 0.0000830 \\ 20000 & 0.0000563 & 0.0000194 & 0.0000287 & 0.0001105 \\ 22000 & 0.0000692 & 0.0000264 & 0.0000327 & 0.0001464 \\ 24000 & 0.0000835 & 0.0000354 & 0.0000364 & 0.0001918 \\ 26000 & 0.0000993 & 0.0000464 & 0.0000397 & 0.0002481 \\ 28000 & 0.0001165 & 0.0000594 & 0.0000429 & 0.0003165 \\ 30000 & 0.0001353 & 0.0000746 & 0.0000459 & 0.0003984 \\ \hline \end{array} \]