Estimating failure-point when failed datapoints are unobservable

Question

I have several datasets of two variables, y and time. Datapoints are objects subject to physical wear, and wear-rates are (roughly) constant over the life of the objects. y is defined as the proportion of the object's mass that has been lost to wear.

At a certain amount of wear, the objects fail and become unobservable. I expect the failure point to be fairly consistent in terms of y - i.e., I expect most objects to fail when they have lost a given proportion of their mass, +/- a negligible amount.

I wish to estimate two parameters as accurately as possible:

1) The mean rate at which objects lose mass to wear

2) The proportion of mass loss at which objects fail and become unobservable

For large datasets, a very coarse method to estimate these parameters would be:

for 1): trim the dataset by time, keeping only arbitrarily 'young' objects. Assume that no wear-related failures have occurred in the population of young objects, and fit a linear regression between mass and time to estimate mean wear-rate.

for 2): estimate the amount of wear at which objects experience wear-rated failure as ~ max(y)

However, this method is likely to return useless estimates for smaller datasets and in datasets where the wear-rate is very slow or the sampled time-range is narrow (i.e., where max(y) is not determined by wear-related failures but by normal wear with few or no failures).

I briefly considered fitting a piecewise regression model with 1) estimated as the slope of the first segment and 2) estimated as the y value at the break-point, but the assumptions of piecewise regression are not met by my datasets - there are not two linear processes forming distinct sections, though the existence of only slow-wearing objects could generate a very similar pattern in the data.

Is there a better way that I could estimate these parameters?

Edit:

I have written a simulation of the processes underlying my dataset for illustration. Here's an example dataset affected by wear-related loss, with a lowess smoother line (f = 0.1) in red:

In this case, wear-related failure occurred when the objects had 55% of their original mass remaining. Because wear-related failure occurs to the fastest-wearing objects first, this creates an elbow in the smoothed line. Before the elbow, the slope is a good approximation of the mean wear-rate for all objects (parameter 1), and after the elbow, the slope flattens fairly dramatically just above the wear-related failure point (parameter 2). This was why I initially considered a piecewise model to estimate the parameters, though the problems with using a piecewise approach are clear in this figure (most notably, the estimates of Parameter 2 would be substantially upwardly-biased).

The following R code generates the data and figure:

points <- 1000

time <- runif(points, 0, 20000)  ## anything up to 20,000 days
wearRates <- rnorm(points, (0.022/365), (0.0022/365))  
## typical wear-rates 
wear <- 1-(time*wearRates) 
losspoint <- 0.55  ## failure when 55% of initial mass remains
observable <- wear > losspoint

allPoints <- data.frame(time, wear, observable)
data <- subset(allPoints, observable == TRUE)

with(data, plot.default(wear ~ time, ylim = c(0,1), ylab="Proportion mass remaining"))
lines(lowess(data$time, data$wear, f = 0.1), type = "l", col=2, lwd = 2)

Answer 1

Before I describe how I approach the problem, I would make a small change in your setup. I would replace $w = 1 - r\,t$ with $w = e^{-r\,t}$, where $w$ is what you call wear and $r$ is the wear rate. (If you don't like that, I think it's easy to follow my approach with your setup.)

The data are of the form $y_{1:n} = (y_1,\ldots,y_n)$ where $y_i = (w_i,t_i)$. You have suggested the unobserved wear rates are distributed $r_i \stackrel{\text{iid}}{\sim} \textsf{N}(\mu,\sigma^2)$. Given $w_i = e^{-r_i\, t_i}$, we have $w_i \stackrel{\text{iid}}{\sim} \textsf{LN}(-\mu\,t_i, \sigma^2\,t_i^2)$, where $\textsf{LN}$ denotes the log-normal distribution.

However, as you point out, we only observe data for which $w_i > b$, where $b$ is the bound you describe. Therefore the actual observations have a truncated log-normal distribution: $\textsf{LN}_{[b,1]}(-\mu\,t_i, \sigma^2\,t_i^2)$. (I have included a truncation from above at 1 as well; it probably doesn't matter.)

Let $\theta = (\mu,\sigma,b)$ denote the unknown parameters. The likelihood is given by $p(y_{1:n}|\theta) = \prod_{i=1}^n p(y_i|\theta)$ where $$ p(y_i|\theta) = \textsf{LN}_{[b,1]}\big(w_i\,\big|\,-\mu\,t_i,\, \sigma^2\,t_i^2\big). $$

Given the likelihood, my approach is Bayesian, in part because the distribution for $b$ is nonstandard and has a sharp cutoff at $\min(w_1,\ldots,w_n)$. The posterior distribution is proportional to the likelihood times the prior distribution $p(\theta)$: $$ p(\theta|y_{1:n}) \propto p(y_{1:n}|\theta)\,p(\theta). $$ Any information you have about the unknown parameters $(\mu,\sigma,b)$ should be incorporated into your prior distribution. As I have no subject-matter expertise, I adopted the flat prior $p(\theta) \propto 1$ to experiment with. This prior is improper, but doesn't seem to produce any problems here. (One has to be careful.)

At this point there are a number of ways to proceed. I chose to adopt a simple random-walk Metropolis sampler to make MCMC draws from the posterior distribution. I use my own code written in Mathematica, so I'm not sure how helpful it would be. The main trick for me is to find a reasonable scale for the proposals, which is not too hard in this case because there are only three parameters and it turns out they are not highly correlated in the posterior

You may be able to code this up in JAGS or Stan and run it from within R. I don't know off hand.

Good luck.

Answer 2

Part 1: Assuming that you are not close to the wear = 0.55 region, the wear at a particular time can be modeled as a normal distribution with non-stationary mean and standard deviation:

$w(t) \sim \mathcal{N}(.|1-st, t\sigma)$

Assuming that your data $(t_i, w_i)$ are independent, you can get MLE estimates for $s$ and $\sigma$ as follows:

$L = \prod_{i=1}^N \mathcal{N}(w_i|1-st_i, t_i\sigma)=\dfrac{1}{(2\pi \sigma^2 \prod_i t_i^2)^{N/2}} \exp \left[ -\dfrac{1}{2\sigma^2} \sum_{i=1}^N \left( \dfrac{w_i-1+st_i}{t_i}\right)^2 \right]$

The loss-function is then:

$\mathcal{L}=-\ln L= \dfrac{N}{2}\ln2\pi+\dfrac{N}{2} \ln \sigma^2 +N\sum_i \ln t_i + \dfrac{1}{2\sigma^2} \sum_{i=1}^N \left( \dfrac{w_i-1+st_i}{t_i}\right)^2$

Partial derivative for the MLE estimates:

$\dfrac{\partial \mathcal{L}}{\partial s}=0 \Rightarrow \hat{s}=\dfrac{1}{N} \sum_i \left( \dfrac{1-w_i}{t_i} \right)$

$\dfrac{\partial \mathcal{L}}{\partial \sigma^2}=0 \Rightarrow \hat{\sigma^2}=\dfrac{1}{N} \sum_{i=1}^N \left( \dfrac{w_i-1+\hat{s}t_i}{t_i}\right)^2$

Here, $\hat{s}$ should be a good estimate of your wear rate. $t\hat{\sigma}$ will give you the spread of wear rates to expect at any time.

(Note: I assume $N$ is large so that you don't have to worry about unbiasing your estimate of $\hat{\sigma})$.

Part 2:

I assume you should have at least a few full data sets where the failures occur.

Say you have $M$ datasets, both small and full. Plot a histogram of the max recorded wears (your $y_{\max}$) over the $M$ datasets. There should be a cluster at the end corresponding to datasets that ran to failure. You can crop out the datasets that don't belong to that cluster and fit a normal distribution to the $y_{\max}$ set of that cluster.

(I say normal distribution since you mentioned:

I expect most objects to fail when they have lost a given proportion of their mass, +/- a negligible amount.

You could fit another, more suitable, distribution that resembles the one you're left with)