HIERARCHICAL GAUSSIAN SPATIO TEMPORAL MODELS

[Jo Lweny]

Hell Stan Team,

Thanks for your previous direction, I got a lot of help with defining my covariance. The forum guys are taking a bit long.

Generally, we are trying to  model an underlying Gaussian Process with binomial sampling distribution(likelihood). Our models

are hierarchical over the parameters, Bk,which I give a joint normal prior and over hyper-parameters of the covariance matrix. We approximate the spatial dependency using the non-stationary,

anisotropic matern covariance ,whose parameters : length-scale(rho) and signal-variance(sigma-squared) are given the inverse_gamma(0.001,0.001) and uniform(0,200) priors

respectively.

Specifically my target model is: "glmm-model.png"- attached in this post.

#Now if fj is an Auto-regressive dynamics- for example, when each observation was made at a yearly time-point, which is very important for our inference- then I prefer to model it using an ar model or a random-walk ,where most appropriate.

 

#My Challenge

I have set up the model as given in the attached script, apparently it gives back results that are not consistent with the expected/known trend. I thought first to diagnose convergence but ended up with inconsistent volatility across the sampling chains. I suspected my initial parameterization values to be a possible cause and have since given different permutations of priors and inits, including specifying flat priors but the convergence is not clearly coming out.

 

I would like you to help me preview my model specification if it is indeed consistent with the literature of gaussian hierarchical models especially as specified in this mail, or even point to an unintentional failure that might have occurred. Your comments are highly regarded! I have attached an example dataset(Sl_ageCorrected*), stan model(sl_multigp.stan)  with prediction raster files: (sl_prec.asc,sl_urb.asc). You can also see the surface prediction I got initially(sl_PR_180615)

 

Kind Regards!

Lukio

 

[Daniel Lee]

Hi Lukio,

Sorry no one's replied. Could you take a few minutes to clean up your question so it's clear?

After reading your post again, it's still unclear what you're asking.



Daniel

On Monday, August 10, 2015 at 11:05:34 AM UTC-4, Jo Lweny wrote:

Hell Stan Team,

Thanks for your previous direction, I got a lot of help with defining my covariance. The forum guys are taking a bit long.

Generally, we are trying to  model an underlying Gaussian Process with binomial sampling distribution(likelihood). Our models

are hierarchical over the parameters, Bk,which I give a joint normal prior and over hyper-parameters of the covariance matrix. We approximate the spatial dependency using the non-stationary,

anisotropic matern covariance ,whose parameters : length-scale(rho) and signal-variance(sigma-squared) are given the inverse_gamma(0.001,0.001) and uniform(0,200) priors

respectively.

Specifically my target model is: "glmm-model.png"- attached in this post.

#Now if fj is an Auto-regressive dynamics- for example, when each observation was made at a yearly time-point, which is very important for our inference- then I prefer to model it using an ar model or a random-walk ,where most appropriate.

 

#My Challenge

I have set up the model as given in the

...

[Jo Lweny]

I am sorry if my question sounded vague. I was simply asking the forum guys to help me check what's wrong with my stan code, given the data I was using and the covariate grids- all attached.

When I call stan() with a list of input data - based on the ".stan" code attached here- as given below, I end up with the an error ,also posted at the end of this reply.

#sl_data is the imported csv that has Pf(The number of people tested positive for malaria parasite)

sl_data_training<- list(D=4,

N_obs=270,

y_obs=sl_data$Pf,

N_trials<- sl_data$Pf,

x_obs=array(data=c(Long=sl_data$Long,sl_data$Lat,sl_data$Urbanizati,sl_data$prec),dim=c(270,4),dimnames=list(NULL,c("Long","Lat","Urbanizati","prec"))),

N_new=133947,

x_new=array(data=c(xy_new[,1],xy_new[,2],urbs,precs),dim=c(133947,4),dimnames=list(NULL,c("Long","Lat","Urbanizati","prec")))

);

#make stanfit object

fit<- stan(model_code=sl_multi_gp,data=sl_data_training,iter=100,chains=2) ;

Error in stanc(file = file, model_code = model_code, model_name = model_name,  : 

  failed to parse Stan model 'c0963fce9d899b62b11e06534aaaa3b9' with error message:

SYNTAX ERROR, MESSAGE(S) FROM PARSER:

ERROR at line 12

 10:     vector[D] x_new[N_new];  //vector of arrays new x, raster grid values and their locations

 11:     int<lower=2> N;

 12:     N<- N_obs + N_new;

               ^

 13:     }

PARSER EXPECTED: <one of the following:

  a variable declaration, beginning with type,

      (int, real, vector, row_vector, matrix, unit_vector,

       simplex, ordered, positive_ordered,

       corr_matrix, cov_matrix,

       cholesky_corr, cholesky_cov

  or '}' to close variable declarations>

I am trying to model the distribution of malaria across a particular region by the covariates : precipitation(attached as sl_prec), urbanization(attached as  sl_urban.asc).

The Number of People Examined is the variable :"Ex" in the  csv Sl_AgeCorrectedEx10_300615, which also contains Number found malaria positive Pf.

My observation is therefore binomial. The data also contain the variables  prec,urbanization, Longitude ,Latitude as each observation was taken at a geo-referenced point.

Am I spcifying my model wrongly? 

[Jonah]

That particular syntax error is because in transformed data you need to declare everything before doing any assignments. For example, you have 

 ...
 int<lower=2> N;

 N<- N_obs + N_new;

 vector[D] x[N];
 ...

but you can't assign a value to N until after declaring x, etc. 

[Jonah Sol Gabry]

Also, that's true in the other blocks that allow assignment too. The exception is that you can make new declarations in local blocks. For example:

model {

 real a;

 real c;

 a <- 2;

 {

  real b;

  b <- 2;

  c <- a + b;

 }

}

is ok, but 

model {

 real a;

 real c;

 a <- 2;

 real b;

 b <- 2;

 c <- a + b;

}

will give you an error.