I’ve been working on a hierarchical model for election data which attempts to incorporate surveys (approx. 60,000 respondents in the entire set, O(1000) per context) and actual election results (approx. 150 million voters, O(3 million) per context). I am modeling both kinds of data as binomial with probabilities parameterized by demographic information.
In practice, I’m rescaling the election result data or else the Stan run takes too long on my laptop. I’m currently able to run in reasonable time after dividing all the election results (eligible voters, votes, votes for candidate) by 100. Which still leaves more than an order of magnitude difference in scale between the surveys and election results.
There have been various challenges! One which remains is trying to understand the validity/goodness of fit. I can do posterior-predictive checks on the various data-sets. I would also like to do some cross-validation but–unsurprisingly I guess?–loo has bad/Inf pareto-k for all the election results because, I think, they are way more influential on the fit than the survey data given the relatively smaller variance of the binomial distribution with such large numbers of trials.
Is there some standard way to deal with this? Look at the loo for each data-set separately or something?
Can you post your actual model, or a sketch of it? I wonder if what’s happening is that you have context-specific parameters but are trying to do loo while leaving out entire contexts?
Yeah, complicated indeed! I don’t have time to wade in, but my guess is that the election results are highly informative for context-specific parameters, such that each parameter is getting informed by just one or a handful of data points. Is this consistent with how this model works?
I don’t think it is important per se that there’s an imbalance between the informativeness of election data versus survey data, except inasmuch as you can’t count on the survey data to be the extra datapoints that add more information to each parameter. If each parameter were informed by numerous election-based datapoints, then the problem would likely disappear.
That’s a good thing to check. I’ll do that soon and reply. But I’m not convinced, particularly because the number of problematic pareto-k’s goes down as I shrink the election results. When I shrink them to an order of magnitude larger than the surveys, I get problematic pareto-k’s for all of them but when I shrink them by another order of magnitude, I only get about 5% of them being problematic. So there is definitely something going on with the spread in scale between the surveys and the election results.
But, for sure, I should demonstrate that more convincingly by just modeling the election results without the surveys and making sure that the pareto-k’s are all fine and it’s not just an issue with the election results themselves as you suggest.
Do you mean that you’re subsampling it? I don’t understand why just dividing numbers would change efficiency. If you can subsample, then your results should be valid (assuming the model can be subsampled).
One thing you can do is replace binomial_lpmf with binomial_lupmf to remove normalizing constants, or switch to sampling notation which does that automatically. Otherwise, you’re computing a ton of gamma functions here which will dominate the computations. Any loops that can be replaced with matrix operations are also a win, but I can’t quit see how to do that with the dot products the way you’ve defined them (I’m not sure it’s possible).
Lines like these
for (k in 1:N_Elections_House) {
ElexT_Ho_ps[k] = ACS_By_CD[Elections_House_CD[k]];
}
can be simplified to one-liners like:
ElexT_Ho_ps = ACS_By_CD[Elections_House_CD];
assuming that Elections_house_CD is N_Elections_House elements long. Overall, I’d suggest pulling all the common computations into functions and trying to dramatically shorten your variable names so the structure of the code is more clear. It’ll also make sure that your code is doing the right thing, which is hard to verify with this much cut-and-paste.
Vectorizing the log likelihoods would be a huge win for efficiency if you don’t need them for x-validation.
I’m only dividing the election results by something. The surveys stay as they are. So what’s changing is the difference in scale between the election results and the CCES and CPS data.
Yep. I was doing that before but the elections are large and so there were numerical issues (underflow?) when I used “binomial_lupmf”. But once I’ve rescaled them, maybe I can switch back to “lupmf”
Yeah. The loops in “transformed parameters” don’t have any simple way to rewrite them as matrix ops because of the inv_logit functions.
Lines like these
That’s a good point. Thanks!
Yeah. That’s tricky. This is being generated by Haskell code (mine) so I can fix some things by improving the code generator, but others are harder than they might seem. It’ s not cut and paste right now so much as the code generator doing the same thing for each data-set. But I can think about how to have it generate functions and then re-use them.
I’m not sure what you mean here. You mean the stuff in the “generated quantities” block? That’s not my efficiency issue. When I run the model itself, that’s not even present. I compute those later via cmdstan and the “generate quantities” mode. And that doesn’t take very long. Or do you mean the sampling statements?
The code is using non-centered parameterization, which can be inefficient for informative likelihood, which would explain the improved efficiency when making the likelihood less informative.’
Sorry for the long delay in answering. I’ve been experimenting…
That’s an interesting point! I find the centering/non-centering to be very confusing. In this model, how informative each context (state) is–which I understand to mean how much the observations in that context depend on being in that context rather than other parameters–varies a lot depending on how many predictors I include. So if I fit just the alphas, using no predictors, the centered parameterization works best and then starts to struggle as I add predictors. I think.
But I hadn’t really thought about it in the context of rescaling the election size.
There are also the contexts of each different data-set, the usual sense of context in a meta-analysis. But I am not treating those hierarchically, at least for now. But maybe I should?
I’ve been thinking more about this and I realize it’s a general point of difficulty. Is there some way in stan to build a vector which acts as a “view” into another vector? Or some way of managing these sorts of loops where something is being in re-indexed and this somehow stops the entire thing from being vectorized?
So I’ve done some of this but it’s hard to do as much as might be best for straightforward code because stan seems unhappy with an expression that evaluates to a vector being given as an argument to a function that expects a vector. Which makes sense, since it’s C++ underneath. I can go further in the direction of functionalization by making functions which are model specific, which is a little painful from the code-generation standpoint but might help clarify the code. I’m still not sure I can remove any loops, though I might be able to combine a few.
Which beings me to another question: Is it worth outsize work to combine loops? What is the relative cost of
vector v[N];
vector w[N];
for (k in 1:N) {
v[k] = f(...)
}
for (k in 1:N) {
w[k] = g(...)
}
vs
vector v[N];
vector w[N];
for (k in 1:N) {
v[k] = f(...)
w[k] = g(...)
}
?
Combining loops is hard from for me from the code-generation side, because I don’t usually have more than one line or loop in hand before I write out the relevant Stan code. But in the long term, I plan to have the code generator work a bit more globally and this would be one good reason. Though to do it optimally, I’d need to keep track of expression dependencies so as to do whatever loops I need before I use things computed in them…
If a is a vector then a[L:U] is the slice of a from L to U. For now, these are being allocated and copied, but that’s not usually the bottleneck because there’s no autodiff overhead there.
If you have a function declared as int foo(vector y), then you should be able to call foo(e) with any expression e that denotes a vector. It’s hard to say more without seeing the error message you’re getting. If e is a vector and foo(e) isn’t accepted, that’s a bug that we would be very grateful if you could report.
There’s almost no difference in speed between
for (j in 1:J) foo[j] = ...;
for (j in 1:J) bar[j] = ...;
and
for (j in 1:J) {
foo[j] = ...;
bar[j] = ...;
}
In fact, the first one can be faster because it often has better memory-locality properties and loops are cheap. It’s even better to not have loops at all, because evaluating all the foo[j] and bar[j] independently (instead of in a shared expression that can share computation and expression graphs in the autodiff library) is suboptimal.
Ah. I think this was my misunderstanding from code generation and whatever makes those dot products hard to vectorize. I will try to make a more coherent statement of what doesn’t work as I might hope. It’s also connected to what I meant by a “view” into a vector–not a slice, but a reindexing via an array[] int. Now that I’ve rewritten the parameter transformations as a function, it’s clear that I can pass any vector expression where a vector is expected, as you indicate.
I’ve included the code again, now with the transformation done mostly as a function–there’s no way to return 2 things from a function, right? Are there tuples of vectors in Stan?–and you’re right, it’s clearer and, I think because of some manual common-subexpression-elimination, quite a bit faster. I’m including it here in case there are more easy opportunities for optimization that I have missed!
I still have my two original questions:
These fits have high (“bad” or worse) pareto-k for the election-data as compared to the survey data. I don’t really understand what that means but my rough understanding is that it indicates that those points are in some sense disproportionately influential. That might make sense here, since the elections have many more “trials” (in the binomial sense) than the surveys. So should I just expect the high pareto-k? Or am I misunderstanding?
What does it mean if I have low N_eff but only in “latent” variables? That is, every variable I use for post-stratification and analysis has reasonable N_eff but some of the variables used in the fit, e.g., some alphaT_raw but not alphaT. I am curious about two things: can I rely on the values of alphaT even though alphaT_raw has low N_eff and, less coherently, might that information indicate specific places to look for useful re-parameterization?
Anyway, here’s the re-factored code. Apologies for not shortening variable names! For me, that makes it harder to follow since there are various data sources with similar data.
Thanks for all the great suggestions and thoughts so far!