Question on multilevel regression with group-level interaction

[David Pedrosa]

Hi all,

I recently started to use Stan for differnt models I have been trying to fit. Basically, I am interested in building a regression model from the scratch adding all kinds of additional information/complexity. The manual, this list and all that can be found especially on github have been incredibly helpful to get things running and get first results. So right now I am a little stuck at fitting a multilevel regression model. It's repeated mesaures data recorded from a a group of 12 subjects and the outcome is the log-transformed amplitude of a signal. This signal (muscle recordings of the arm) was recorded while we put an electrode at very specific parts of the brain, so that the only predictor is the distance towards this target in mm. Everything works fine when I code model with complete pooling or another one with partial pooling, that is the one with the slope/intercept varying per subject. The problem arises when I want to incorporate a group-level predictor (in this case a single number between 1-5 which describes the amount of electrodes used at the brain for the specific subject "s"). I basically copied (and modified) the code from the pystan example with the radon contamination (originally from Gelman and Hill 2006). This is the code:

data{
int<lower=0> N; //number of total observations
int<lower=1> J; //number of subjects
int<lower=1,upper=J> s[N]; // suject ID
int<lower=1> K; // sum intercept + slope(s)
vector[J] e; // group-level interactions, number of electrodes

vector[N] y; // outcome variable
vector[N] x; // predictor(s)
}

transformed data{
// Data are z-transformed
vector[N] x_std;
vector[N] y_std;
x_std = (x - mean(x)) / sd(x);
y_std = (y - mean(y)) / sd(y);
}

parameters{
real betamu_std[K]; // group averages for intercept[1] and slope[2]
real<lower=0> sigma_y_std; // total error
vector[K] e_betaj[J];
vector[K] sigma_j;
vector[K] gamma; // gamma is the random part of the group-level predictor for the slope
}

transformed parameters{
vector[N] yhat; //
vector[J] beta_std[K]; // number of beta vectors with group-level interaction
vector[J] beta_raw[K]; // number of beta vectors without group-level interaction
vector[J] beta[K]; // number of (back-transformed) beta vectors

real betamu[K]; // group averages for intercept[1] and slope[2]
real<lower=0> sigma_y; // (back)-transformed error

for(j in 1:J) {
for(k in 1:K) {
beta_raw[k,j] = betamu_std[k] + sigma_b[k] *e_betaj[j,k];
beta_std[k,j] = beta_raw[k,j] + e[s[n]] * gamma[k];
}
}

for(n in 1:N) {
yhat[n] = beta_std[1,s[n]] + beta_std[2,s[n]] * x_std[n];
}

// Transform parameters back to original scale
for(j in 1:J) {
beta[1,j] = sd(y) * (beta_std[1,j] - beta_std[2,j] * mean(x) / sd(x)) + mean(y);
beta[2,j] = beta_std[2,j] * sd(y) / sd(x);
y_sig = sd(y) * sigma_y;
}

betamu[2] = betamu_std[2] * sd(y) / sd(x);
betamu[1] = betamu_std[1] * sd(y) + mean(y) - betamu_std[2] * mean(x) * sd(y) / sd(x);
}

model{
//Priors
sigma_y ~ cauchy(0, 2.5);
sigma_b ~ cauchy(0, 2.5);
betamu_std ~ normal(0,5);
gamma ~ normal(0,1);

for(k in 1:K) {
for(j in 1:J) {
e_betaj[j,k] ~ normal(0,1);
}
}

y_std ~ student_t(4, yhat, sigma_y);
}

"

The problem is that the group-level predictor "flips" the direction of the slope from negative to positive and the checks I did to assess the model fit look really bad. I know this may be just because these group-level predictors do not add anything to the model fit but because I don't trust (yet) my model definition I wanted to know if the way I am coding especially for beta_raw and beta_std make sense at all. I am grateful for every comment on the code. And thanks for the great work you guys do.
Best


David

[Bob Carpenter]

I don't quite follow what you're doing with the code for
multilevel regression (we'd usually put a multivariate prior
on the two variables per group---is that what your transforms are
doing?) , but there are lots of places to improve
its efficiency.

For a start, take all those expensive computations on the data
like mean(x) and sd(y) and even sd(y) / sd(x)
and do them in transformed data and store the result.

You can vectorize almost every statement in the program after
that, which will also give you speedups.

for (n in 1:N)

yhat[n] = beta_std[1,s[n]] + beta_std[2,s[n]] * x_std[n];

as

yhat = beta_std[1, s] + beta_std[2, s] .* x_std;

Won't help with efficiency, but hopefully easier to read and digest.

And you can vectorize this:

for (k in 1:K)
for (j in 1:J) {
e_betaj[j,k] ~ normal(0,1);

by declaring e_betaj as a matrix as:

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

This

for(j in 1:J) {
for(k in 1:K) {
beta_raw[k,j] = betamu_std[k] + sigma_b[k] *e_betaj[j,k];
beta_std[k,j] = beta_raw[k,j] + e[s[n]] * gamma[k];
}
}

can be coded as

for (j in 1:J) {
beta_raw[ , j] = betamu_sd + sigma_b .* e_betaj[j, ];
beta_std[ , j] = ...

Is K = 2? I find it hard to follow when a bound's used as a declared
data variable and hard coded. It means the program won't even work for
K != 2 and you shouldn't have K as data.

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

[David Pedrosa]

Hi Bob,

first of all, please apologise the somehow very bulky piece of code I provided in a first instance and thank you very much for your suggestions with the vectorization. I have tried to tidy up the unnecessary parts (transformation of the data and back transformation to obtain z-scores) and to apply the steps of vectorisation you proposed (although for clarity I have transformed beta to beta0 and beta1, I will probably vectorise this when I'm sure the model is right) . This is the code now:

data{
    int<lower=0> N;                     //number of total observations
    int<lower=1> J;                      //number of subjects
    int<lower=1,upper=J> s[N];     // suject ID

    vector[J] e;                            // group-level interactions (number of electrodes)

    vector[N] y;                           // outcome variable
    vector[N] x;                           // predictor(s)   
}

parameters{
    real beta_mu[2];                                    // group averages for intercept[1] and slope[2] 
    real<lower=0> sigma_y;                         // total error
    matrix[2,J] e_betaj;
    vector<lower=0,upper=100>[2] sigma_b;   
    vector[2] gamma;                                  // random part of group-level predictors (electrode count) (?)
}
  
transformed parameters{
    vector[N] y_hat;
    vector[J] beta0_raw;
    vector[J] beta1_raw;
    vector[J] beta0;
    vector[J] beta1;

    beta0_raw = beta_mu[1] + sigma_b[1] * e_betaj[1]';
    beta1_raw = 0.1 * beta_mu[2] + sigma_b[2] * e_betaj[2]';         // why multiplying with 0.1 ?!?

    beta0 = beta0_raw + e * gamma[1];
    beta1 = beta1_raw + e * gamma[2];
   
    for (i in 1:N)
        y_hat[i] = beta0[s[i]] + beta1[s[i]] * x[i];

}
  
model {
    //Priors
    sigma_y ~ cauchy(0, 2.5);
    sigma_b ~ cauchy(0, 2.5);

    beta_mu ~ normal(0, 5);
    gamma ~ normal(0,1);
    to_vector(e_betaj) ~ normal(0, 1);
   
    y ~ student_t(4, y_hat, sigma_y);
}

Unsurprisingly, it runs much faster than before, but the results remain the same. Without the group-level predictor ("e") the slope is negative while after including it, the slope is reversed and the fit of the model is a disaster. I haven't found too much on varying intercepts and slopes with group-level predictors but my hunch would be that it is not valid to include the "e * gamma"-term when there is one value of "e" per subject and it ranges from 1-5 the way I'm doing it. And it may not make much sense but I thought I run the model like this and when this works I extend it incorporating multivariate priors according to chapter 6.12 of the rstan manual as you also suggested (Although if the way I'm doing it is nonsense and it makes more sense to do it in one step, it would be of course helpful to know as well). I know these are rather basic questions and I appreciate the time you are spending. Thanks in advance.

Best
David 

[Bob Carpenter]

At least that means you refactored it properly :-)

> Without the group-level predictor ("e") the slope is negative while after including it, the slope is reversed and the fit of the model is a disaster.

> I haven't found too much on varying intercepts and slopes with group-level predictors but my hunch would be that it is not valid to include the "e * gamma"-term when there is one value of "e" per subject and it ranges from 1-5 the way I'm doing it.

By ranges from 1 to 5 do you mean that there are five unrelated values? If
so, you probably don't want to multiply but instead index, so that e is declared
as an integer and you use:

gamma[e[j]]

where one of the gamma is set to 0, say gamma[1] = 0, for identifiability.

If it's ordinal, you can sometimes get away with treating it as a scalar.

> And it may not make much sense but I thought I run the model like this and when this works I extend it incorporating multivariate priors according to chapter 6.12 of the rstan manual as you also suggested (Although if the way I'm doing it is nonsense and it makes more sense to do it in one step, it would be of course helpful to know as well). I know these are rather basic questions and I appreciate the time you are spending. Thanks in advance.

You're debugging this the right way starting from simple things that work and
then adding bits until something breaks so you know where the problem is (or at
least that it's a combination of the new bit and whatever was there before).

- Bob

[David Pedrosa]

Yes e is just the number of electrodes used for this recording at patient j and so unrelated between the subjects. I will try your suggestion with the index but just to get it clear, you mean extending the gamma to a vector of five (with the first one being 0) and replacing:

beta0 = beta0_raw + e * gamma[1] with

with

beta0 = beta0_raw * gamma[e,1]

David

[Bob Carpenter]

No, I meant just having beta0_raw + gamma[e[j]] where e in 1:5 is
the index into the gamma coefficients.

- Bob

[Daniel Myall]

I've fit similar sort of models on similar types of data.

Do you have plots of the raw data (i.e., distance towards target versus signal amplitude with a linear fit for each subject, with separate facets for each count of electrodes? If not, ggplot2 is great for this). This will help considerably in deciding how the number of electrodes should be incorporated into the model. With 12 subjects you will likely have difficulty in estimating the effect of number of electrodes if estimating for each case of 1 to 5 and will need to include appropriate priors. Depending on the plot of the data, it may be appropriate to do partial pooling for the effect of the number of electrodes on intercept and slope. Alternatively, if the plot shows there is an approximately linear influence of number of electrodes on intercept and/or slope then can make some modifications of your model to hopefully get reasonable estimates.

Numerous other improvements are possible but first need to determine the appropriate structure of the model given the data.

Daniel

[David Pedrosa]

Hi Daniel,

thanks for your suggestions. I have plotted the amplitude vs. the distance to target but I'm not sure if I understand what you are suggesting with adding the count of electrodes (although admittedly, I have not used ggplot2 too much). Do you have an example figure somewhere to understand what you mean? I know it is not much data which may lead to not very representative results but I would like to give it a try and I would be really interested to find out how you solved this.
Best

David

[Daniel Myall]

Below is an example where there is a geom_line() linking up data points for each subject over time and separate panel for each level of a factor (in this case cognitive classification, in your case would be number of electrodes). To get this you'd add '+facet_wrap(~n_electrodes)' to your plot.

In this case the group-level predictor of cognitive group has a strong effect on the intercept and a modest effect on the slope. 

I appreciate how difficult it is collecting this type of data. You'll still be able to fit a model that will give useful results, although the choice of priors and model structure become more important.

Regards,

Daniel

[David Pedrosa]

Alright, so I made a (scatter)-plot of the data just to give you a an idea of what I was talking about. You can see on the y-Axis the amplitude of the signal and on the x-axis the distance to the target. In my case I would say it is the opposite to what you were showing: a big amount of intercept variability in all 5 groups and if at all the number of electrodes having an influence on the slop. As you can also see, there is not much data, in fact for 2 electrodes it is just one subject. It may be thus not very promising but I would like to give it a try. In the worst case I learn more about regression and bayesian regression. What would your approach be to implement the number of electrodes as group-level predictor for the slope?

Thanks in advance

Cheers

David

--
You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/pLbnxj1pxXs/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+...@googlegroups.com.

[Daniel Myall]

Thanks, the plot makes it much easier to see the form of an appropriate model.

Firstly I'd demean the log amplitude: df$logamps <- df$logamp - mean(df$logamp)

Then try this model which has some minor changes to above:

data{
    int<lower=0> N;                     //number of total observations
    int<lower=1> J;                      //number of subjects
    int<lower=1,upper=J> s[N];     // suject ID
    vector[J] e;                            // group-level interactions (number of electrodes)

    vector[N] y;                           // outcome variable
    vector[N] x;                           // predictor(s)   
}

parameters{
    real beta_mu[2];                                    // group averages for intercept[1] and slope[2] 
    real<lower=0> sigma_y;                         // total error
    matrix[2,J] e_betaj;

    vector<lower=0>[2] sigma_b;    
    real gamma; // effect of electrode count on slope
}
  
transformed parameters{

    vector[N] y_hat;
    vector[J] beta0_raw;
    vector[J] beta1_raw;
    vector[J] beta0;
    vector[J] beta1;

    beta0_raw = beta_mu[1] + sigma_b[1] * e_betaj[1]';

    beta1_raw = beta_mu[2] + sigma_b[2] * e_betaj[2]';

    beta0 = beta0_raw; // No evidence of an effect of electrode count on Intercept and given limited data exclude estimation
    beta1 = beta1_raw + e * gamma;

   
    for (i in 1:N)
        y_hat[i] = beta0[s[i]] + beta1[s[i]] * x[i];

}
  
model {
    //Priors
    sigma_y ~ cauchy(0, 2.5);
    sigma_b ~ cauchy(0, 2.5);

    beta_mu ~ normal(0, 1); // Need tighter prior

    gamma ~ normal(0,1);
    to_vector(e_betaj) ~ normal(0, 1);
   
    y ~ student_t(4, y_hat, sigma_y);
}

Would be great if you could post the resulting inferences for sigma_y, sigma_b, beta_mu, and gamma. If still getting results that aren't sensible then making the change to a multivariate prior for the intercept and slope should hopefully fix it.

On Thursday, 6 October 2016 12:08:07 UTC+13, David Pedrosa wrote:

Alright, so I made a (scatter)-plot of the data just to give you a an idea of what I was talking about. You can see on the y-Axis the amplitude of the signal and on the x-axis the distance to the target. In my case I would say it is the opposite to what you were showing: a big amount of intercept variability in all 5 groups and if at all the number of electrodes having an influence on the slop. As you can also see, there is not much data, in fact for 2 electrodes it is just one subject. It may be thus not very promising but I would like to give it a try. In the worst case I learn more about regression and bayesian regression. What would your approach be to implement the number of electrodes as group-level predictor for the slope?

Thanks in advance

Cheers

David

[David Pedrosa]

Hi Daniel,

first of all thanks for all the suggestions. So I have tried both, using indexed gamma values for every of the five electrode counts (as Bob suggested) and I have also tried multiplying gamma with the number of electrodes. Results were both similar. But what actually changed a lot was using multivariate priors for intercept and slope (all according to chapter 6.12 in the manual) and especially when I restricted the group-level interaction to the slope (see the attached model). This is the output:

sd    2.5%    50%   97.5% n_eff Rhat betamu[1]   0.16  -3.76    -3.46  -3.13    3166    1 betamu[2]   0.01  -0.06    -0.04  -0.01    4014    1 sigma_y     0.07   0.01     0.07   0.09    4922    1 gamma_ec 0.02  0.00      0.01   0.06    3988    1 beta[2,1]    0.01 -0.03     -0.01   0.01    5000    1 beta[2,2]    0.01 -0.04     -0.02   0.00    4898    1 beta[2,3]    0.02 -0.09     -0.04   0.00    4648    1 beta[2,4]    0.01 -0.03     -0.02   0.00    5000    1 beta[2,5]    0.01 -0.03     -0.01   0.01    5000    1 beta[2,6]    0.01 -0.02      0.00   0.01    4854    1 beta[2,7]    0.02 -0.08     -0.04   0.01    4321    1 beta[2,8]    0.01 -0.03     -0.01   0.01    5000    1 beta[2,9]    0.01 -0.03     -0.01   0.01    5000    1 beta[2,10]  0.01 -0.04     -0.02   0.00    5000    1 beta[2,11]  0.01 -0.08     -0.06  -0.04    4740    1 beta[2,12]  0.01 -0.06     -0.04  -0.01    5000    1

 

Looking at the predictive posterior check and generating data with the values improves the fit and reduced the residual. How would you interpret the value of gamma? Does it even make sense to look at it or is enough to plot the results (e.g. slopes estimates vs. number of electrodes).

And another question. You introduced a comment saying that "beta_mu [...] // Need tighter prior" than normal(0, 1). Why do you think so?

 

Cheers
David

[Daniel Myall]

Hi David,

Re "beta_mu [...] // Need tighter prior" than normal(0, 1)" -- it was previously normal(0,5) which is quite wide, normal(0,1) was a change I'd made and the comment was just noting why. Tighter prior helps with sampling, although need to demean data so intercept becomes closer to 0.

The model is on the right track, and can see how useful a multivariate prior is. Overall mean negative intercept and negative slope for the relationship between signal and distance matches the data.

One thing is that gamma is the opposite sign to what it should be (and can't be any other way with the constraint that gamma has a lower bound of zero ;-). As you increase the number of electrodes the plot indicates that the slope becomes more negative. Hence as I assume ec[j] is positive, gamma should be negative. Thus change to real<upper=0> for gamma. A negative gamma can be interpreted as the amount that the slope becomes more negative per extra electrode, and definitely makes sense to look at. I assume the beta_raw's will be compensating for the current constraint, hence the other variable that would be useful to see is beta_raw as can match with plot to check that estimates are correct.

The results are looking promising, a few more tweaks to the model and I think will be close.

Regards,

Daniel

[David Pedrosa]

Hi Daniel,

I think I don't understand the point with the gamma. To be honest, this was unintended because I didn't mean to restrict the influence of the number of electrodes at all. Is it reasonable to do so because I can see the negative slope in the plots? If I want the influence of the number of electrodes why not allowing any values in a first instance?

I have been trying the two combinations with very different results for the group level averages of the slope (betamu[1,2]). While as I posted before the restriction to lower=0 makes the results comparable to when I do not include the group-level predictor, with real<upper=0> gamma_ec; the (group-average) slope turns positive. I have the suspicion that there is something wrong in my definition of betamu, because the individual slopes are comparable when I restrain gamma to either of the two. But I don't see any obvious problem there. Again, any help is appreciated.

Thanks in advance
Regards

David