If we were working with normal distributed data wouldn’t be any problem. But now there is a non linear transformation that in multivariate data can be something to pay attention.
The exp family DLM can be programmed in Stan, I think I have some preliminary codes somewhere, and can be used to compare VAR models.
Great. Can we apply your code in the data I uploaded here to see if it works?
Thanks
@asael_am
Hello,
I ran your code.
A couple of questions:
1- What I understood is the code is a code for a repeated measure for one person ( there was no looping for the number of individuals).
2- The data is in long format, right?
So with this token, for a generic data as below
id time e hr bon
1 1 1 2 8
1 2 2 3 11
1 3 3 3 11
1 4 5 4 10
offset <- 2 # column number we ignore to start with in data.frame
nth <- as.integer(11)
data_bon <-list (
N = nrow(data_csv),
D = 2,
y = as.array(data_csv[, "bon"]),
p = 3,
nthres = nth,
x = as.matrix(data_csv[, (1 + offset):(offset + 2)], rownames.force = NA)
)
data {
int<lower=2> nthres; // number of thresholds
int<lower=0> N; // number of data items
int<lower=1> D; // number of regressor variables
int<lower=0> p; // number of autoregressive coefficents
matrix[N,D] x; // regressor data in a matrix
int<lower=1,upper=nthres> y[N];// outcome time series
}
parameters {
real theta0;
vector[D] beta;
ordered[nthres] Intercept;
vector<lower=-1,upper=1>[p] phi;
}
transformed parameters{
// form parameter;
vector[N] theta;
theta = x * beta;
for (i in 1:N) {
// constant in autorregresive
theta[i] += + theta0;
// Autoregressive part
if(p > 0) for(j in 1:p) if(i > j) theta[i] += y[i-j]*phi[j];
}
}
model {
target += normal_lpdf(phi|0,0.5);
target += normal_lpdf(beta| 0,1);
target += normal_lpdf(Intercept | 0,10);
target += ordered_logistic_lpmf(y|theta, Intercept);
}
and the results are
mean se_mean sd 2.5% 25% 50%
theta0 6.05083793 0.18218977 4.7225481 -2.45601200 2.6989507 5.88226463
beta[1] 0.05891928 0.03595626 0.9972350 -1.91102536 -0.6256622 0.07992466
beta[2] -0.22489832 0.03413903 0.9243351 -1.93378050 -0.8606048 -0.20425814
Intercept[1] -15.82924254 0.26063241 5.6535578 -28.29973037 -19.3954881 -15.45379435
Intercept[2] -10.56605088 0.12927146 4.4399404 -19.31396617 -13.3723420 -10.45228431
Intercept[3] -7.29034039 0.10717417 4.1156654 -15.43819591 -10.0767843 -7.38588243
Intercept[4] -4.67818672 0.10796824 3.8041731 -12.19126058 -7.3074708 -4.56105732
Intercept[5] -2.36494641 0.09727661 3.6234013 -9.15883921 -4.8871701 -2.30541060
Intercept[6] -0.03352899 0.10584216 3.6771769 -7.15484537 -2.4607416 0.00188459
Intercept[7] 2.08107561 0.11679291 3.7372556 -5.11887969 -0.3678580 2.09622741
Intercept[8] 5.68207277 0.13567748 3.8958760 -1.65207017 2.7898973 5.67823170
Intercept[9] 6.98947730 0.12545235 3.9298268 -0.28491536 4.1641930 6.84170598
Intercept[10] 9.55464230 0.12401111 4.0633199 2.00991320 6.5827207 9.63606982
Intercept[11] 17.32928191 0.18058059 5.8664985 6.73634647 13.4729556 16.93757762
phi[1] 0.56543382 0.01177349 0.2685494 -0.02187386 0.3972104 0.59387129
phi[2] -0.07473985 0.01222218 0.3013170 -0.67846314 -0.2695174 -0.07868584
phi[3] -0.21309316 0.01581580 0.3642518 -0.86264280 -0.4797770 -0.22058142
theta[1] 5.65996056 0.13791274 4.1647586 -2.03755919 2.6111322 5.60531775
theta[2] 10.18343109 0.12429771 4.1346106 2.41715481 7.1571874 10.06207776
theta[3] 11.28181377 0.13723356 4.3630684 3.25439068 8.1563613 11.13605511
theta[4] 8.90305231 0.12764299 4.3760935 0.56697446 5.9668145 8.82654073
lp__ -46.82198511 0.21402505 3.5702248 -54.76456947 -48.8048916 -46.33826929
75% 97.5% n_eff Rhat
theta0 9.24229201 15.6236864 671.9003 1.0024450
beta[1] 0.75413159 2.0703808 769.2118 0.9996928
beta[2] 0.40909431 1.5709185 733.0887 1.0001862
Intercept[1] -11.79663717 -6.2759799 470.5295 0.9999744
Intercept[2] -7.60552494 -1.9692638 1179.6385 0.9994127
Intercept[3] -4.31153205 0.4554698 1474.6869 0.9995625
Intercept[4] -2.13647482 2.7373154 1241.4479 0.9991600
Intercept[5] 0.03317037 4.8361339 1387.4457 0.9990185
Intercept[6] 2.48364431 6.6368272 1207.0121 0.9994116
Intercept[7] 4.60798435 9.2595222 1023.9356 0.9997249
Intercept[8] 8.35263168 13.4646941 824.5077 0.9997787
Intercept[9] 9.71929549 15.1773088 981.2716 0.9999141
Intercept[10] 12.19482859 17.7019915 1073.5959 0.9990707
Intercept[11] 20.91592393 29.8830983 1055.3969 1.0003892
phi[1] 0.77509512 0.9750913 520.2810 1.0076334
phi[2] 0.12745142 0.5161664 607.7849 1.0009013
phi[3] 0.02433008 0.5359634 530.4212 0.9992130
theta[1] 8.40285139 13.6804988 911.9497 1.0010187
theta[2] 12.99233049 18.4470858 1106.4784 0.9990081
theta[3] 14.26037077 20.4522761 1010.7957 0.9990170
theta[4] 11.74867580 17.9808806 1175.3826 0.9991229
lp__ -44.26685907 -41.1199757 278.2671 1.0035558
If this is correct, then I need to add a loop for more individuals, and another loop to change the response variables. Right?
Should the data be in a long format?
Thanks,
But if there are only 4 times per individual, then a time series model is not suitable (you need more times) is better to put each time as regressors variables. And then N will be the total of individuals (id).
But seeing the code just now, a multilevel model for each e will do a better job than a time series with super informative priors (because of the low times).
A hierarchical regression or multilevel models can be done in brms with no problems. For this kind of problems I am not the right person to ask, but all Stan Forum can help you super fast (sorry I do time series).
Can you send me a mail with the problem you wanna analyze with this data structure? Seeing your code now make me 70% sure that a time series structure might not do what you really want. :)
Thanks for your response.
I know hoe to do multilevel models.
However, I am using lagged model because I want to know how lagged features predict responses.
You are right about time-series, but I still think a dynamic model helps as 24 features at 4 time points would be response variables. Maybe a dynamic network will work. I want to know how features predict one another as lagged model.
Exactly, similar to the diagram in the first post:
Exactly, similar to the diagram in the first post.
e4 ~ e3 + e2 + e1
h4 ~ e3 + e2 + e1
e4 ~ h3 + h2 + h1
bon4 ~ bon3 + bon2 + bon1
e4 ~ bon3 + bon2 + bon1
bon4 ~ e3 + e2 + e1
bon4 ~ h3 + h2 + h1
bon4 ~ h3 + h2 + h1+ e3 + e2 + e1
In my real data I have 24 features and 4 time points that i need to develop these models, all these correlations includes autoregressive and lagged associations with other features.
a shortened data with 3 variables is in a wide format as :
id e1 hr1 bon1 e2 hr2 bon2 e3 hr3 bon3 e4 hr4 bon4
1 1 2 8 1 2 11 1 2 11 1 4 10
2 2 1 9 2 1 10 2 1 12 2 1 10
3 2 1 8 2 1 10 2 1 10 2 2 10
4 1 1 8 1 1 9 1 1 8 2 1 9
5 5 2 7 5 2 7 5 2 8 5 2 8
6 3 1 6 3 1 5 3 1 6 3 1 6
I hope this is clear.Well in that case adding the loop for id will solve the univariate part!
So now the multivariate part is adding the lag matrix parameters and define the inverse link function.
I tried to follow the whole thread, but I am not sure I understand everything that was proposed, so sorry, if I am repeating something that was already there, but I think what I propose is a slightly different angle than what others were suggesting. Not sure it is better, but wanted to share it.
I also hope I understood the diagram and subsequent posts correctly. But here’s my interpretation (for a single subject):
t, there is some unknown quantity related to eye and body (f_eye[t], f_body[t]).f_eye[t] and f_body[t] via ordinal variables Y_eye[t] and Y_body[t]f_eye[t] is a function of f_eye[t-1], f_body[t - 1] and some other predictors (say a vector covariate X[t, ]) plus potentially some unaccounted other influences that we model as additional gaussian noisef_eye[1] and f_body[1]. We have multiple options: we might give them an explicit prior or we might compute them from the covariates assuming f_eye[0] = f_body[0] = 0. I’ll chose the first option in the following.So if this is correct, you can translate this very directly into a Stan program (just a sketch, didn’t try to compile/run this, but hope it shows the approach well):
data {
int<lower=1> N_time;
int<lower=2> K; //number of observable categories
int<lower=0> G; //number of covariates
int<lower=1,upper=K> Y_eye[N_time];
int<lower=1,upper=K> Y_body[N_time];
matrix[N_time - 1, G] X; // We don't have predictors for the first time point
}
parameters {
vector[G] beta_eye;
vector[G] beta_body;
// Effects from previous state to new state
real b_eye_eye;
real b_body_eye;
real b_eye_body;
real b_body_body;
vector[N_time] f_eye;
vector[N_time] f_body;
real<lower = 0> sigma_eye;
real<lower = 0> sigma_body;
vector[K - 1] c ; // cut-points - assuming same for eye and body
}
model {
//Prior for coefficients - choose your own
beta_eye ~ normal(0, 1);
beta_body ~ normal(0, 1);
b_eye_eye ~ normal(0, 1);
b_body_eye ~ normal(0, 1);
b_eye_body ~ normal(0, 1);
b_body_body ~ normal(0, 1);
sigma_eye ~ normal(0, 1);
sigma_body ~ normal(0, 1);
c ~ student_t(3, 0, 1);
//Prior for initial state - choose your own
f_eye[1] ~ normal(0, 5);
f_body[1] ~ normal(0, 5);
//Propagate the state information
for(t in 2:N_time) {
f_eye[t] ~ normal( X[t - 1, ] * beta_eye + b_eye_eye * f_eye[t - 1] + b_body_eye * f_body[t - 1], sigma_eye);
f_body[t] ~ normal( X[t - 1, ] * beta_body + b_eye_body * f_eye[t - 1] + b_body_body * f_body[t - 1], sigma_body);
}
// Likelihood
Y_eye ~ ordered_logistic(f_eye, c);
Y_body ~ ordered_logistic(f_body, c);
}
If that’s what you are after and it works for a single subject, you can hopefully generalize to multiple subjects (add more dimensions to X, f_, Y_) .
This would be computationally very expensive (you have a new parameter for each observation), but if your dataset isn’t big, it should be workable.
Technical aside: If the distribution for y_eye and y_body was Gaussian, we could get rid of treating f_eye and f_body as parameters using the Kalman filter - there are some implementations in Stan floating around (e.g. Kalman filter in Stan · GitHub). I however don’t think the Kalman filter can be easily extended for non-Gaussian likelihood.
Hope that illuminates more than confuses.
I would slightly disagree with the certainty of the message. While fitting a time series to 4 times can really be problematic in many cases, it is IMHO not implausible. A lot depends on the exact data you have - if they give you a lot of information about the state of the process, it might work.
Best of luck with your model!
Omg! I think you have it @martinmodrak (insert a celebration emoji here!)
Well, what I tried to propose just solve for one dimension (one side of the diagram) cause I didn’t know how to generalized for vectors (how to apply the inverse loglik to multivariate data). This approach of the likelihood to every separate variable looks good and practical.
My only question is here:
for(t in 2:N_time) {
f_eye[t] ~ normal( X[t - 1, ] * beta_eye + …
f_body[t] ~ normal( X[t - 1, ] * beta_body + …;
}
It will not be :
for(t in 2:N_time) {
f_eye[t] ~ normal( X[t , ] * beta_eye + …
f_body[t] ~ normal( X[t , ] * beta_body + …;
}
instead?
Yes this approach is easier to handle.
I don´t have a reason to say no. I will actually trust you more! Thank you for taking your time :)
[quote=“martinmodrak, post:27, topic:15155”]
I tried to follow the whole thread, but I am not sure I understand everything that was proposed, so sorry, if I am repeating something that was already there, but I think what I propose is a slightly different angle than what others were suggesting. Not sure it is better, but wanted to share it.
[/quote
Dear @martinmodrak and @asael_am,
Thank you so much for your input. Your suggestion is very clear and straightforward. The time in my data is not equal for each individual, range from 2 to 8 years intervals.
The cut-point for each variable is different. I assume I should loop over the cut-point values or define it for each variable.
What is your advice?
Thanks a lot for your help.
Shabnam
Oh, that’s tough and makes the problem much harder - it forces you to be more explicit about the process you think about - e.g. how does the influence of previous observations change with increasing intervals?
If the only thing that’s missing are the observations, but you have the covariates(predictors) for each year (or can assume they don’t change between measurements), you can treat that the years you have no observations for as missing data. So you compute f_eye and f_body for each year, but only those that were observed enter the likelihood part. But once again this increases the computational burden substantially and might make the model intractable (if it was tractable before).
Yes, just have two sets of cutpoints - either in [K-1,2] matrix or as two separate variables. (note all the _eye and _body alternatives can be treated as vectors/matrices to simplify the code, I kept it separate to make the code more instructive, but you certainly don’t need to do so).
Yes and no. The way I’ve written the model, we ignore predictors for the first time point, so we have one fewer predictors than times and the X array is shifted by one. But you could build the model so that f_eye[1] depends on X[1], and then you would use all the predictors.
Yes, I am with @martinmodrak not sure if the current model is tractable!
Thanks for your hints and explanation. To be honest, I am not concerned about the variable time-intervals as that is because of the study design not missingness. Thanks again for the code. It is very helpful. I am working on incorporating various cut-points and other covariates as I have 25 variables to loop on.
Thnaks,
Just a quick note - it is usually advisable to start with a small model and see how that behaves, either you can use a subset of the full dataset or (even better) simulated data where you know the true values of all parameters. Only once you have the small model working you should add more complexities. Otherwise it is all too easy to end up with a few hundred lines of code that have an issue somewhere, but it is impossible to debug.
Also in this specific case, testing whether the model doesn’t take ages to fit a dataset of your size is probably sensible before you spend a lot of time on it.
Best of luck!
@martinmodrak
Thanks for all hints. I started with simple code. Your code worked for one person and one variable.
I am adding more variables for one person and still working around looping parts.
I will put the results here to update. Thanks again @martinmodrak and @asael_am for your help.
Hi,
We developed the model with two covariates and three time-points as below, as time-intervals between three time-points are varied for each individual, I also defined time as variable “yrs” as a covariate.
functions {
//this function gets a matrix input all vars except one which is outcome
matrix select_columns(matrix input, int nRows, int nCols, int dropColIndex) {
matrix[nRows, nCols - 1] ret;
for(i in 1:nRows){
for(j in 1:nCols){
if(j < dropColIndex) ret[i, j] = input[i, j];
else if (j > dropColIndex) ret[i, j - 1] = input[i, j];
}
}
return ret;
}
// personIndex is the person in the data with its position
matrix select_person(matrix input, int nTime, int nCols, int personIndex) {
matrix[nTime - 1, nCols] ret;
int startIndex = (nTime * (personIndex - 1)) + 1;
// we need data from startIndex by nTime-1 (lagged) and also subtract the startIndex to move to next person
for(r in startIndex:(startIndex + nTime - 2)){
for(c in 1:nCols){
ret[r - startIndex + 1, c] = input[r, c];// e-startIndex=0 then
}
}
return ret;
}
}
data {
int<lower=1> N_time;
int< lower=1 > N; //total population;
int<lower=0> G; //number of covariates, excluding bodn
int<lower=2> K_G[G]; //number of observable categories
int<lower=2> K_Sum; //sum of observable predictors excluding body
int<lower=2> K_body; //number of observable categories
int<lower=1> Y[N_time * N, G];
int<lower=1,upper=K_body> Y_body[N_time * N];
real<lower=0> G_years[N_time * N];
matrix[N_time * N, G] X; // We don't have predictors for the first time point
}
transformed data {
int<lower=1> Y_tr[N_time, G, N];
int<lower=1> Y_body_tr[N_time, N];
for(t in 1:N_time){
for(n in 1:N){
for(g in 1:G){
Y_tr[t, g, n] = Y[t * n, g];
}
Y_body_tr[t, n] = Y_body[t * n];
}
}
}
parameters {
vector[G - 1] beta[G];
vector[G] beta_body;
real beta_years;
real b[G];
real b_body[G];
real b_body_body;
vector[N_time] f[G, N];
vector[N_time] f_body[N];
real<lower = 0> sigma[G];
real<lower = 0> sigma_body;
// with K_Sum we produce a sum of cut-point sof all G variables and then each time we subtract
// G cutpoints from the sum
ordered[K_Sum - G] c;
ordered[K_body - 1] c_body;
}
model {
int position = 0; //for cut-point selection
for (i in 1:G){
beta[i] ~ normal(0, 1);
b[i] ~ normal(0, 1);
b_body[i] ~ normal(0, 1);
sigma[i] ~ cauchy(0,0.5);
for(n in 1:N){
f[i, n][1] ~ normal(0, 5);
}
}
beta_body ~ normal(0, 1);
beta_years ~ normal(0, 1);
b_body_body ~ normal(0, 1);
sigma_body ~ normal(0, 1);
for(n in 1:N){
f_body[n][1] ~ normal(0, 5);
}
//Propagate the state information
for(t in 2:N_time) {
for (g in 1:G){
for(n in 1:N){
f[g, n][t] ~ normal( select_person(select_columns(X, N_time * N, G, g),
N_time, G - 1, n)[t - 1,] * beta[g] +
b[g] * f[g, n][t - 1] +
b_body[g] * f_body[n][t - 1] + beta_years * G_years[n * t], sigma[g]);
}
}
for(n in 1:N){
f_body[n][t] ~ normal( select_person(X, N_time, G, n)[t - 1, ] * beta_body +
b_body_body * f_body[n][t - 1] + beta_years * G_years[n * t], sigma_body);
}
}
for (g in 1:G){
vector[K_G[g] - 1] cp = c[(position + 1):(position + K_G[g] - 1)];
position += position + K_G[g] - 1;
cp ~ student_t(3, 0, 1); // prior for cut-point
for(n in 1:N){
Y_tr[, g, n] ~ ordered_logistic(f[g, n], cp);
}
}
c_body ~ student_t(3, 0, 1);
for(n in 1:N){
Y_body_tr[, n] ~ ordered_logistic(f_body[n], c_body);
}
}
I appreciate your feedback on this.
The output looks like below:
beta_e_hr 0.050395909 0.6672185 -1.24885134 1.25779763
beta_hr_e 0.512087322 0.6372881 -0.52652329 1.89903187
beta_body_e 0.755302199 0.7433808 -0.51206939 2.34644795
beta_body_hr 0.376309329 0.5188516 -0.55336491 1.42843827
beta_years 0.055748186 0.1765089 -0.25565982 0.37867082
b_e_e 0.247697540 0.4390430 -0.59364752 1.08782187
b_hr_hr 0.341984254 0.4071307 -0.46687935 1.09842227
b_e_body -0.835428520 0.5538218 -2.02667535 0.23969166
b_hr_body -0.520789651 0.7067315 -1.88771354 0.78589185
b_body_body -0.000964786 0.6562512 -1.16369311 1.25138771
f_e_n1_t1 -4.309393434 2.6136138 -9.88533566 0.50273342
f_hr_n1_t1 1.864706667 1.8866074 -1.23764746 6.20032112
f_e_n2_t1 -1.802023365 1.7879316 -5.43283475 1.81219151
f_hr_n2_t1 4.418231434 2.8910659 -0.65202836 10.88769634
f_e_n3_t1 -0.538218374 1.8782298 -3.76102981 3.80766239
f_hr_n3_t1 4.015153954 2.6181371 -0.49278353 9.70408960
f_e_n4_t1 -4.720660029 2.6004952 -10.07332466 -0.38849620
f_hr_n4_t1 -3.671408473 3.1290634 -10.67296449 1.57139916
f_e_n1_t2 -1.034875571 1.5588196 -4.16177197 2.09016662
f_hr_n1_t2 2.735759315 2.8921214 -0.76002867 11.48318328
f_e_n2_t2 -4.418297951 2.5775154 -9.43458532 -0.34621993
f_hr_n2_t2 -0.354805627 2.0865423 -4.47516567 3.68655972
f_e_n3_t2 -3.232843959 2.5141742 -9.53302301 0.35976341
f_hr_n3_t2 1.306214148 1.4669948 -1.44666155 4.39781256
f_e_n4_t2 -1.917104784 1.4182010 -5.17641266 0.66321080
f_hr_n4_t2 -1.554761259 2.1645973 -6.69465708 1.84760814
f_e_n1_t3 -1.122146605 1.5466718 -4.29318127 2.15248253
f_hr_n1_t3 2.906326979 3.0378513 -1.18167060 11.07972257
f_e_n2_t3 -3.072469728 2.9283578 -10.85647339 1.54530070
f_hr_n2_t3 1.181675570 1.5492420 -1.91342947 4.21035962
f_e_n3_t3 -2.100910942 1.7284324 -5.70743849 1.38905440
f_hr_n3_t3 -0.032990289 1.6439964 -3.48733198 2.87561716
f_e_n4_t3 -0.297312014 1.7325890 -3.66683649 3.16950395
f_hr_n4_t3 1.375329315 1.8175420 -2.11799206 5.13158080
f_body_n1_t1 0.409607207 1.4789244 -2.75101440 3.59299201
f_body_n2_t1 3.656919862 2.0129381 0.42547093 8.42773947
f_body_n3_t1 2.249041015 1.5623943 -0.88685460 5.23952809
f_body_n4_t1 1.920366698 1.7185490 -0.44175115 6.77178836
f_body_n1_t2 1.602497719 1.5366898 -1.51546581 4.53498282
f_body_n2_t2 1.490072507 2.1342403 -2.07665200 5.71590120
f_body_n3_t2 3.004104444 1.2065585 0.88768375 5.40344697
f_body_n4_t2 1.006368311 1.1745803 -1.92303071 3.13812846
f_body_n1_t3 2.935316123 1.5960795 0.12325283 6.08925689
f_body_n2_t3 3.207557836 1.4854862 0.65255579 6.09422554
f_body_n3_t3 1.825944726 1.4909447 -0.75891584 5.28242235
f_body_n4_t3 4.493960648 2.0041912 1.13737928 8.48069745
sigma_e 1.163451959 1.3358080 0.05728349 4.57031299
sigma_hr 1.097063883 1.0925747 0.04970679 3.92241782
sigma_body 1.105672017 0.4808567 0.33428851 2.18500816
c_e_1 -2.590221808 0.9329944 -4.56878218 -0.89342448
c_e_2 -0.736207512 0.5757990 -1.93693628 0.30522640
c_e_3 0.233032975 0.5397744 -0.81525620 1.34900019
c_hr_1 0.559681445 0.6149305 -0.56188720 1.85774239
c_hr_2 2.351823838 1.0160333 0.64800451 4.61528595
c_body_1 -3.319405639 2.2604276 -9.87711969 -0.74233403
c_body_2 -1.774732378 1.0330290 -4.45708587 -0.34920066
c_body_3 -1.147109826 0.6979263 -2.84662593 -0.07616401
c_body_4 -0.782233757 0.5568800 -2.01395094 0.16376616
c_body_5 -0.546834694 0.4918412 -1.57788496 0.37887587
c_body_6 -0.269197560 0.4437017 -1.16517754 0.58749350
c_body_7 -0.016720799 0.4336173 -0.80421459 0.81480846
c_body_8 0.818056044 0.5055725 -0.07427233 1.92403837
c_body_9 1.453760796 0.6194285 0.37590679 2.80068910
c_body_10 1.763228035 0.6651253 0.57346015 3.28159577
c_body_11 3.661580482 1.0968645 1.83732355 5.86169802
lp__ -89.908971873 14.7324540 -117.60427823 -63.69930211
@martinmodrak, @asael_am
Now, my questions are
1- whether the beta for the ordinal predictor is for the entire time-series or should be time-specific?
2-Can we say the beta coefficients are the maximum beta for the ordinal variable or I should include mo function? What about I add the monotonic effect to this syntax? @paul.buerkner
3- There are 1200 people with three time-points, but if I want to add the fourth time-point, 400 people do not have the data as they either died or have not been observed yet. Can this be modeled as truncation because the whole data does not exist for 400 people?
Thanks a lot for your inputs.
I don’t know how to answer any of the 3 questions! I will see your code and try to answer 1. And also tag @martinmodrak so he got a notification.
In question 2? Are you talking about the estimator? (Talking of MAP)
@asael_am, in the literature, it is written that despite forecast analyses., in time-series the beta coef indicates the beta for the entire time-series. I wonder if this is the case in Stan as well?
And If I should create a matrix for coefficients to get time-specific coef?
Thanks
@asael_am: yes, I meant the estimator.
Well If you are doing a MCMC method (default in Stan) the estimator is the posterior mean,
So I don’t think you got a maximum
That is true regarding the mean posterior, in general. But what I meant was different. I refer you to the https://onlinelibrary.wiley.com/doi/abs/10.1111/bmsp.12195 for what I meant by maximum estimate.