I am working on an example in which I have obtained parameter estimates using Stan.

In a real life scenario, I would receive more data every week. The data are covariates of a product which are used as predictors to estimate the failure-time of the product. I would also like to put more weighting on the local data since it is a better indicator of the current use-rate/use-environment of the product than the past data.

I would like to update my parameters every week to make inferences about the failure-time for the product. I have read that I have two options:

(1) Approximate the MCMC posterior in a mathematical form and then use this as the prior for the next update.

(2) Run the MCMC code with all of the data $$\{y_0, y_1\}$$ and use a power prior, where $$y_0$$ is the past data and $$y_1$$ is the data obtained one week from now. The likelihood can be written as $$L(\theta \mid y_0)^aL(\theta \mid y_1)$$, where $$a$$ is a parameter to be estimated.

Option (1) would be very difficult with the number of correlated parameters in the model.

Option (2) seems like an okay option. However, suppose I now have the following data $$\{y_0, y_1,\dots, y_n\}$$. I would, once again, like the most recent observations to have more weighting in the likelihood function, but I do not want to give the same weight to all new observations as $$n$$ gets sufficiently large. Instead, I am imagining the likelihood function would be

$$\begin{equation} L(\theta \mid y_0)^{a_1}L(\theta \mid y_1)^{a_2}\dots L(\theta \mid y_{n-1})^{a_n}L(\theta \mid y_n). \end{equation}$$

This way I would be able to assign weights (where the weights are parameters to be estimated) to the observations as they arrive each week. But, this would mean that I have to constantly change my code and add more parameters each week.

I could instead write the likelihood as

$$\begin{equation} L(\theta \mid y_{old})^{a}L(\theta \mid y_{new}). \end{equation}$$

Fixing the number of parameters in the model. But how would I decide when a "new" observation becomes old?

Can anyone recommend an approach that they would take?

Is option 1 feasible with approximation methods?

Edit: The likelihood function is given by

$$\begin{equation} L(\theta_T, \theta_X \mid \text{Data}) = L(\theta_T \mid \text{Failure-time Data, Covariate History}) \times L(\theta_X \mid \text{Covariate History}), \end{equation}$$

where $$\theta_T = (\mu_0, \sigma_0, \beta)$$ and $$\theta_X = (\eta, \sigma_1, \sigma_2, \rho, \sigma)$$.

The expression of the failure-time data conditional on the covariate history is

\begin{equation} \begin{aligned} L&(\theta_T \mid \text{Failure-time Data, Covariate History}) \\ = &\prod_{i=1}^n\{\exp(\beta x_i(t_i)]f_0(u[t_i;\beta, x_i(t_i)], \theta_0)\}^{\delta_i} \times \{1 - F_0(u[t_i;\beta, x_i(t_i)], \theta_0)\}^{1-\delta_i}, \end{aligned} \end{equation}

where $$f_0$$ and $$F_0$$ are the pdf and cdf of a Weibull distribution with shape $$1/\sigma_0$$ and scale $$\exp(\mu_0)$$.

The likelihood function of the covariate history is

$$\begin{equation} L(\theta_X \mid \text{Covariate History}) = \prod_{i=1}^n\int_{w_i}\bigg\{\prod_{t_{ij} \leq t_i} f_{\text{NOR}}[x_i(t_{ij}) - \eta - Z_i(t_{ij})w_i; \sigma^2]\bigg\} \times f_{\text{BVN}}(w_i; \Sigma_{w})dw_i, \end{equation}$$

where, $$f_{\text{NOR}}( \cdot ; \sigma^2)$$ is the pdf of a univariate normal distribution with mean 0 and variance $$\sigma^2$$, and $$f_{\text{BVN}}( \cdot; \Sigma_{w})$$ is the pdf of a bivariate normal distribution with mean 0 and variance-covariance matrix $$\Sigma_{w}$$.

The following model for $$X_i(t_{ij})$$ is used

$$\begin{equation} X_i(t_{ij}) = \eta + Z_i(t_{ij})w_i + \epsilon_{ij}, \end{equation}$$

where $$\eta$$ is the mean, $$Z_i(t_{ij}) = [1, \log(t_{ij})]$$, $$w_i = (w_{0i}, w_{1i})' \sim N(0,\Sigma_{w})$$, $$\epsilon_{ij} \sim N(0, \sigma^2)$$, and

$$\begin{equation*} \Sigma_{w} = \begin{pmatrix} \sigma^2_1& \rho\sigma_1\sigma_2 \\ \rho\sigma_1\sigma_2 & \sigma^2_2 \end{pmatrix}. \end{equation*}$$

Here, $$w_i$$ is the random effect for modelling variability in a unit’s covariate process over time.

Since the model for the covariate process is normal, I coded it using the built in normal log probability density in Stan (I am assuming this is correct but the likelihood function for the covariate history is not quite a product of normal densities, so I may be wrong?).

The code I used for the likelihood is given below:

target += weibull_lpdf(u_obs| 1/sigma0, exp(mu0));
target += Beta*x_DFD;
target += weibull_lccdf(u_cens|1/sigma0, exp(mu0));
target += normal_lpdf(y_tt|Mu, sig);


Here $$u_{obs}$$ are the values of $$u$$ for the products that failed, $$x_{DFD}$$ is the value of the covariate at the failure time, $$u_{cens}$$ are the values of $$u$$ for the products that are still working, and $$y_{tt}$$ is the value of the covariate $$x_i(t_{ij})$$.

I hope this makes sense. Further description can be found in https://www.tandfonline.com/doi/abs/10.1080/00401706.2013.765324

I also note and acknowledge that the model choice was made by the authors. I am using this model for future work.

• Depending upon the likelihood function, it may be that $\alpha$ can be incorporated into the model cleanly (consider an exponential smoothing model, which effectively downweights past observations in a manner similar to what you are suggesting). Can you be more specific in this regard? Aug 2, 2019 at 18:29
• @jbowman Thank you for your comment. I have edited by post.
– JLee
Aug 2, 2019 at 22:04

From my understanding, option 1 or using a power prior are very similar in the sense that the power prior uses a weighted version of the previously obtained posterior as a new prior. Therefore, the power prior requires that your posterior as a nice mathematical form so that you can approximate it with a known distribution, but this is not always the case unless you are using conjugate priors.

I am not entirely convinced with option 2 since you plan to train the model on all your data and therefore you would lose the apparent benefit of using a power prior. Moreover, if we forget the power prior for a minute, you would expect your model to be valid whether you are using 10 or 20 observations, and by that I mean that you implicitly assume stationarity, i.e. the parameters shouldn't change. In this context, to me, $$a$$ is more a decision parameter and not something to be inferred.

If you want to use all the data, I would instead suggest that you explicitly model that the parameters are time-varying, then you could use a simple exponential smoothing and infer the corresponding smoothing factor.

For instance, let's consider the parameter $$\theta(t)$$ that you would define in the model block and it's smoothed version $$\theta_s(t)$$ that you would define in the transformed parameters block as: $$\theta_s(t) = \alpha \theta(t) + (1 - \alpha) \theta_s(t - 1)$$ Where $$\alpha$$ is the smoothing factor, to be inferred, which you can relate to the time constant $$\tau$$ of the learning process by $$\alpha = 1 - e^{\frac{-\Delta T}{\tau}}$$. You can model the evolution of $$\theta(t)$$ in the model block with: $$\theta(t + 1) \sim \mathcal{N}(\theta_s(t), \sigma^2)$$

Using this kind of approach, you will not only be able to infer when "a parameter becomes old" (cf. $$\alpha$$ or $$\tau$$) but also how much "a parameter is expected to change (cf. $$\sigma$$).

• I am treating $a$ as an unknown parameter. This means that I need to evaluate a normalisation constant. This requires numeric integration. Do you have any suggestions about how to calculate this? The normalisation constant is described here: dx.doi.org/10.1002/sim.3722
– JLee
Aug 3, 2019 at 18:45
• Thanks for the link, I will edit my answer in light of this. Aug 5, 2019 at 8:17