Re: [stan-users] Just first iteration with Hierarchical model responses and response times

[Mike Lawrence]

Several issues:

- You'll need to explain your data. I don't understand what you're doing with NA data. It looks like you're only ever asserting sampling when the data is not NA, in which case there's no need to include the NA cases in the first place. Including NA data is usually done when you want to use the non-NA data to infer what the NA data might have been. 

- I also suspect you might be able to get a speedup by using a non-centered parameterization on paramP and paramQ and thereby avoiding looping over them to assert the sampling of their rows individually and instead simply asserting to_vector(paramP_raw) ~ normal(0,1), etc.

- why are you using inv_sqrt(0.001) in the specification of the scale for the normal priors on muP & muQ? Why not just put 30? (Stan gurus: if 30ish is indeed desired as the scale of these priors, should we furthermore expect improved efficiency from a non-centered parameterization where, for example, muP_raw is defined as normal(0,1) and muP is defined in the transformed parameters section as muP_raw*30 ?

- You are computing diag_pre_multiply(sigmaP,L_OmegaP) for every iteration of the outter loop over 1:P, and, worse, computing diag_pre_multiply(sigmaQ,L_OmegaQ) on every inner loop over 1:Q. You should just compute both once, storing the result to a variable that is then referenced inside the loops. 

- You also might try to move all computations you can to the transformed parameters block. For example, here's your code in a cleaned-up format (please try to use consistent spacing) with those cholesky-decomposed covariance matrices computed in the transformed parameters block:

  data {

    int P ; // # Individus

    int Q ; // # Items

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

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

    int ItemsNA[P,Q] ; // Responses

  }

  parameters {

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

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

    vector[2] muP ; // mean of multivariate normal for parameters (theta, tau) of persons

    cholesky_factor_corr[2] L_OmegaP ; 

    real<lower=0> sigmaP2 ; 

    vector[2] muQ ; // mean of multivariate normal for parameters (betaR, betaT) of items

    cholesky_factor_corr[2] L_OmegaQ ; 

    vector<lower=0>[2] sigmaQ ;

    vector<lower=0>[Q] sigmaT ; // variance of Log times

  }

  transformed parameters{

    matrix LP ;

    matrix LQ ;

    {

      // I create vector sigmaP to fix sigmaP[1] = 1

      vector[2] sigmaP ;

      sigmaP[1] <- 1 ;

      sigmaP[2] <- sigmaP2 ;

      LP <- diag_pre_multiply(sigmaP,L_OmegaP)

      LQ <- diag_pre_multiply(sigmaQ,L_OmegaQ)

    }

  }

  model {

    // Priors

    muP ~ normal(0,inv_sqrt(0.001)) ;

    sigmaP2 ~ cauchy(0,2.5) ;

    L_OmegaP ~ lkj_corr_cholesky(2) ;

    muQ ~ normal(0,inv_sqrt(0.001)) ;

    L_OmegaQ ~ lkj_corr_cholesky(2) ;

    sigmaQ ~ cauchy(0,2.5) ;

    sigmaT ~ cauchy(0,2.5) ;

    for(p in 1:P){

      paramP[p] ~ multi_normal_cholesky(muP,LP) ;

      for(q in 1:Q){

        paramQ[q] ~ multi_normal_cholesky(muP,LQ) ;

        if(IndNA[p,q]==0){

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

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

        }

      }

    }

  }

[Bob Carpenter]

> On Jun 7, 2016, at 11:34 AM, Mike Lawrence <mike....@gmail.com> wrote:
>
> Several issues:
>
> - You'll need to explain your data. I don't understand what you're doing with NA data. It looks like you're only ever asserting sampling when the data is not NA, in which case there's no need to include the NA cases in the first place. Including NA data is usually done when you want to use the non-NA data to infer what the NA data might have been.

I also couldn't follow this.

> - I also suspect you might be able to get a speedup by using a non-centered parameterization on paramP and paramQ and thereby avoiding looping over them to assert the sampling of their rows individually and instead simply asserting to_vector(paramP_raw) ~ normal(0,1), etc.

If there's not a lot of data, non-centering is a win.

The "centering" in question is of both the mean and the log scale
around zero.

>
> - why are you using inv_sqrt(0.001) in the specification of the scale for the normal priors on muP & muQ? Why not just put 30? (Stan gurus: if 30ish is indeed desired as the scale of these priors, should we furthermore expect improved efficiency from a non-centered parameterization where, for example, muP_raw is defined as normal(0,1) and muP is defined in the transformed parameters section as muP_raw*30 ?

Yes, because then muP_raw, which is what's actually being sampled, is
on the unit scale, so we don't have to estimate that scale during
warmup. So warmup should go faster and be more accurate.

> - You are computing diag_pre_multiply(sigmaP,L_OmegaP) for every iteration of the outter loop over 1:P, and, worse, computing diag_pre_multiply(sigmaQ,L_OmegaQ) on every inner loop over 1:Q. You should just compute both once, storing the result to a variable that is then referenced inside the loops.

What you really need to do is vectorize these. That way, even the
Cholesky factor is only factored once.

> - You also might try to move all computations you can to the transformed parameters block. For example, here's your code in a cleaned-up format (please try to use consistent spacing)

+1 on consistent spacing (and no tabs!).

> with those cholesky-decomposed covariance matrices computed in the transformed parameters block:

Computing in the transformed parameters block saves the parameters
in the output, so only define there if you want to save. Otherwise,
you can use a local variable in the model block. Either way, it's
better to just vectorize. Then you can just write the expression directly
inside the call to the density.

- Bob

> --
> 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,

Thanks for you help.

When submitting a computerized test limited in time to the subjects, we can record two types of informations, responses and responses times that suejts put to answer the various questions of the test.

Items: matrix (P × Q)

Time: matrix (P × Q)

Furthermore, it is known that parameters such topics as competence (theta) which can be estimated from the Items and speed (tau) that can be estimated from the responses times.

It is the same for the items where we can estimate the difficulty of each item with items (betaT) and the time required to answer an item (betaR).

We know that the items (0 and 1) follow a parameter Bernoulli  P = 1 / (1 + exp (theta-betaR)).

It is known that the logarithmic transformation of the response times follow a normal distribution (tau-betaT, sigma)

But the parameters of the subjects (theta and tau) follow a multivariate normal (muP,sigmaP)

and item parameters (BetaR and BetaT) a multivariate normal. (muQ,OmegaQ)

muP = (mu_theta, mu_tau)

sigmaP =(sigma_theta, cov_thetatau, cov_thetatau, sigma_tau)

muQ = (mu_betaR, mu_betaT)

OmegaP =(1, covvar_thetatau, covvar_thetatau, sigma_betaT)

paramP = (theta,tau) , paramQ = (betaR, betaT).

In this data, there are missing data NA. As Stan does not support NA, I have identified the IndNA data.

So I made different cal estimates of observed data by ignoring NA missing data.

I hope to bring explicitations

Bests regards.

[Mike Lawrence]

Ah, so the inclusion of the NA indicator is to maintain a PxQ data matrix amidst missing data. Note that the alternative solution (that would get rid of the need for the if() check in the model) to this is to reformat the data to long format and toss the NA entries, after which you'd have a data section akin to:

  data {

     int P ; // # Individus

     int Q ; // # Items

     int N ; // total # of non-NA data points

     vector[N] LogTimes ; // Log responses times

     vector[N] Responses ; // Responses

     vector[N] P_label ; // P labels associated with each of the N data points

     vector[N] Q_label ; // P labels associated with each of the N data points

   }

You then need to re-code your model to account for the new data structure. A similar (though only including one random effect where you have two) hierarchical model with long data can be seen here: https://gist.github.com/mike-lawrence/4989099

If I were you, I'd try our earlier suggestions first, and if you're still sampling slowly, try reformatting to long and tossing the NA data. 

(Stan gurus: would you expect the conversion to long-format to avoid the if() check in the model to be worth it? My intuition on what is more efficient is still sketchy...)

--

[Eric Innocents EBOULET]

Hi Mike,

I try and I come back to you.

Best regards.

[Eric Innocents EBOULET]

Hi Mike,

I modified the model with your suggestion and it works. However the convergence of the estimations are not possible.

To do this, I generate 2 datas, responses and responses times given parameters. So, I tried to find them. But I can not.

If you put files in the same environment, you can use only StanExcution.R to create datas and to do estimations.

File StanExecution.R Generate the datas and excute the model to estimate the parameters

File S1M_TimeResponseDependent.stan contains the model

File FunctionItemsLogTimesData.R contains different functions used in file StanExecution.R 

File FunctionsThesisStart.R contains different functions used in file StanExecution.R 

data {

    int P; // # Individus

  int Q; // # Items

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

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

  int ItemsNA[P,Q]; // Responses

 

}

parameters {

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

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

  

 vector[2] muP;

 cholesky_factor_corr[2] L_OmegaP;

 real<lower=0> sigmaP2; 

 vector[2] muQ;

 cholesky_factor_corr[2] L_OmegaQ;

 vector<lower=0>[2] sigmaQ;

 

    vector<lower=0>[Q] sigmaT;

}

transformed parameters{

    matrix[2,2] LP;

    matrix[2,2] LQ;

  

    vector[2] sigmaP;

    sigmaP[1] <- 1;

    sigmaP[2] <- sigmaP2;

  

    LP <- diag_pre_multiply(sigmaP,L_OmegaP);

    LQ <- diag_pre_multiply(sigmaQ,L_OmegaQ);

}

model {

  

  // Priors

  muP ~ normal(0,1);

    sigmaP2 ~ invgamma(0.001,0.001);

  L_OmegaP ~ lkj_corr_cholesky(1);

   

    muQ ~ normal(0,1);

    L_OmegaQ ~ lkj_corr_cholesky(1);

    sigmaQ ~ invgamma(0.001,0.001);

    

    sigmaT ~ invgamma(0.001,0.001);

    

    

// 1st level

 for(p in 1:P)

  {

     for(q in 1:Q)

     {

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

       {

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

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

     }

   }

  }

 

  // 2nd level

  paramP ~ multi_normal_cholesky(muP,LP);

  paramQ ~ multi_normal_cholesky(muQ,LQ);

}

Bests regards

[Mike Lawrence]

First, your stan code doesn't compile because invgamma doesn't exist. It needs to be inv_gamma.

Second, I'm not clear on what you mean by "convergence of the estimations are not possible". With a iter=1000 run (which takes about 7mins on my machine), most parameters have a good ESS of about 500, with muP and muQ being the exceptions, so it looks like you'd need more iterations for those. Or are you saying that the posterior on one or more of the parameters doesn't mesh with the known true values (since you simulated the data)? Again, only muP and muQ look off from their true values in my iter=1000 run, which isn't surprising given the low ESS. Do you have a high iteration run that achieves good ESS/rhat for muP & muQ and still observe inconsistency between the posterior and known true values?

[Eric Innocents EBOULET]

Hi Mike,

Thank you a lot.

I try with a high iteration run - 20 000. But the graphicals diagnostics of convergences of parameters  with fonctions stan_rhat and

stan_ess on Rhat and n_eff are below :

To interpretation with the graphicals are most parameters have not a good ESS and rhat. 

If I finish this convergence, I'll concentrate me if estimations of parameters are near the known true values (since you simulated the data).

What am I doing change in my models.

Best regards.