SBC: Summarizing ranks

After SBC, we have a set of ranks for each parameter, and for a given parameter one can visualize the ECDF or dECDF as described in the excellent recent Graphical Tests… manuscript. If one has many parameters it might be desirable to summarize the ranks in some manner, and the original SBC manuscript describes (though with note on its poor performance) the use a simple \chi^2 quantity to achieve this. Before reading of that approach, I’d independently used a correlation-based quantity instead, and while this latter doesn’t have (at least to my admittedly analytics-naive knowledge) the nice analytic p-value derivation as the \chi^2, it’s straightforward to simulate from the “null” of uniform ranks to derive a reference distribution of correlations within which the observed correlation’s percentile can be obtained. Subsequently learning of the \chi^2 approach led me to be curious if the two approaches differed at all in their output (again, being rather analytics-naive I had a slight suspicion they might be different formulae for the same quantity), and if they did, how they might differ in their interpretation/use.

Here’s code demonstrating that they do indeed differ, which I’m most folks would have known they would, but included for explicitness in the quantities to which I’m referring (including note a third quantity that logit-transforms both sets of ranks before computing the correlation; something I had simply vague feelings might be better than the correlation on the raw ranks somehow):

library(GGally) #for easy pairs plot
library(tidyverse) #for all that is good and holy

#define a function to simulate some observed-vs-expected rank summaries
get_sims = function(
	n_ranks = 1e3 #number of SBC-derived ranks to simulate
	, n_sims = 1e4 #number of simulations to run
){
	expected_rank = (1:n_ranks)/n_ranks - 1/n_ranks/2
	out = tibble(
		i_sim = 1:n_sims
		, chisq = NA
		, r = NA
		, qr = NA
	)
	for(i in 1:n_sims){
		observed_rank = sort(runif(n_ranks))
		out$chisq[i] = sum(((observed_rank-expected_rank)^2)/expected_rank)
		out$r[i] = cor(observed_rank,expected_rank)
		out$qr[i] = cor(qlogis(observed_rank),qlogis(expected_rank))
	}
	return(out)
}

#compute the sims
sims = get_sims(1e3,1e4)

#visualize
(
	sims
	# apply some normalizing transforms
	%>% mutate(
		chisq = log(chisq)
		, r = qlogis(r)
		, qr = qlogis(qr)
	)
	# pairs plot
	%>% select(-i_sim)
	%>% GGally::ggpairs(
		#fix for font rendering
		upper = list(
			continuous = GGally::wrap("cor", family="sans")
		)
		#add alpha to points
		, lower = list(
			continuous = GGally::wrap("points", alpha = .01)
		)
	)
)

And here’s the viz that is spit out at the end (readers skipping the code: note there’s normalizing log and atanh transforms applied to the \chi^2 and correlation quantities, respectively):

So in posting here I’m mostly looking for folks’ thoughts on my half-baked alternatives to the \chi^2 approach. Are they outright nonsense?

3 Likes