Skip to content

Ulam's pairs method does not filter pars parameters. [Easy fix?] #455

@jr-free

Description

@jr-free

Background Details:

OS: Windows 11
R version 4.4.3
cmdstanr version 2.36.0
Rethinking ver. 2.42

Problem Statement

The pairs method for ulam fails to correctly limit output to variables specified in the pars argument. Rather, regardless of the value supplied to pars, the resulting plot displays output for all model parameters (this gets quite ugly as you can imagine).

Minimum Reproducible Example

Consider the first example from the ulam man page:

library(rethinking)
data(chimpanzees)

m1 <- ulam(
    alist(
        pulled_left ~ bernoulli(theta),
        logit(theta) <- a + bp[condition]*prosoc_left  ,
        a ~ normal(0,4),
        bp[condition] ~ normal(0,1)
    ) ,
    data=d, chains=2, cores=1 , sample=TRUE )

pairs(m1)

The last line generates the following plot, as expected (note that pars is not specified):

Image

Supplying an argument to pars, however, does not change the output. Consider below:

# an idiot check to make sure the name(s) is(are) correct before going to pairs
m1@pars 
# [1] "a"  "bp"

# This *should* constrain the pairs plot to include only bp[1] and bp[2] -- it doesn't.
pairs(m1, pars='bp')

This code snippet yields the same plot as above. See below for a screenshot from my R session.

📝NOTE
This screen shot was taken using a cleaned R session with all previous plots discarded.

Image

Solution and Root Cause Analysis

The solution to the problem appears to be straight forward. There is a block in the ulam-class.R file (below) that is commented out which appears to correspond to filtering the relevant pars from the posterior samples. When this code is uncommented, the pairs plot functions as expected. Hence, uncommenting this block of code will solve the problem I am outlining here. However, this probably treats the symptom and not the cause.

rethinking/R/ulam-class.R

Lines 204 to 215 in ac1b3b2

setMethod("pairs" , "ulam" , function(x, n=200 , alpha=0.7 , cex=0.7 , pch=16 , adj=1 , pars , ...) {
#require(rstan)
if ( missing(pars) )
posterior <- extract.samples(x)
else
posterior <- extract.samples(x,pars=pars)
#if ( !missing(pars) ) {
# # select out named parameters
# p <- list()
# for ( k in pars ) p[[k]] <- posterior[[k]]
# posterior <- p
#}

By my uneducated guess-timation, it appears that the true offender here is the extract.samples method seen on line 209. From what I gather, extract.samples also is failing to recognize the pars parameter when it is called in the else branch in 209. Consider the screenshot below comparing the behavior of extract.samples against the output of the commented code block:

Image

Admittedly, this is about as deep as I went down this rabbit hole. I haven't fully diagnosed yet what the issue could be with extract.samples, but I'd like to think that it's minor (he said hopefully). This could potentially be an easy one to knock down.

Obligatory Disclaimer

I have the wonderful character trait of being exceptionally naive and absent-minded at times. There is the very distinct possibility that I completely misunderstood or overlooked something REALLY obvious (let me tell you about how I get on with minus signs....). So, apologies in advance if this issue is a false positive. I'm just trying to do my See Somethin', Say Somethin' part and contribute in my own little way.

Postscript / Edit

For some much needed catharsis, see the screenshot below where I manually updated and re-set the pairs method for ulam and re-issued the previous call. The plot is as expected. Everything is back in balance....except maybe those warnings 🤔

Image

Post-postscript

I believe this is the offending block of code from extract.samples. Again, unless I'm gaslighting myself something fierce, there should be logic somewhere before line 66 to handle the pars filtration. Looks like it just got missed (pars made it in the else block).

extract_post_ulam <-
function(object,n,clean=TRUE,pars,...) {
#require(rstan)
if ( missing(pars) & clean==TRUE ) pars <- object@pars
if ( !is.null(attr(object,"cstanfit")) ) {
# use posterior to extract draws and convert to array format
pr <- as_draws_rvars( attr(object,"cstanfit")$draws() )
p <- list()
for ( i in 1:length(pr) )
p[[ names(pr)[i] ]] <- draws_of( pr[[i]] )
} else
# assume old rstan fit
p <- rstan::extract(object@stanfit,pars=pars,...)
# get rid of dev and lp__

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions