Sampling and Optimizing returned very different results

[Xulong Wang]

Hi Stan,

In my model, generalized linear mixed model, optimizing() and sampling() functions returned very different estimation on parameter sigma. optimizing() estimated a tiny value (0.00173), while sampling() estimated a much bigger one (1.3).

Could you help me understand this?

data {

  int<lower=1> N;  // Sample

  int<lower=1> K;  // Category

  int<lower=1> D;  // Covariate

  row_vector[D] cov[N];  // Covariate

  cholesky_factor_cov[N] L;  // Covariance

  vector[N] g;  // genotype

  int<lower=1,upper=K> Ad[N];  // phenotype

}

parameters {

  real p;  // variant effect

  vector[D] beta;  // covariate effect

  real<lower=0.0001> sigma;

  ordered[K-1] c;  // cutpoints

  vector[N] z;

}

transformed parameters {

  vector[N] u;

  u <- L * z;

}

model {

  c ~ normal(3, 1);

  p ~ normal(0, 1);

  beta ~ normal(0, 1);

  sigma ~ inv_gamma(2, 1);

  z ~ normal(0, sigma);

  for (n in 1:N)

    Ad[n] ~ ordered_logistic(p * g[n] + cov[n] * beta + u[n], c);

}

[Andrew Gelman]

Don’t do real<lower=0.0001> sigma.  It should be real<lower=0>.  Also I recommend not using the inverse gamma prior distribution on sigma.  It’s non-intuitive, i.e. if you use it you could well be including prior information you don’t really want to include.  Also I don’t know how many cutpoints you have but it seems to me kinda weird to give them a common normal(3,1) prior.

--
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 Goodrich]

On Saturday, February 6, 2016 at 5:54:02 PM UTC-5, Xulong Wang wrote:

In my model, generalized linear mixed model, optimizing() and sampling() functions returned very different estimation on parameter sigma. optimizing() estimated a tiny value (0.00173), while sampling() estimated a much bigger one (1.3).

Could you help me understand this?

Plain optimization does not generally work well for what you call a generalized linear mixed model. The mode (assuming it is found) can be one number and the posterior mean or median a very different number. That is why the frequentist estimators of such models integrate out the group specific terms from the likelihood, which is then maximized over the common parameters.

Ben

[Xulong Wang]

Andrew, there are 3 cutpoints in my models. The normal(3, 1) prior was set after trying multiple methods, and I found the estimations are faster and reasonable this way. 

Also, change real<lower=0.001> to real<lower=0> did not resolve the issue. Other priors, such as uniform tends to push the estimation to the boundaries, and a normal(0, 1) prior tends to generate non-finite gradient issue more often. 

[Xulong Wang]

Ben, does not mean that I have to use sampling() in Stan for a generalized linear mixed model? 

I also found estimations of other parameters that I am mostly interested in, p and beta, are similar between optimizing() and sampling(). 

[Ben Goodrich]

On Sunday, February 7, 2016 at 5:16:14 AM UTC-5, Xulong Wang wrote:

Ben, does not mean that I have to use sampling() in Stan for a generalized linear mixed model? 

I also found estimations of other parameters that I am mostly interested in, p and beta, are similar between optimizing() and sampling(). 

Plain maximum likelihood doesn't have much justification for GLMMs because of the incidental parameter problem where the number of parameters to estimate increases as the number of levels of the grouping factors increases. For fixed levels, as the number of observations increases, maximizing the likelihood could be a consistent estimator, but otherwise only the invariance property of MLEs holds. That is why lme4::glmer and similar functions integrate out the group-specific unknowns.

Ben

[Xulong Wang]

Ben, in my model, I found the problem came from estimation of sigma, z ~ normal(0, sigma). What about simply put z ~ normal(0, 1) as a prior of z? The model will estimate sigma implicitly anyway while estimating z, just that I would not know what the exact value is.