pairs() in RStan?

[Bob Carpenter]

Is there a way to have it plot all of the chains above
the diagonal? I just want the full set of correlation plots.
I'd prefer below the diagonal to just be empty.


When I try:

|
> pairs(fit, pars=c("k1","k2","alpha21","alpha12","gamma"), condition = list(1:4,1:4))
|

I get:

|
Error in if (!missing(bandwidth) && min(bandwidth) <= 0) stop("'bandwidth' must be strictly positive") :
missing value where TRUE/FALSE needed
In addition: Warning messages:
1: In min(x) : no non-missing arguments to min; returning Inf
2: In max(x) : no non-missing arguments to max; returning -Inf
3: In min(x) : no non-missing arguments to min; returning Inf
4: In max(x) : no non-missing arguments to max; returning -Inf
|


It works fine when I split the chains in the list, as in

condition = list(1:2, 3:4)

I also thought I'd give this example a spin:

|
> pairs(fit, pars=c("k1","k2","alpha21","alpha12","gamma"), condition = lp > median(lp))
|

but it says:

|
Error in pairs.stanfit(fit, pars = c("k1", "k2", "alpha21", "alpha12", :
object 'lp' not found
|

I couldn't get any form of the comparison working, using lp__ or the other variables,
so I'm probably missing something fundamental about R and this example.

- Bob

[Richard McElreath]

Here's my modified pairs method that I use with Rstan fits. It doesn't plot in the lower diagonal, but rather just displays the correlation value.

It does however require extracting the samples first, e.g. post <- extract(fit), then you can pass whichever samples you like to the function. You could also do the conditions on the samples first, before passing them in.

You'll see in the code how each portion of the plot is defined separately. If you don't want anything in the lower diagonal, just set lower.panel=NULL in the last line of the function.

mcmcpairs <- function( posterior , n=500 , cex=0.7 , pch=16 , adj=1 , ...) {

    panel.dens <- function(x, ...) {

        usr <- par("usr"); on.exit(par(usr))

        par(usr = c(usr[1:2], 0, 1.5) )

        h <- density(x,adj=adj)

        y <- h$y

        y <- y/max(y)

        abline( v=0 , col="gray" , lwd=0.5 )

        lines( h$x , y )

    }

    panel.2d <- function( x , y , ... ) {

        i <- sample( 1:length(x) , size=n )

        abline( v=0 , col="gray" , lwd=0.5 )

        abline( h=0 , col="gray" , lwd=0.5 )

        dcols <- densCols( x[i] , y[i] )

        points( x[i] , y[i] , col=dcols , ... )

    }

    panel.cor <- function( x , y , ... ) {

        k <- cor( x , y )

        cx <- sum(range(x))/2

        cy <- sum(range(y))/2

        text( cx , cy , round(k,2) , cex=2*exp(abs(k))/exp(1) )

    }

    pairs( posterior , cex=cex , pch=pch , upper.panel=panel.2d , lower.panel=panel.cor , diag.panel=panel.dens , ... )

}

[Ben Goodrich]

On Tuesday, November 11, 2014 1:12:05 AM UTC-5, Bob Carpenter wrote:

Is there a way to have it plot all of the chains above
the diagonal?  I just want the full set of correlation plots.
I'd prefer below the diagonal to just be empty.


When I try:

|
> pairs(fit, pars=c("k1","k2","alpha21","alpha12","gamma"), condition = list(1:4,1:4))
|

I get:

|
Error in if (!missing(bandwidth) && min(bandwidth) <= 0) stop("'bandwidth' must be strictly positive") :
  missing value where TRUE/FALSE needed

I think it is not prepared for overlapping sets of chains.

It works fine when I split the chains in the list, as in

  condition = list(1:2, 3:4)

I also thought I'd give this example a spin:

|
> pairs(fit, pars=c("k1","k2","alpha21","alpha12","gamma"), condition = lp > median(lp))
|

but it says:

|
Error in pairs.stanfit(fit, pars = c("k1", "k2", "alpha21", "alpha12",  :
  object 'lp' not found
|

Need to create lp first?

lp <- simplify2array(get_logposterior(fit))
 

I couldn't get any form of the comparison working, using lp__ or the other variables,
so I'm probably missing something fundamental about R and this example.

The only ones that it knows how to do automatically are the things in get_sampler_params()

Ben

[Bob Carpenter]

Totally my fault --- I missed the definition of lp above
the last example.

Would it be hard to just get it to accept everything for
above the diagonal (or below --- either way)?

I can try Richard's code --- that'll probably work for what
I need --- I don't need anything fancy now, I'm just trying
to reproduce the output from another package (Mueller and Sierra's
SoilR).

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