hi all.
I would like to extend the LONGITUDINAL RANDOM EFFECTS MODELS WITH MEASUREMENT ERROR 1) to a Bayesian framework.
In the following master thesis 2), it is actually done and I want to reproduce it, but it does not converge well.
Could you please help me? please…
-
Buonaccorsi, John, et al. “ESTIMATION IN LONGITUDINAL RANDOM EFFECTS MODELS WITH MEASUREMENT ERROR.” Statistica Sinica , vol. 10, no. 3, 2000, pp. 885–903. JSTOR , www.jstor.org/stable/24306753. Accessed 8 July 2021.
http://www3.stat.sinica.edu.tw/statistica/oldpdf/A10n311.pdf -
Yuan, Zheng . Linear Mixed Effects Models with Measurement Error: Bayesian Approach.
https://prism.ucalgary.ca/bitstream/handle/11023/935/ucalgary_2013_yuan_zheng.pdf;jsessionid=1F0C071F0C076887AA5EB52BFA599B58?sequence=6
data {
int N; // num of obs (instance * year)
int J; // num of instance
int Y; // num of year
int year[N]; // year vector
vector[N] w; // error in variable
vector[N] y; // outcome
matrix[5, 2] A; // desigh matrix
}
parameters {
real beta;
real alpha;
vector[1] epsiron;
real<lower=0> s_epsiron;
vector[Y] x;
vector[2] alpha_x;
vector<lower=0>[2] phy_x;
corr_matrix[2] corr_phy;
vector[1] me;
real <lower=0> s_me;
}
transformed parameters {
cov_matrix[2] cov_phy;
matrix[Y, Y] x_cov;
matrix[Y, Y] w_cov;
matrix[Y, Y] y_cov;
cov_matrix[Y] me_cov;
cov_matrix[Y] epsiron_cov;
vector[2] x_variables;
cov_phy = quad_form_diag(corr_phy, phy_x);
x_cov = A * cov_phy * A';
w_cov = x_cov + me_cov;
y_cov = beta^2 * x_cov + epsiron_cov;
x_variables = alpha_x + phy_x;
me_cov = s_me * diag_matrix(rep_vector(1,Y));
epsiron_cov= s_epsiron * diag_matrix(rep_vector(1,Y));
}
model {
corr_phy ~ lkj_corr(2);
phy_x ~ cauchy(0, 2.5);
me~ normal(0, s_me);
epsiron~ normal(0, s_epsiron);
x ~ multi_normal(A * x_variables, x_cov); //exposure model
w ~ multi_normal(x[year], w_cov[year, year]); //measurement error model
y ~ multi_normal(alpha + x[year] * beta, y_cov[year, year]); //disease model
}
df = read.csv('sample_1d_data.csv')
a_vec = c(1, 1, 1, 1, 1, 0, 1, 2, 3, 4)
A = matrix(a_vec, nrow=5)
stanfile <- 'stancode.stan'
data <-list(
N=length(df$id),
J=max(df$id),
w=df$w,
year=df$year,
y=df$y,
A=A,
Y=max(df$year),
)
error code
Chain 1: Rejecting initial value:
Chain 1: Error evaluating the log probability at the initial value.
Chain 1: Exception: multi_normal_lpdf: LDLT_Factor of covariance parameter is not positive definite. last conditional variance is -6.66134e-16. (in 'model354c455069e1_stancode_matrix' at line 49)
Chain 1:
Chain 1: Initialization between (-2, 2) failed after 100 attempts.
Chain 1: Try specifying initial values, reducing ranges of constrained values, or reparameterizing the model.
[1] "Error in sampler$call_sampler(args_list[[i]]) : Initialization failed."