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.