Hierarchical modeling joined response and response times.

[Eric Innocents EBOULET]

Hi Ratan's users,

I'm writing to make sure that my code I write reflects the hierarchical model defined below.

1. To first level, Responses are defined with ranch model and log response times to normal distribution
2. To second level, parameters of persons are multivariate normal distribution. Also parameters of items.

If applicable, what is the best way to write the priors for variances covariances matrix and mu because I doubt of this in my code.

data {

  int P; // # Individus

  int Q; // # Items

  real logTwNA[P,Q]; // Log responses times

  real IndwNA[P,Q]; // NA index  (1 if NA, 0 ifelse)

  int RwNA[P,Q]; // Responses

}

parameters {

vector[2] paramP[P]; // Parameters of persons (theta,tau)

vector[2] paramQ[Q]; // Parameters of items (betaR,betaT)

        // Parameters's distributions of parameters of persons

  vector[2] muP;

  corr_matrix[2] rhoP;

  vector<lower=0,upper=100>[2] sigmaP;

  

        // Parameters's distributions of parameters of items

  vector[2] muQ;

  corr_matrix[2] rhoQ;

  vector<lower=0,upper=100>[2] sigmaQ;

  

        // Sd of log response time of each item 

  vector<lower=0,upper=100>[Q] sigmaT;

}

transformed parameters {

cov_matrix[2] omegaP;

cov_matrix[2] omegaQ;

        // Covariance matrix defined for parameters of persons

for(i in 1:2)

{

omegaP[i,i] <- sigmaP[i]*sigmaP[i]*rhoP[i,i];

}

omegaP[1,2] <- sigmaP[1]*sigmaP[2]*rhoP[1,2];

omegaP[2,1] <- omegaP[1,2];

        // Covariance matrix defined for parameters of items

for(i in 1:2)

{

omegaQ[i,i] <- sigmaQ[i]*sigmaQ[i]*rhoQ[i,i];

}

omegaQ[1,2] <- sigmaQ[1]*sigmaQ[2]*rhoQ[1,2];

omegaQ[2,1] <- omegaQ[1,2];

}

model {

  

  // Priors omegaP ~ inversewishart(2,matrix(c(1,0,0,1,ncol=2)), omegaQ ~ inversewishart(2,matrix(c(10,1,1,10,ncol=2))

  // muP ~ MVN(c(0,0),omegaP), muQ ~ MVN(c(0,0),omegaQ)

  matrix[2,2] omegaQ0;

  vector[2] muP0;

  vector[2] muQ0;

  

  muP0 <- rep_vector(0,2);

  mu00 <- rep_vector(0,2);

  omegaQ0 <- diag_matrix(rep_vector(10,2));

  omegaQ0[1,2] <- 1;

  omegaQ0[2,1] <- 1;

  

  omegaP ~ inv_wishart(2,diag_matrix(rep_vector(1,2)));

  omegaQ ~ inv_wishart(2,omegaQ0);

  

  muP ~ multi_normal(muP0,omegaP);

  muQ ~ multi_normal(muQ0,omegaQ);

    

  // 1st level

  for(p in 1:P)

  {

    for(q in 1:Q)

    {

      if(IndwNA[p,q]==0)

      {

        RwNA[p,q] ~ bernoulli_logit(paramP[p,1]-paramQ[q,1]);

        logTwNA[p,q] ~ normal(paramQ[q,2] - paramP[p,2],sigmaT[q]);

      }

    }

  }

  

  // 2nd level

  for(p in 1:P)

  {

    paramP[p] ~ multi_normal(muP,omegaP);

  }

  

  for(q in 1:Q)

  {

    paramQ[q] ~ multi_normal(muQ,omegaQ);

  }

}

Thank you to understand me, because I try to have even more insurance in the writing Bayesian method with Rstan.

Bests regards.

Eric I. E.

[Bob Carpenter]

Check out the section in the regression chapter of the manual
on multivariate priors to see how to do this robustly and
efficiently in terms of Cholesky factors of correlation matrices
plus scales.

You may want to read our other prior advice, because we strongly
recommend against uniform priors on an interval.

- Bob

P.S. I think you have the spelling corrector turned up too high.

> --
> You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
> To post to this group, send email to stan-...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

[Eric Innocents EBOULET]

Hi Bob,

Thanks a lot.

I look.

Best.

[Eric Innocents EBOULET]

Hi Rstan's

With the suggestion of Bob, I don't already arrived to run my model.

I write hierarchical multimodel where in first level, I write Rasch model for responses and 

lognormal for responses times.

In second the parameters of persons are bivariate normale. Also parameters of Items.

With bivariate priors in manuel of rstan, I write the script below :

data {

  int P; // # Individus

  int Q; // # Items

  real logTwNA[P,Q]; // Log responses times

  real IndwNA[P,Q]; // NA index  (1 if NA, 0 ifelse)

  int RwNA[P,Q]; // Responses

}

parameters {

matrix[2,P] zP;

  vector[2] muP;

  cholesky_factor_corr[2] L_rhoP;

  vector<lower=0>[2] sigmaP;

  

  matrix[2,Q] zQ;

  vector[2] muQ;

  cholesky_factor_corr[2] L_rhoQ;

  vector<lower=0>[2] sigmaQ;

  

  vector<lower=0>[Q] sigmaT;

}

transformed parameters{

        // to define second level -  constraints of parameters (In looking in manuel)

matrix[2,P] paramP; // Parameters of persons (theta,tau)

matrix[2,Q] paramQ; // Parameters of items (betaR,betaT)

paramP <- to_matrix(muP) + (diag_pre_multiply(sigmaP,L_rhoP) * zP);

paramQ <- to_matrix(muQ) + (diag_pre_multiply(sigmaQ,L_rhoQ) * zQ) ;

}

model {

  // Priors

  sigmaT ~ cauchy(0,2.5);

  

  muP ~ normal(0,1);

  sigmaP ~ cauchy(0,2.5);

  L_rhoP ~ lkj_corr_cholesky(2);

  to_vector(zP) ~ normal(0,1); 

  

  muQ ~ normal(0,1);

  sigmaQ ~ cauchy(0,2.5);

  L_rhoQ ~ lkj_corr_cholesky(2);

        to_vector(zQ) ~ normal(0,1);

    

  for(p in 1:P)

  {

    for(q in 1:Q)

    {

      if(IndwNA[p,q]==0)

      {

        RwNA[p,q] ~ bernoulli_logit(paramP[1,p]-paramQ[1,q]);

        logTwNA[p,q] ~ normal(paramQ[2,q] - paramP[2,p],sigmaT[q]);

     }

    }

}

  

  // 2nd level

  //paramP ~ multi_normal(muP,(diag_pre_multiply(sigmaP,L_rhoP)*zP));

//paramQ ~ multi_normal(muQ,(diag_pre_multiply(sigmaQ,L_rhoQ)*zQ));

}

So, there is an error when I run this code.

Best regards

[Eric Innocents EBOULET]

Hi Rstan's

With the suggestion of Bob, I don't already arrived to run my model.

I write hierarchical multimodel where in first level, I write Rasch model for responses and 

lognormal for responses times.

In second the parameters of persons are bivariate normale. Also parameters of Items.

With bivariate priors in manuel of rstan, I write the script below :

data {

   int P; // # Individus

   int Q; // # Items

   real logTwNA[P,Q]; // Log responses times

   real IndwNA[P,Q]; // NA index  (1 if NA, 0 ifelse)

   int RwNA[P,Q]; // Responses

}

parameters {

matrix[2,P] zP;

    

   for(p in 1:P)

   {

     for(q in 1:Q)

     {

       if(IndwNA[p,q]==0)

       {

         RwNA[p,q] ~ bernoulli_logit(paramP[1,p]-paramQ[1,q]);

         logTwNA[p,q] ~ normal(paramQ[2,q] - paramP[2,p],sigmaT[q]);

      }

      }

}

  

   // 2nd level

   //paramP ~ multi_normal(muP,(diag_pre_multiply(sigmaP,L_rhoP)*zP));

//paramQ ~ multi_normal(muQ,(diag_pre_multiply(sigmaQ,L_rhoQ)*zQ));

}

So, there is an error when I run this code.

Best regards

Le mardi 23 février 2016 22:08:59 UTC+1, Bob Carpenter a écrit :

[Bob Carpenter]

The error message says that on line 27
of your program, you are operating on matrices or
vectors with sizes that don't match.

- Bob

[Eric Innocents EBOULET]

Thanks Bob,

I understood this, However, I don't know how to do this operation in figure R in Rstan.

In R, it is possible to do,

a = c(1,2)

b = matrix(1:6,ncol=3)

a+b.

You can see in file joined.

In rstan How to write this ?

Best for help.

[Bob Carpenter]

Would you mind just cutting-and-pasting rather than
using screen shots? I can't cut-and-paste from your
screen shots and I have to blow them up to read in my
mailer.

R's just filling column major:

> c(10, 100, 1000) + matrix(1:6, 2, 3)
[,1] [,2] [,3]
[1,] 11 1003 105
[2,] 102 14 1006

So if you literally want to do what R's doing in
terms of broadcasting (filling):

matrix add_broadcast(vector x, matrix y) {
int pos;
matrix xpy[rows(y), cols(y)];
if (num_elements(x) == 0) return y;
pos <- 1;
for (n in 1:cols(y)) {
for (m in 1:rows(y)) {
if (pos > num_elements(x)) pos <- 1;
xpy[m, n] <- xpy[m, n] + x[pos];
pos <- pos + 1;
}
}
return xpy;
}

But you probably want the simpler thing of taking a column vector
and adding it to each column, which you can do with a much simpler
loop.

- Bob

> On Mar 1, 2016, at 11:55 AM, Eric Innocents EBOULET <ericinnoce...@gmail.com> wrote:
>
> Thanks Bob,
>
> I understood this, However, I don't know how to do this operation in figure R in Rstan.
> In R, it is possible to do,
> a = c(1,2)
> b = matrix(1:6,ncol=3)
> a+b.
>
> You can see in file joined.
>
>
>

[Eric Innocents EBOULET]

Hi Bob, 

Thanks for your help.

Best

[Eric Innocents EBOULET]

Hi Bob,

When I looked Multivariate Priors for hierarchical Models and your suggestions with factor cholesky in the manuel (page 75),

transformed parameters {

matrix[J,K] beta;

beta <- u * gamma + (diag_pre_multiply(tau,L_Omega) * z)';

}

u is a matrix JxL

gamma is a matrix LxK

(diag_pre_multiply(tau,L_Omega) * z)' is a JxK because of transpose '

You explained it means beta ~ MultiNormal(ul*gamm, Sigma)

I want to write the same.

Theta is matrix 2xN ~ MultiNormal (u,Sigma) (A)

Sigma is matrix 2*2

u is vector with length equal 2

So Theta <- u + (diag_pre_multiply(tau,L_Omega) * z); where 

vector[2] tau;

cholesky_factor_corr[2] L_Omega;

matrix[2,N] z

My problems are :

1. Is it correct to define z is matrix 2xN to obtain Theta. Ifelse how to define z to obtain matrix 2xN of Theta (A)
   
2. If 1 is True, I don't calculate Theta because u is vector of length is 2 and drag_pre_multiply(tau, L_omega)*z is matrix 2xN, How to do this ?

Bests regards.

[Bob Carpenter]

> On Mar 1, 2016, at 8:09 PM, Eric Innocents EBOULET <ericinnoce...@gmail.com> wrote:
>
> Hi Bob,

I'm not the only one here!

> When I looked Multivariate Priors for hierarchical Models and your suggestions with factor cholesky in the manuel (page 75),
>
> transformed parameters {
>
> matrix[J,K] beta;
>
> beta <- u * gamma + (diag_pre_multiply(tau,L_Omega) * z)';
>
> }
>
>
>
> u is a matrix JxL
>
> gamma is a matrix LxK
>
> (diag_pre_multiply(tau,L_Omega) * z)' is a JxK because of transpose '
>
>
>
> You explained it means beta ~ MultiNormal(ul*gamm, Sigma)

That wasn't me --- Ben did that example in the manual.

Otherwise, I couldn't quite follow the question here. It's
just ordinary linear algebra rules for matrices and all the
signatures for functions are in the back of the manual. We
don't really have time to debug users' programs step by step.

- Bob

[Eric Innocents EBOULET]

Thank you for your availability. So, I am the fact of addressing you directly offended 

you .

Bests

[Jonah Gabry]

Hi Eric, 

Can you explain your question about the multivariate priors section again? I didn't quite understand. 

Best,

Jonah

[Eric Innocents EBOULET]

Hi Jonah,

Thanks !!

I want to write hierarchical model including 2 datas, responses (R) and responses times (RT). R and RT are matrix NxQ

In 1st level, each data is modeled with a different model independent

For Responses : the model is Rasch model with parameters theta and betaR

logit(P) = theta - betaR

For Responses Times : the model is lognormal with parameter tau, betaT and sigmaT.

log(RT) ~ Normal(tau - betaT, sigmaT)

In 2nd level,

paramP is a vector (theta and tau) and paramI is vector (betaR and betaT).

So paramP ~ MultiNormal (muP,sigmaP) and paramI ~ MultiNormal(muI, sigmaI)

One suggestion I got was to rely on what has been done one page 75 of the manuel to write 2nd level.

Using script with my understanding, I write the model below

data {

   int P; // # Individus

   int Q; // # Items

   real logTwNA[P,Q]; // Log responses times

   real IndwNA[P,Q]; // NA index  (1 if NA, 0 ifelse)

   int RwNA[P,Q]; // Responses

}

parameters {

matrix[2,P] zP; // using example in page 75

   vector[2] muP; // parameters of persons (theta and tau)

   cholesky_factor_corr[2] L_rhoP; (Inferior matrix)

   vector<lower=0>[2] sigmaP; 

  

   matrix[2,Q] zQ;

   vector[2] muQ;

   cholesky_factor_corr[2] L_rhoQ;

   vector<lower=0>[2] sigmaQ;

  

   vector<lower=0>[Q] sigmaT;

}

// I'm using transformed parameters because of explication in page 75 to obtain parmaP ~ MultiNormal(muP,sigmaP)

and paramI ~ MultiNormal(muI,sigmaI).

transformed parameters{

        // to define second level -  constraints of parameters (In looking in manuel)

matrix[2,P] paramP; // Parameters of persons (theta,tau)

matrix[2,Q] paramQ; // Parameters of items (betaR,betaT)

paramP <- to_matrix(muP) + (diag_pre_multiply(sigmaP,L_rhoP) * zP);

paramQ <- to_matrix(muQ) + (diag_pre_multiply(sigmaQ,L_rhoQ) * zQ) ;

}

model {

   // Priors

   sigmaT ~ cauchy(0,2.5);

  

   muP ~ normal(0,1);

   sigmaP ~ cauchy(0,2.5);

   L_rhoP ~ lkj_corr_cholesky(2);

   to_vector(zP) ~ normal(0,1); 

  

   muQ ~ normal(0,1);

   sigmaQ ~ cauchy(0,2.5);

   L_rhoQ ~ lkj_corr_cholesky(2);

        to_vector(zQ) ~ normal(0,1);

    

   for(p in 1:P)

   {

     for(q in 1:Q)

     {

       if(IndwNA[p,q]==0)

       {

         RwNA[p,q] ~ bernoulli_logit(paramP[1,p]-paramQ[1,q]);

         logTwNA[p,q] ~ normal(paramQ[2,q] - paramP[2,p],sigmaT[q]);

      }

}

        }

}

So my first question :

1. Is it correct my script to define the hierarchical model
2. If it is tree, when I run model, an error is displayed that it is not running because the definition of parameterP is not correct. 

In the parameterP, I want to do  in Rstan the same :

In R, it is possible to do :

a = 1:2

b = matrix(1:6, ncol=3)

a+b is possible to appear the result. In Rstan, it is not possible because, a is vector and b is matrix.

So how to define a and b to do a+b in Rstan.

Otherwise, if it's possible to give an idea to reorganize my script, It'll be very useful.

Bests regards.

[Bob Carpenter]

> On Mar 2, 2016, at 10:00 AM, Eric Innocents EBOULET <ericinnoce...@gmail.com> wrote:
>
> Thank you for your availability. So, I am the fact of addressing you directly offended
> you .

I wasn't offended --- I just want other people to answer questions.

- Bob

[Jonah Sol Gabry]

In R, it is possible to do :
a = 1:2
b = matrix(1:6, ncol=3)
a+b is possible to appear the result. In Rstan, it is not possible because, a is vector and b is matrix.
So how to define a and b to do a+b in Rstan.

Stan is stricter than R when it comes to arithmetic with matrices and vectors. R will cheat for you so that a + b is valid when a is a vector and b is a matrix, but Stan requires you to do proper matrix arithmetic. If you want to add a to each column of b you can use rep_matrix(a, 3) + b, which will make a 2x3 matrix where each column is a copy of a and then it will add that matrix to b.  

[Eric Innocents EBOULET]

Hi Bob,

Thanks you for your understanding.

Best.

[Eric Innocents EBOULET]

Hi Jonah,

Thanks a lot.

I appreciated.

Bests regards,

Eric

[Eric Innocents EBOULET]

Hi Rstan's Users,

Could you help me dot debug my problem please ? I don't know what this error in file joined means 

Best regards.

[Jonah Gabry]

Hi Eric, 

I've never seen this type of message before. Does this happen when running the same model you were asking about in this thread or is it a different model? What happens after you the messages stop?

Jonah

[Eric Innocents EBOULET]

Thanks Jonah.

It is the same model. Now I think the models runs. Mostly with your suggestion with rep_matrix.

Nothing is happening after this error. I stopped and I restart. The error still occurs.

I joined the model. I try also in internet if I find something for this.

Bests.

[Eric Innocents EBOULET]

It's apple's problem. I'trying to understand and I'll share you the information.

Bests Jonah.

[Eric Innocents EBOULET]

Hi Ratan's Users ?

Is it normal to my model run after 13h. (Files joined)

Also, all parameters are estimated very good except the matrix L_rhoP defined cholesky_factor_corr[2] L_rhoP

[Eric Innocents EBOULET]

Hi Ratan's Users ?

Is it normal to my model run after 13h. (Files joined)

Also, all parameters are estimated very good except the matrix L_rhoP defined cholesky_factor_corr[2] L_rhoP (Files joined)

Bests

Le samedi 5 mars 2016 00:48:27 UTC+1, Eric Innocents EBOULET a écrit :

Hi Ratan's Users ?

Is it normal to my model run after 13h. (Files joined)

Also, all parameters are estimated very good except the matrix L_rhoP defined cholesky_factor_corr[2] L_rhoP

[Jonah Sol Gabry]

Hi Eric,

The NaN for the rhat for those elements of the cholesky are ok because a cholesky factor is a triangular matrix. Those elements aren't estimated. 

Regarding the time it takes to fit the model there's no "normal". It will depend on how many parameters there are, the posterior geometry, the amount of data, how efficient your code is, etc. 

Jonah

[Eric Innocents EBOULET]

Thanks Jonah !