Parameter Updates in Stan

Question

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.

Answer 1

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$).