Gaussian variance as a linear function of the mean

[jonatho...@gmail.com]

Hello,

I recently looked at some data where there was an obvious linear relationship between the mean and variance in a nested Gaussian model.  My initial attempts to model this relationship resulted in an enormous increase to warmup time and horrific ESS and Rhat figures.  For example the following simple model works well, 

data{

  int<lower=0> N_ ;   //number of data points e.g. N_ = 5000

  int<lower=0> n_b ;  //number of biologicaly unique proteins, e.g. n_b = 250

  int<lower=0, upper=n_b> bioID[N_] ; // Prot ID

  

  real lr[N_] ; //  observations

}

parameters{

  real<lower = 0> sigma; 

  real beta[n_b] ;

}

model{

  for( i in 1:n_b){

    beta[i] ~ normal(0,10);

  }

  

  sigma ~ normal(0,5)

  for(i in 1:N_){

    lr[i] ~ normal(beta[bioID[i]] , sigma ) ; 

  }

} // end model

However, just changing the distribution for the observations to 

      lr[i] ~ normal(beta[bioID[i]] ,  fabs(beta[bioID[i]])) ; 

    

or ideally to 

  lr[i] ~ normal(beta[bioID[i]] ,  fabs(beta[bioID[i]]) + sigma) ; 

causes substantial problems with the sampling.  

    

Does anyone have any suggestions for how this can be done more efficiently?  

Thanks,

Jonathon

[Ben B]

Howdy Jonathon

My first question is, do you have a reason to believe your means and standard deviations should be linked like this in your model? Because if there isn't some higher level parameter you know about that you want to estimate, I think it'd be better to not code that in (cause you could be forcing a questionable model on your data). Rather, just give every group a mean and a sigma, and then look at the posteriors to figure out if there is a link or not.

My second is, if you estimate the means and variances of each group, you'll probably want to do some sort of hierarchical model so things don't go crazy. Have you seen these examples? There's the 8 schools one in BDA3, or Bob has some baseball examples (https://youtu.be/qQFF4tPgeWI?t=1h33m43s, https://cran.r-project.org/web/packages/rstanarm/vignettes/pooling.html). 

Ben

[jonatho...@gmail.com]

Hi Ben,

Thanks for the suggestion!

In general, I agree that forcing the relationship will not be desirable.  However, I ran into a specific dataset where an exploratory analysis (just plotting group means against group sd's) was highly suggestive of a linear relationship.  So, I decided to play around and see if I could model it.  

Creating separate variance components for each group is problematic because the data is highly unbalanced with small sample sizes for many of the groups.  In general I'm happy with complete pooling, but a partial pooling model of the variance component is an interesting idea.  I'll give that a try.  

Do you expect that modeling the variance as a linear function of the mean is doomed to cause computational problems?

Thanks,

Jonathon  

[Ben B]

I'm not super familiar with modeling this, so maybe someone else can jump in. This sounds like something someone would have experience with.

Probly the easiest way to diagnose if a model is neurotic is some way (without digging into the math) is to just simulate data from the thing you're trying to fit, and then trying to fit it. Like, it's easy to generate data from a distribution normal(a/b, 1.0), but then trying to estimate 'a' and 'b' separately in a model is going to be hard (cause there are many 'a's and 'b's that are gonna give you good fits).

If Stan fits the simulated data easily, then the issue is probably something about the model not explaining your data very well that's causing your issues.

Hopefully that helps :D

Ben

[jonatho...@gmail.com]

That's a good point.  The convergence problems very well might have been because it was a bad idea to begin with.  I took a closer look at those plots and realized that the high variance was strongly associated with small sample sizes.  

I also took your advice and tried fitting a hierarchical model on the variance components, but I'm getting a lot of divergent transitions.  Does this specification look OK?

parameters{

  real beta_b[n_b] ;

  real<lower = 0> cScale ; //Cauchy scale parameter 

  real<lower = 0> sigmaP_b[n_b] ;

}

model{

    for(i in 1:n_b){

      beta_b[i] ~ normal(0, 10) ;

      sigmaP_b[i] ~ cauchy(0, cScale) ;

    }

    for(i in 1:N_){

      lr[i] ~ normal(beta_b[bioID[i]] , sigmaP_b[bioID[i]]) ;

      }

} // end model

I also tried putting a half-normal prior on the cScale parameter, but that didn't make a difference.  

Any advice would be greatly appreciated!

Thanks,

Jonathon

[Michael Betancourt]

http://mc-stan.org/documentation/case-studies/divergences_and_bias.html

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

[Ben B]

Yeah, build a non-centered form of the model like that case study by Michael shows you how to do. Your problem sounds exactly like what's going on there.

Make sure and send a note even if it just works and all your problems are solved. I'm curious. Good luck!

Ben

[andre.p...@googlemail.com]

It can get difficult to estimate individual mean-variance relations, it might be not 

the case, when you group the variances together.

Andre

[jonatho...@gmail.com]

Michael, thanks for sending me your exploration of divergence, it was very interesting.

Unfortunately, I'm still having a lot of trouble achieving convergence.  I tried specifying half-Cauchy(0,1) distributions and multiplying by cScale in the model.  Even with adapt_delta at .99 there were hundreds of divergent transitions. I then tried switching from half-Cauchy distributions to half-Normals.  This actually made things worse in terms of rhat and ess.   

So, I started to suspect that the problem might go beyond the geometry of the log-likelihood.  It turns out that the model with no pooling on the variance components also fails to converge.  Presumably, this is because many of the groups have only 1 or 2 data points making it impossible to estimate individual standard deviations.  In these cases I'm happy for the prior variance to be used, but maybe the difference in the scale of the prior is causing convergence problems?   

My current plan is to get the model with no pooling working before trying the hierarchical model again.  I'm hoping this can be done with the right selection of informative priors, but I'd be happy to hear any suggestions.  

Thanks again,

Jonathon

[Ben B]

I think you're trying the right sorts of things, but I don't have enough experience to say off the top of my head what's going wrong.

If you don't have restrictions on posting your model, a little example data, and a runfile, I wouldn't mind playing around with this a little myself.

Ben

[andre.p...@googlemail.com]

Jonathon,

you could avoid Stans normal-distribution which requires the absolute value and write

the normal-log-likelihood explicitly, eg. a squared version. The square should be more

easy to process by Stan autodiff than abs.

Andre

[jonatho...@gmail.com]

Hi Andre,

Thanks for the suggestion.  To clarify, are you suggesting the use of an explicit normal-log-likelihood for the mean parameters, variance components or both?

Jonathon

[jonatho...@gmail.com]

I would love to have some more eyes on this problem!  

I've attached the Rdata from the session I was working in, along with the Stan model that converges.  In this model variance is completely pooled with a handful of different variance components based on a categorical variable.  Separate means are estimated for each bioID.  So far, every attempt I've made to give each bioID its own variance component has resulted in convergence problems.  Please let me know if you find anything interesting.  

Thanks,

Jonathon

[andre.p...@googlemail.com]

Jonathon,

in your earlier post you use fabs:

 lr[i] ~ normal(beta[bioID[i]] ,  fabs(beta[bioID[i]]) + sigma)

I first about using the log-likelihood of the normal distribution direct to used

a squared version, instead take the absolut.

https://en.wikipedia.org/wiki/Normal_distribution#Maximum_entropy

The mean parameter remains the same, your sigma has to be squared, eg. the variance.

You would avoid the fabs. Also the variances of the normal can be added,

in your example you added the parameter interpreted as SD. At least I would square and add

them and then take the sqrt.

 lr[i] ~ normal(beta[bioID[i]] ,  sqrt(square(beta[bioID[i]]) + square(sigma)))

A problem might be also the multi-modality, eg. square(sigma) and square(-sigma) leading to same.

Thus I think the constraint on sigma remain necessary. But what about square(beta[bioID[i]]))?

Andre

[jonatho...@gmail.com]

Ah, I see what you're saying.  Thanks Andre, I'll give that a try.  

Regarding the hierarchical model I finally got something to work.  Using a non-centered parameterization with either a half-Cauchy or half-Normal didn't seem to help (actually it might have made things worse).  However, with a non-centered Inverse-Gamma distribution (shape =1), I was able to reduce the divergent transitions to 7.  Increasing adapt-delta got rid of the rest of them.  

I still never managed to get the unpooled model to run without convergence problems.  I know it's a bad idea to begin with since many levels are highly dependent on the priors, but I'm surprised that it causes so much trouble.

Jonathon

[Ben B]

Do the results make sense?

I just took the set of bioIDs defined for tagID == 1 and fit a little unpooled model there. I see divergent transitions but I wasn't able to make heads or tails of anything. All I managed to do was learn some new R commands haha.

Ben

[Ben B]

I'm not super sold on fitting a bunch of normals to this data. Maybe we got too caught up on trying to estimate the standard deviations and did too much path forking haha.

I attached two plots, one (bioIDs.png) is a bunch of boxplots showing the distribution of lrs for different bio/tag IDs. Hopefully it's not too confusing (I put the stuff you gave me in a dataframe). It seemed like there were an awful lot of data points sitting a long way away from the boxes.

The second (hist_lr_bioID_222) is the histogram of lr for one of the bioIDs (222).

For this bioID, I compute the mean as 0.8 and the standard deviation as 2.455. Now I get these are point estimates or whatever, but doesn't seem like a normal is fitting this. 0.8 is to the right of the first lump in the histogram. It seems like there's more to this than just a normal. Maybe this is where the problems in the model are coming from (vs. a sampling thing)? Does this make sense?

Ben

These are the commands I used to do the vis:

(bioIDs.png)

df %>%
  filter(bioID > 200 & bioID < 240) %>%
  ggplot(aes(bioID, lr)) +
  geom_boxplot(aes(group = bioID, fill = factor(tagID)))

(hist_lr_bioID_222.png)

df %>%
  filter(bioID == 222) %>%
  ggplot(aes(lr)) +
  geom_histogram()

[jonatho...@gmail.com]

You're right that a Gaussian model won't fit the data very well.  But I don't expect than anything will.  What you can see quite clearly in your plot is that a subset of the data appears to be coming from a totally different process.  Finding these clumps of data is one of the objectives we hope to achieve with our modeling efforts.  

I can't comment on whether or not this is contributing to the convergence issues.  It sounds like a reasonable guess.  

Thanks for taking a look!

Jonathon  

[andre.p...@googlemail.com]

You're right that a Gaussian model won't fit the data very well.  But I don't expect than anything will.  What you can see quite clearly in your plot is that a subset of the data appears to be coming from a totally different process.  Finding these clumps of data is one of the objectives we hope to achieve with our modeling efforts.  

I would choose a skewed t-distribution, the skewness parameter to be positive with low degree of freedom

and another student-t and having a mixture distribution together with an mixture parameter for each bioID. 

 

Added:

This fits your group mean:

qqplot(tapply(lr,bioID,mean),rnorm(10000,0,0.5))

and skewed t-distribution fits your errors:

lrN <- lr-tapply(lr,bioID,mean)[bioID]
fit<-gamlss(lrN ~ 1,family=ST5())

qqplot(lrN,rST5(1000,-.0712,exp(-0.97),0.007,exp(0.28)))

Then we have:

means ~ normal(0, 0.5);
lr ~ ST5(means[bioID] -.0712,exp(-0.97),0.007,exp(0.28));

The ST5 is the skewed t-dist. from this forum.

Andre

[jonatho...@gmail.com]

Hi Andre, 

Thanks for the suggestion!  I'm reading the Stan related paper on model evaluation right now http://www.stat.columbia.edu/~gelman/research/unpublished/loo_stan.pdf .  Once I have a good understanding of these ideas, I'll start playing around with optimizing the model fit and will definitely try out your changes.

Jonathon