Using Laplace approximation in user-defined functions

[Bayesian Baba]

Dear list:

My question is about using Laplace approximation of the GLMM likelihood in user-defined function.  

My question is very similar to the problem/solution posed/answered here:

https://github.com/stan-dev/stan/wiki/Max-Marginal-Optimization-(lmer)-Design

Following Bob's convention, if I assume 'alpha' to be the random effect, 'phi' to be the fixed effect, 'p(y | phi, alpha)' to be a member of the exponential family in

p(phi | y) = INTEGRAL p(phi,alpha | y) d.alpha,

then I can interpret 'log p(phi | y)' as the contribution of an individual to the log-likelihood in GLMMs.  Bob proposes a very general approach to solve this problem in Stan, but my problem is very structured.  To make the discussion concrete, in my case (with notational abuse)

p(phi,alpha | y)  = p1(y | alpha, phi) p2(alpha | Sigma),

where p1 is the Bernoulli likelihood with 

prob_success = 1 / [1 + exp{ - (x_t' * phi + z_t' * alpha)}]

and p2 is Multivariate Normal density with mean 0 and covariance matrix Sigma.  If I was in R, then for a given value of 'phi', 'Sigma', and 'y' I know how to find the Laplace approximation for the INTEGRAL p(phi,alpha | y) d.alpha:

(based on McCulloch & Searle 2001; pp 282)

1. calculate the gradient 

2. calculate the Hessian

3. use Fisher scoring or optim

4. plug the optimum in the expression for Laplace approximation, which is analytically tractable.

I can also do this in Rstan by using the 'optimizing' function; however, my needs are a bit different.

I want to generate random variables with the likelihood obtained after Laplace approximation. Specifically, I want define a function

functions {

  // takes care of the modified likelihood for one individual; assume I have sample size 1 currently.

  real laplace_approx_log (vector y, vector phi, matrix Sigma, vector startVal) {

        // return the result of 4. (i.e., analytic expression for Laplace approximation evaluated at the optimum with 'startVal' as starting value)

  } 

and use it later as

model {

   ...

   y ~ laplace_approx(phi, Sigma, startVal)  // I would like this function to be free of startVal, but this would work.

}

I see that the 'Point Estimation' section in the manual deals with general optimization problems.  

Question: I was wondering if there is a way to use the 'optimizing' function to solve my special problem using Hessians and gradients?  Unless I am missing something obvious, I could not find any documentation related to 'optimizing' in the manual. (I could also implement my own Newton-Raphson in the functions, but it would be amazing if I can use Stan's function.)

Thanks for your time and patience,

B.

[Seth Flaxman]

See page 15 of the rstan manual--there's a description of the function "optimizing". You can tell it what algorithm to use (algorithm = LBFGS, BFGS, or Newton) and whether or not it should calculate the hessian (hessian = TRUE).

(Did that answer your question?)

Seth

--
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.
For more options, visit https://groups.google.com/d/optout.

[Bayesian Baba]

Hi Seth:

Thanks for your quick reply.

I am not sure, if your suggestion would solve my problem.

I would like to use optimizing in the ".stan" script that would be sourced to RStan.  (May be I don't know how do use RStan functions like 'optimizing' inside the .stan files)

Here is a bit more concrete version:

I want to sample from an analytically intractable posterior with the log posterior of the following form:

real laplace_approx_log (vector y, vector phi, matrix Sigma, vector startVal) {

         // optimize a known function 'g' of (y, phi, Sigma, alpha) with respect to 'alpha'; say the optimum is alpha* = argmax g(y, phi, Sigma, alpha)

         // evaluate another known function 'f' at 'alpha*'; i.e., f(alpha*)  

        // return f(alpha*)

} 

The nice thing about Stan is that if I want to sample from a posterior with log-posterior 'laplace_approx_log', then all I need to do is the following in the .stan file

functions {

  real laplace_approx_log (vector y, vector phi, matrix Sigma) {

  // copy paste the above description ...

  } 

and sample from it later as

model {

   ...

   y ~ laplace_approx(phi, Sigma)  

}

I would like to optimize using Newton's method inside 'laplace_approx_log.'  Then my question is

is there a way to use 'optimizing'-like function in .stan file and inside the 'functions' chunk?

I don't think I can use an R function inside .stan file. 

Thanks again.

[Ben Goodrich]

On Wednesday, June 17, 2015 at 1:03:32 PM UTC-4, Bayesian Baba wrote:

I would like to optimize using Newton's method inside 'laplace_approx_log.'  Then my question is

is there a way to use 'optimizing'-like function in .stan file and inside the 'functions' chunk?

Not currently. We've been talking about it for a while. It would be less difficult to implement than the ODE stuff was but the ODE stuff had higher priority.
 

I don't think I can use an R function inside .stan file. 

Also not possible but any attempt along these lines would have to account for the derivatives, which an R function cannot do with Stan's custom scalar types.

Ben

[Bayesian Baba]

Hi Ben:

If I understand you correctly, then I need to implement the Fisher scoring/Newton-Raphson inside the 'functions' chunk if I were to do Laplace approx. for a multidimensional integral.

Thanks for all the help.

Best,

B.

[Ben Goodrich]

Something like that and even still, I don't think it would get the implicit differentiation right.

Ben

[Bayesian Baba]

Hi Ben:

I went ahead an programmed the laplace approx. using Newton-Raphson.

While this approx. is useful in a more general context, let's focus on logistic mixed effects model with normal random effects.  I want to weight the GLMM likelihood by 'nrep', so I am using laplace approx. as a method to marginalize out the random effects.  After some effort, I could do this in Stan.

However, I still have a few questions:  

- Why is a GLMM with 30 observations so slow? [I checked and found that Newton-Raphson iterations are pretty fast, so laplace approx is pretty fast] 

- This approx also fails sometimes with the sampler complaining that it is unable to find a good move (probably this is what Ben was referring to?)

- Am I expecting too much from stan?

- Or may be I am missing something obvious ...

I have attached the stan script. You can use the following command and script data

simList = structure(list(nobs = 30L, nfixef = 3L, nranef = 2L, ngroup = 4L, 

    nrep = 1, xmat = structure(c(-1, -1, -1, 1, 1, -1, -1, 1, 

    1, -1, -1, 1, -1, 1, -1, 1, -1, 1, 1, -1, 1, 1, 1, -1, 1, 

    1, -1, 1, -1, -1, 1, 1, 1, -1, -1, 1, 1, 1, 1, -1, -1, 1, 

    1, 1, 1, 1, -1, -1, 1, 1, 1, -1, -1, -1, 1, 1, 1, -1, 1, 

    1, -1, -1, 1, -1, 1, 1, 1, 1, -1, 1, -1, 1, 1, 1, 1, 1, -1, 

    1, 1, 1, 1, -1, -1, 1, 1, 1, -1, -1, 1, 1), .Dim = c(30L, 

    3L), .Dimnames = list(c("3", "13", "16", "17", "24", "1", 

    "2", "5", "6", "10", "14", "19", "20", "26", "29", "30", 

    "8", "9", "11", "12", "18", "23", "4", "7", "15", "21", "22", 

    "25", "27", "28"), NULL)), zmat = structure(c(-1, 1, 1, 1, 

    1, 1, -1, 1, 1, 1, -1, -1, 1, 1, -1, -1, 1, -1, 1, -1, -1, 

    1, 1, -1, 1, 1, -1, 1, 1, 1, 1, 1, -1, -1, -1, 1, -1, 1, 

    -1, 1, 1, 1, 1, 1, 1, -1, -1, -1, 1, 1, -1, -1, -1, 1, -1, 

    -1, 1, 1, -1, -1), .Dim = c(30L, 2L), .Dimnames = list(c("3", 

    "13", "16", "17", "24", "1", "2", "5", "6", "10", "14", "19", 

    "20", "26", "29", "30", "8", "9", "11", "12", "18", "23", 

    "4", "7", "15", "21", "22", "25", "27", "28"), NULL)), group = c(1L, 

    1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 

    3L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L), 

    yvec = c(1, 1, 1, 0, 0, 0, 1, 0, 0, 0, 1, 1, 0, 0, 0, 1, 

    1, 0, 1, 1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 1), pos1 = c(1, 6, 

    17, 23), pos2 = c(5L, 16L, 22L, 30L)), .Names = c("nobs", 

"nfixef", "nranef", "ngroup", "nrep", "xmat", "zmat", "group", 

"yvec", "pos1", "pos2"))

stanCode <- readChar("mod_lik.stan", file.info("mod_lik.stan")$size)

mdl <- stan(model_code = stanCode, data = simList, iter = 2000, warmup = 1000, chains = 1, thin = 5)

[Bob Carpenter]

I'm confused as to why you can't just code the "mixed effects"
logistic regression directly in Stan:

To make the discussion concrete, in my case (with notational abuse)

p(phi,alpha | y) = p1(y | alpha, phi) p2(alpha | Sigma),

where p1 is the Bernoulli likelihood with

prob_success = 1 / [1 + exp{ - (x_t' * phi + z_t' * alpha)}]

and p2 is Multivariate Normal density with mean 0 and covariance matrix Sigma.

What's the statistical operation you want to perform where
you need a Laplace approximation? Is it just because
you want a max marginal likelihood estimate? If so, you're probably
better of sticking with lme4 until we get all that working in Stan.

- 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.
> For more options, visit https://groups.google.com/d/optout.

> <mod_lik.stan>

[Bayesian Baba]

Hi Bob:

Sorry to be unclear.

"What's the statistical operation you want to perform where 

you need a Laplace approximation?  Is it just because 
you want a max marginal likelihood estimate?  If so, you're probably 
better of sticking with lme4 until we get all that working in Stan."

I would like to find the likelihood of 'phi' (fixed effects) and 'Sigma' (the covariance matrix for random effects).  As you correctly point, this is just the likelihood after marginalizing the random effects; however, unlike Lme4, I am not interested in maximizing the likelihood of 'phi' and 'Sigma'.  I am interested in sampling when the likelihood is modified as 

myLogLikelihood = nrep * marginalLogLikelihood // marginal LogLikelihood is what lme4 maximizes and would use with a usual GLMM...

In a Gaussian LMM, the marginaLogLikelihood is also a normal, so I was able to get this to work with Stan [very nicely actually; see code below].  For the non-Gaussian case, one needs to approximate the integral and I chose laplace approximation as quadrature would be slow. The laplace approx. part was attached in the stan file earlier. I think it works correctly, but I am not sure how to increase the speed.  

May be I am doing something ridiculously stupid and I am not able to catch my error.

Thanks again for your time.

Best,

B.

functions {

  // takes care of the modified likelihood

  real marg_approx_log (vector y, vector mu, real stdErr, matrix covMat, matrix zmat, real nrep) {

    matrix[num_elements(y), num_elements(y)] covY; 

    int tmp1; int tmp2;

    covY <- (zmat * covMat) * zmat' + stdErr * diag_matrix(rep_vector(1.0, rows(zmat)));

    return (nrep * multi_normal_log(y, mu, covY));

  } 

}

....

model {  

  // handle prior for rand. eff. prior and fixed eff.

  for (ii in 1:ngroup) {

    segment(yvec, pos1[ii], pos2[ii] - pos1[ii] + 1) ~ marg_approx(segment(mu, pos1[ii], pos2[ii] - pos1[ii] + 1), 

                                                                   stdErr, covRanef,

                                                                   block(zmat, pos1[ii], 1, pos2[ii] - pos1[ii] + 1, nranef),

                                                                   nrep);

  }

}

[Bob Carpenter]

What I don't understand is why you don't just code the
model directly in terms of all the parameters and then just
ignore the parameters you don't care about in the posterior
sample. The ease of marginalization is one of the big selling
points of MCMC.

- Bob

[Bayesian Baba]

Hi Bob:

I totally agree with you, but may be I am trying to do this because I don't know how to use Stan well.  This is what I want to do:

1. Generate random error:  stdErr ~ cauchy(0, 2.5);

2. Generate random effect covariance matrix:

     matrix[nranef, nranef] covRanef;

      sclRanef ~ cauchy(0, 2.5);

      corrRanef ~ lkj_corr(2);

       covRanef <- quad_form_diag(corrRanef, sclRanef);

3. Generate random effects:

  for (ii in 1:ngroup) {

    ranef[ii] ~ multi_normal(meanRanef, covRanef); 

  }

  

4. Generate fixed effects:

  stdErrFixef ~ cauchy(0, 2.5);  

  fixef ~ normal(meanFixef, stdErrFixef);

Steps 1-4 are usual GLMM steps, but the last step that samples data is different.

5. I want to draw data 'y' from a normal / bernoulli distribution raised to the power 'nrep'. For normal everything is tractable, so I can compare the results of Stan 

  real yHat[nobs];

  for (ii in 1:nobs) {

    yHat[ii] <- xmat[ii] * fixef + zmat[ii] * ranef[group[ii]]; // individual level mean 

  }

  // sample data

  yvec ~ normal(yHat, stdErr / nrep)  // is dividing by nrep correct (?) if my likelihood is normal likelihood raised to a power .. This seems to be true for the conditional 

and my own version based on laplace approx.  I am a bit surprised myself to see that Stan and Laplace approx. results don't match .... Laplace approx. matches the truth better.

Best,

S.

[Bob Carpenter]

> On Jun 18, 2015, at 8:24 PM, Bayesian Baba <resear...@gmail.com> wrote:
>
> Hi Bob:
>

> I totally agree with you, but may be I am trying to do this because I don't know how to use Stan well. This is what I want to do:
>
> 1. Generate random error: stdErr ~ cauchy(0, 2.5);
> 2. Generate random effect covariance matrix:
> matrix[nranef, nranef] covRanef;
> sclRanef ~ cauchy(0, 2.5);
> corrRanef ~ lkj_corr(2);
> covRanef <- quad_form_diag(corrRanef, sclRanef);
> 3. Generate random effects:
> for (ii in 1:ngroup) {
> ranef[ii] ~ multi_normal(meanRanef, covRanef);
> }
>
> 4. Generate fixed effects:
> stdErrFixef ~ cauchy(0, 2.5);
> fixef ~ normal(meanFixef, stdErrFixef);
>
> Steps 1-4 are usual GLMM steps, but the last step that samples data is different.
>
> 5. I want to draw data 'y' from a normal / bernoulli distribution raised to the power 'nrep'. For normal everything is tractable, so I can compare the results of Stan
>
> real yHat[nobs];
> for (ii in 1:nobs) {
> yHat[ii] <- xmat[ii] * fixef + zmat[ii] * ranef[group[ii]]; // individual level mean
> }
>
> // sample data
> yvec ~ normal(yHat, stdErr / nrep) // is dividing by nrep correct (?) if my likelihood is normal likelihood raised to a power .. This seems to be true for the conditional

Dividing by nrep is not correct. If you want to raise the normal
density to a power, you have to multiply on the log scale,
as in:

increment_log_prob(nrep * normal_log(yvec, yHat, stderr)

If some of the values yvec, yHat, or stderr are unknown, just
declare them as parameters.

In this way, we can do any kind of weighted regression (or other
model) you want, right in Stan.

Not quite sure why you'd want to do this other than if you're doing
something like EM or have importance weights or something from the
outside.

- Bob