Limited density variations: general question

[Cole Monnahan]

I'm wondering if there is any general advice for fixing models that have really limited density variations and are demonstrating random walk behavior. Michael lays out the problem really well in this paper, but RMHMC isn't available yet.

I ask because I have a GLMM with ~500 random effects, where both the centered and noncentered form have really low Neff/N for one hypervariance in particular. The NUTS criteria look good -- no divergences at 0.8 and reasonably low tree depths. The issue appears to be that a trajectory can only traverse a really small portion of the marginal distribution for the hypervariance.

Is there any general advice for dealing with a situation like this? Or maybe something else is going on in the model that I'm missing (see attached if interested).

Thanks,

Cole

[Ben Goodrich]

On Wednesday, August 17, 2016 at 2:29:45 PM UTC-4, Cole Monnahan wrote:

Is there any general advice for dealing with a situation like this? Or maybe something else is going on in the model that I'm missing (see attached if interested).

We have had really good success with reparameterizing the X matrix of a GL{M}M in terms of its QR factors; i.e. if X = Q * R run the model with Q in place of X and then multiply the coefficients on Q by the inverse of R after-the-fact. This is semi-automatic in rstanarm, just change

fit <- glmer(toF ~ State + Pods + (1 | PlantID) + (1 | Year) + (Pods-1 | PlantID),
             family = binomial, data = flowers)

to

library(rstanarm)
options(mc.cores = parallel::detectCores())
fit <- stan_glmer(toF ~ State + Pods + (1 | PlantID) + (1 | Year) + (Pods-1 | PlantID), family = binomial,
                  data = flowers, QR = TRUE, prior = NULL)

or do it yourself with rstan.

Ben

[Michael Betancourt]

This reparameterization resolves slow sampling due to linear correlations,

not the nonlinear correlations that can limit density variation, but it has

worked wonders for me on models with non identified effects and is definitely

worth trying.

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

[Cole Monnahan]

This is really interesting. I fixed a bug in my noncentered application and it is mixing MUCH better now.

I ran a quick comparison against the rstanarm code you posted and got this:

model	minESS	time	ESS.per.time
noncentered	70.1	135	0.52
rstanarm	92.2	455	0.202

 So it's mixing better but taking much longer. Is this what you'd expect? I had to switch prior=NULL with prior=normal() due to an error.

I haven't looked any deeper and would be willing to if you think it's worth it. The code used to make table is below. 

Thanks,

Cole

noncentered.fit <-

  stan(file='wildflower_nc.stan', data=data, init=inits, pars=pars, chains=1)

noncentered.sims <- extract(noncentered.fit, permuted=FALSE)

noncentered.ess <- monitor(noncentered.sims, print=FALSE)[, 'n_eff']

noncentered.time <- sum(get_elapsed_time(noncentered.fit))

noncentered.df <- data.frame(model='noncentered', minESS=min(noncentered.ess),

           time=noncentered.time,

           ESS.per.time= min(noncentered.ess)/noncentered.time)

## rstanarm code from Ben Goodrich; modified slightly

## https://groups.google.com/d/msg/stan-users/LRnibbozy_w/mVE64492AAAJ

flowers <- as.data.frame(data[c('year', 'plant', 'stage', 'Pods', 'toF')])

rstanarm.fit <- stan_glmer(toF ~ stage + Pods + (1 | plant) + (1 | year) + (Pods-1 | plant), family = binomial,

                  data = flowers, QR = TRUE, prior = NULL, chains=1)

rstanarm.fit2 <- rstanarm.fit$stanfit

rstanarm.sims <- extract(rstanarm.fit2, permuted=FALSE)

rstanarm.ess <- monitor(rstanarm.sims, print=FALSE)[, 'n_eff']

rstanarm.time <- sum(get_elapsed_time(rstanarm.fit2))

rstanarm.eff <- min(rstanarm.ess)/rstanarm.time

rstanarm.df <- data.frame(model='rstanarm', minESS=min(rstanarm.ess),

           time=rstanarm.time,

           ESS.per.time= min(rstanarm.ess)/rstanarm.time)

print(rbind(noncentered.df, rstanarm.df), digits=3)

[Ben Goodrich]

On Thursday, August 18, 2016 at 1:39:31 PM UTC-4, Cole Monnahan wrote:

This is really interesting. I fixed a bug in my noncentered application and it is mixing MUCH better now.

I ran a quick comparison against the rstanarm code you posted and got this:

model	minESS	time	ESS.per.time
noncentered	70.1	135	0.52
rstanarm	92.2	455	0.202

 So it's mixing better but taking much longer. Is this what you'd expect? I had to switch prior=NULL with prior=normal() due to an error.

stan_glmer() is slower than it should be because the csr_matrix_times_vector() function that evaluates the Z * b part of the model is slower than a sliderule. But the sampling should be good. What was the error with prior = NULL?

Ben

 

[Bob Carpenter]

> On Aug 18, 2016, at 7:39 PM, Cole Monnahan <mon...@uw.edu> wrote:
>
> This is really interesting. I fixed a bug in my noncentered application and it is mixing MUCH better now.

That's an example of Andrew's "folk theorem" in action. :-)

- Bob

[Cole Monnahan]

Ben: When I install rstanarm and choose to compile from source I get a series of errors (see below). So I used the CRAN version. When I specify prior=NULL I get this error:

Error in if (!prior_dist %in% unlist(ok_dists)) { : 

  argument is of length zero

In addition: Warning message:

In is.na(prior_df) : is.na() applied to non-(list or vector) of type 'NULL'

This happened on two computers. I can send you session info if you want. I don't use the package at the moment so it's not an issue for me.

Thanks,

Cole

rstanarm compiler errors:

* installing *source* package 'rstanarm' ...

** package 'rstanarm' successfully unpacked and MD5 sums checked

** libs

*** arch - i386

g++ -m32 -I"c:/PROGRA~1/R/R-32~1.4RE/include" -DNDEBUG -I"`"c:/PROGRA~1/R/R-32~1.4RE/bin/i386/Rscript" --vanilla -e "cat(system.file('include', 'src', package = 'StanHeaders'))"`" -DBOOST_RESULT_OF_USE_TR1 -DBOOST_NO_DECLTYPE -DBOOST_DISABLE_ASSERTS -DSUPPORT_OPENMP   -I"C:/Users/Cole/Rlibs/StanHeaders/include" -I"C:/Users/Cole/Rlibs/rstan/include" -I"C:/Users/Cole/Rlibs/BH/include" -I"C:/Users/Cole/Rlibs/Rcpp/include" -I"C:/Users/Cole/Rlibs/RcppEigen/include" -I"d:/RCompile/r-compiling/local/local323/include"     -O2  -mtune=core2 -c Modules.cpp -o Modules.o

In file included from C:/Rtools/mingw_64/x86_64-w64-mingw32/include/c++/bits/stl_algo.h:59:0,

                 from C:/Rtools/mingw_64/x86_64-w64-mingw32/include/c++/algorithm:62,

                 from C:/Users/Cole/Rlibs/Rcpp/include/RcppCommon.h:63,

                 from C:/Users/Cole/Rlibs/Rcpp/include/Rcpp.h:27,

                 from Modules.cpp:1:

C:/Rtools/mingw_64/x86_64-w64-mingw32/include/c++/cstdlib:178:10: error: expected unqualified-id before '__int128'

   inline __int128

          ^

In file included from C:/Users/Cole/Rlibs/Rcpp/include/RcppCommon.h:65:0,

                 from C:/Users/Cole/Rlibs/Rcpp/include/Rcpp.h:27,

                 from Modules.cpp:1:

C:/Rtools/mingw_64/x86_64-w64-mingw32/include/c++/limits:1405:35: error: template argument 1 is invalid

     struct numeric_limits<__int128>

                                   ^

C:/Rtools/mingw_64/x86_64-w64-mingw32/include/c++/limits:1479:44: error: template argument 1 is invalid

     struct numeric_limits<unsigned __int128>

                                            ^

make: *** [Modules.o] Error 1

Warning: running command 'make -f "Makevars" -f "c:/PROGRA~1/R/R-32~1.4RE/etc/i386/Makeconf" -f "c:/PROGRA~1/R/R-32~1.4RE/share/make/winshlib.mk" -f "C:\Users\Cole\Dropbox\emacs-lisp/.R/Makevars" SHLIB_LDFLAGS='$(SHLIB_CXXLDFLAGS)' SHLIB_LD='$(SHLIB_CXXLD)' SHLIB="rstanarm.dll" OBJECTS="Modules.o"' had status 2

ERROR: compilation failed for package 'rstanarm'

* removing 'C:/Users/Cole/Rlibs/rstanarm'

--
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/LRnibbozy_w/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+unsubscribe@googlegroups.com.

[Jonah Gabry]

On Thursday, August 18, 2016 at 5:27:43 PM UTC-4, Cole Monnahan wrote:

When I specify prior=NULL I get this error:

Error in if (!prior_dist %in% unlist(ok_dists)) { : 

  argument is of length zero

In addition: Warning message:

In is.na(prior_df) : is.na() applied to non-(list or vector) of type 'NULL'

This happened on two computers. I can send you session info if you want. I don't use the package at the moment so it's not an issue for me.

I think that's a bug. I'll get it fixed on GitHub.  

[Jonah Gabry]

Should be fixed on GitHub now 

https://github.com/stan-dev/rstanarm/commit/51a188106bf61b51e6e968456747eb4f2d5aa119