Accurate higher-order likelihood inference on $P(Y<X)$

Abstract

The stress-strength reliability $R=P(Y<X)$, where $X$ and $Y$ are independent continuous random variables, has obtained wide attention in many areas of application, such as in engineering statistics and biostatistics. Classical likelihood-based inference about $R$ has been widely examined under various assumptions on $X$ and $Y$. However, it is well-known that first order inference can be inaccurate, in particular when the sample size is small or in presence of unknown parameters. The aim of this paper is to illustrate higher-order likelihood-based procedures for parametric inference in small samples, which provide accurate point estimators and confidence intervals for $R$. The proposed procedures are illustrated under the assumptions of Gaussian and exponential models for $(X,Y)$. Moreover, simulation studies are performed in order to study the accuracy of the proposed methodology, and an application to real data is discussed. An implementation of the proposed method in the R software is provided.

Appendices

Appendix A: Computation of statistics for the Gaussian and exponential models

Under the assumption of Gaussian distributions with equal variances for $Y$ and $X$ as described in Sect. 3.1, in order to compute the Wald statistic, the profile observed information is obtained from (5). Its expression is

$$\begin{aligned} j_p(\hat{\psi }) = \frac{\hat{G}}{(n\, \hat{\lambda }_2)^2 \, +\, n_1\, n_2\, (\hat{\lambda }_1 \hat{\lambda }_2\, -\, D(\hat{\theta }))^2} \, , \end{aligned}$$

(17)

where $\hat{G} =2 n n_1 n_2 \hat{\lambda }_2^2 \hat{d}^2$, with $\hat{d} = \partial \varPhi ^{-1}(\hat{\psi })/\partial \hat{\psi }= 1/\phi (\varPhi ^{-1}(\hat{\psi }))$ and $\phi (\cdot )$ standard normal probability density function.

In the same example, for computing the adjustment term $q$ of the $r_p^*$ statistic, given in (7), we have

$$\begin{aligned} \frac{ |\ell _{\lambda ;\hat{\lambda }}(\hat{\theta }_\psi )|}{(|j(\hat{\theta })| |j_{\lambda \lambda }(\hat{\theta }_\psi )|)^{1/2}} = \frac{2 \hat{\lambda }_2^2 \left(2 n_1 n_2 \hat{B}^2 +(n \hat{\lambda }_2)^2 + n_1 n_2 \hat{B} (\tilde{\lambda }_1 \tilde{\lambda }_2 - D(\hat{\theta }_{\psi })) \right) }{\tilde{\lambda }_2^2\, \left( \hat{G}\, (12 n_1 n_2 \hat{B}^2 + 8\hat{C}^2 - 2 n^2 (\tilde{\lambda }_2^2 - 3\hat{\lambda }_2^2) - 8n \hat{C}) \right)^{1/2} } \, , \end{aligned}$$

with $\hat{B} = \hat{\lambda }_1 \hat{\lambda }_2 -D(\hat{\theta })$, $\hat{C} = n_1 \hat{\lambda }_1 \hat{\lambda }_2 + n_2 D(\hat{\theta })$ and $\hat{F} = n_1\tilde{\lambda }_1 \tilde{\lambda }_2 \hat{\lambda }_1 \hat{\lambda }_2 + n_2 D(\hat{\theta }) D(\hat{\theta }_{\psi })$. The tilde above the parameters indicates an estimate constrained to a fixed $\psi$, while the hat indicates an MLE estimate. The sample space derivatives needed for computing the $q$ term are

$$\begin{aligned} \ell _{;\hat{\psi }}(\hat{\theta })=0, \ \quad \ell _{;\hat{\lambda }}(\hat{\theta })=\left( \ 0, \ -\frac{n}{\hat{\lambda }_2} \ \right), \ \quad \ell _{;\hat{\psi }}(\hat{\theta }_\psi )= \frac{2 n_2 \hat{\lambda }_2}{\tilde{\lambda }_2^2} K [ \tilde{\lambda }_2 b - D(\hat{\theta })], \end{aligned}$$

with $K= \partial \varPhi ^{-1} (\hat{\psi }) / \partial \hat{\psi }$.

The sample space derivatives $\ell _{;\hat{\lambda }}(\hat{\theta }_\psi ) = ( \ell _{;\hat{\lambda }_1}, \ \ell _{;\hat{\lambda }_2} )$ are computed as

$$\begin{aligned}&\displaystyle \ell _{;\hat{\lambda }_1}(\hat{\theta }_\psi ) = \frac{-2 \hat{\lambda }_2}{\tilde{\lambda }_2^2} \left[n_1 \hat{\lambda }_1 \hat{\lambda }_2 + n_2 D(\hat{\theta }) - n_1 \tilde{\lambda }_1 \tilde{\lambda }_2 - n_2 \tilde{\lambda }_2 b\right],\\&\displaystyle \ell _{;\hat{\lambda }_2}(\hat{\theta }_\psi ) = \frac{-1}{\tilde{\lambda }_2^2} \left[n \hat{\lambda }_2 + 2 n_1 \hat{\lambda }_1^2 \hat{\lambda }_2 + \frac{2 n_2 D(\hat{\theta })^2}{\hat{\lambda }_2} - 2 n_1 \tilde{\lambda }_1 \tilde{\lambda }_2 \hat{\lambda }_1 - \frac{2 n_2 \tilde{\lambda }_2 b D(\hat{\theta })}{\hat{\lambda }_2}\right], \end{aligned}$$

while $\ell _{\lambda ;\hat{\psi }}(\hat{\theta }_\psi ) = ( \ell _{\lambda _1;\hat{\psi }}, \ \ell _{\lambda _2;\hat{\psi }} )$ are given as

$$\begin{aligned} \ell _{\lambda _1;\hat{\psi }}(\hat{\theta }_\psi )=\frac{2 n_2 \hat{\lambda }_2 K}{\tilde{\lambda }_2}, \ \quad \ell _{\lambda _2;\hat{\psi }}(\hat{\theta }_\psi )= \frac{2 n_2 \hat{\lambda }_2 K}{\tilde{\lambda }_2} \left[\frac{2 D(\hat{\theta })}{\tilde{\lambda }_2} - b\right], \end{aligned}$$

with $b = \varPhi ^{-1} (\psi ) + \tilde{\lambda }_1$.

The quantity $\ell _{\lambda ;\hat{\lambda }}(\hat{\theta }_\psi )$ is equal to the $(d-1)$ square matrix

$$\begin{aligned} \begin{bmatrix} \ell _{\lambda _1;\hat{\lambda }_1}&\quad \ell _{\lambda _1;\hat{\lambda }_2} \\ \ell _{\lambda _2;\hat{\lambda }_1}&\quad \ell _{\lambda _2;\hat{\lambda }_2}. \\ \end{bmatrix} \end{aligned}$$

The sample space derivatives therein can be computed from the full likelihood function

$$\begin{aligned} \ell (\theta )&= \ell (\psi ,\lambda _1,\lambda _2)\\&= -n log \lambda _2 - \frac{1}{\lambda _2^2} \left[\frac{n \hat{\lambda }_2^2}{2} + n_1 (\lambda _1 \lambda _2 - \hat{\lambda }_1 \hat{\lambda }_2)^2 + n_2 (D(\hat{\theta }) - D(\theta ))^2 \right], \end{aligned}$$

which only depends on the MLE estimates and the parameter $\theta$. Thus, we have

$$\begin{aligned}&\displaystyle \ell _{\lambda _1;\hat{\lambda }_1}(\hat{\theta }_\psi )= \frac{2n \hat{\lambda }_2}{\tilde{\lambda }_2} \ \qquad \qquad \ell _{\lambda _1;\hat{\lambda }_2}(\hat{\theta }_\psi )= \frac{2}{\tilde{\lambda }_2} \left[n_1 \hat{\lambda }_1 + n_2 \frac{D(\hat{\theta })}{\hat{\lambda }_2}\right]\\&\displaystyle \ell _{\lambda _2;\hat{\lambda }_1}(\hat{\theta }_\psi )= \frac{2 \hat{\lambda }_2}{\tilde{\lambda }_2^3} \left[2 n_1 \hat{\lambda }_1 \hat{\lambda }_2 + 2 n_2 D(\hat{\theta }) - n_2 b \tilde{\lambda }_2 - n_1 \tilde{\lambda }_1 \tilde{\lambda }_2 \right]\\&\displaystyle \ell _{\lambda _2;\hat{\lambda }_2}(\hat{\theta }_\psi )= \frac{2}{\hat{\lambda }_2 \tilde{\lambda }_2^3} \left[n \hat{\lambda }_2^2 + 2n_1 \hat{\lambda }_1^2 \hat{\lambda }_2^2 + 2n_2 D(\hat{\theta })^2 - n_1 \tilde{\lambda }_1 \tilde{\lambda }_2 \hat{\lambda }_1 \hat{\lambda }_2 - n_2 b \tilde{\lambda }_2 D(\hat{\theta }) \right]. \end{aligned}$$

Appendix B: Implementation in R

We present the R code for the estimation of $R$ under a parametric setting. Continuous quantities of interest can be assumed as exponential variables or Gaussian variables with either equal or different variances. The R package ProbYX can be downloaded from the webpage http://homes.stat.unipd.it/gcortese under the heading ‘Software’.

The two data sets from two different populations are called, respectively, ydat and xdat. MLEs for the parameter of interest $\psi$ and the nuisance parameter $\lambda$ can be obtained from the following code

MLEs(ydat,xdat,distr),

where distr can be set equal to either "exp", "norm_EV" or "norm_DV", according as the distributions assumed for the continuous markers are exponential or Gaussian with equal or unequal variances, respectively. The log-likelihood can be computed for a given set of values of $\psi$ and $\lambda$ by means of the function

loglik(ydat,xdat,lambda,psi,distr),

where lambda and psi are, respectively, $\lambda$ and $\psi$. For the case of Gaussian distributions with different variances the following simpler reparameterisation has been used

$$\begin{aligned} \psi = R(\theta ) = \varPhi \left( \frac{\mu _2 - \mu _1}{\sqrt{\sigma _1^2 + \sigma _2^2}} \right) \qquad \text{ and} \qquad \lambda =(\lambda _1,\lambda _2,\lambda _3) = (\mu _1, \sigma _1^2, \sigma _2^2) \, , \end{aligned}$$

with $\sigma _1^2$, $\sigma _2^2$ variances of the Gaussian variables $Y$ and $X$, respectively.

Point estimates and confidence intervals for $R$ are obtained by using the R function Prob as follows:

Prob(ydat,xdat,distr,method,level),

where the argument method, set to be equal to either “Wald”, “RP” or “RPstar”, allows to choose confidence intervals based on the Wald, $r_p$ or $r_p^*$ statistic, respectively (see formula (3) and (4) and Sect. 3). When the methods “Wald” or “RP” are selected, point estimate for $\psi$ corresponds to the MLE, while the method “RPstar” yields the $r_p^*$-based point estimate for $\psi$ as presented in Sect. 3. The confidence level $(1-\alpha )$ can be decided by setting level $=\alpha$.

The Wald, $r_p$ and $r_p^*$ statistics can also be computed for a given value of $\psi$ by applying, respectively, the R functions

• ${\mathtt {w <- wald(ydat,xdat,psi,distr)} },$

• ${\mathtt {r <- rp(ydat,xdat,psi,distr)} },$

• ${\mathtt {r\_star <- rpstar(ydat,xdat,psi,distr)} }$.

The steps 4–6, given at the end of Sect. 3, for application of the higher-order procedures based on the signed log-likelihood ratio statistic $r_p^*$, can be implemented by means of the following R commands:

• ${\mathtt {smoother <- smooth.spline(V.r\_star,r\_star.range)} }$

• ${\mathtt {psi1 <- predict(smoother,z1)\$y} }$

• ${\mathtt {psi2 <- predict(smoother,z2)\$y} }$

• ${\mathtt {hatpsi <- predict(smoother,0)\$y} }$

where V.r_star is the vector of values of $r_p^*$, r_star, calculated by applying the R function rpstar repeatedly on an appropriate range r_star.range of values of $\psi$, psi1 and psi2 are the limits $\psi ^*_1$ and $\psi ^*_2$ of the confidence interval corresponding to the percentiles z1 $=z_{1-\alpha /2}$ and z2 $=z_{\alpha /2}$, and hatpsi is the point estimate $\hat{\psi }^*$.