Andrew and I have both written here about our Software Validation paper with Don Rubin. The last thing to add on the topic is that my website now has newly updated software to implement our validation method (go down to Research Software and there are .zip and .tar versions of the R package). If the software you want to validate isn’t written in R, you can always write an R function that calls the other program. Feel free to contact me (cook@stat.columbia.edu) with questions, comments, etc.

I’m not one to go around having philosophical arguments about whether the parameters in statistical models are fixed constants or random variables. I tend to do Bayesian rather than frequentist analyses for practical reasons: It’s often much easier to fit complicated models using Bayesian methods than using frequentist methods. This was the case with a model I recently used as part of an analysis for a clinical trial. The details aren’t really important, but basically I was fitting a hierarchcal, nonlinear regression model that would be used to impute missing blood measurements for people who dropped out of the trial. Because the analysis was for an FDA submission, it might have been preferable to do a frequentist analysis; however, this was one of those cases where fitting the model was much easier to do Bayesianly. The compromise was to fit a Bayesian model with a vague prior distribution.

Sounded easy enough, until I noticed that making small changes in the parameters of what I thought (read: hoped) was a vague prior distribution resulted in substantial changes in the resulting posterior distribution. When using proper prior distributions (which there are all kinds of good reasons to do), even if the prior variance is really large there’s a chance that the prior density is decreasing exponentially in a region of high likelihood, resulting in parameter estimates based more on the prior distribution than on the data. Our attempt to fix this potential problem (it’s not necessarily a problem if you really believe your prior distribution, but sometimes you don’t) is to perform preliminary analyses to estimate where the mass of the likelihood is. A vague prior distribution is then one that is centered near the likelihood with much larger spread.

We estimate the location and spread of the likelihood by capitalizing on the fact that the posterior mean and variance are a combination of the prior mean and variance and the “likelihood” mean and variance. Consider the model for multivariate normal data with known covariance matrix, and a multivariate normal prior distribution on the mean vector:

y| μ, Σ ~ N(μ, Σ)
μ ~ N(μ₀, Δ₀).

The posterior distribution of μ (where n represents the number of observations) is:

μ|y, Σ ~ N(μ_n, Δ_n), where

μ_n = (Δ₀^-1 + nΣ^-1)^-1(Δ₀^-1μ₀ +nΣ^-1y)

Δ_n^-1 = Δ₀^-1 + nΣ^-1.

Here y and Σ/n represent what I’m calling the likelihood mean and variance of μ. If we were unable to calculate them directly, we could do so by solving the above two equations for y and Σ/n, obtaining

Σ/n = (Δ_n^-1 – Δ₀^-1)^-1

y = (Σ/n) Δ₀^-1 (μ_n – μ₀ ) + μ_n.

A vague prior distribution for μ could then be something like N(y, Σ) or N(y, 20Σ/n). For more complicated models you could do the same thing. Let (μ₀, Δ₀) and (μ_n, Δ_n) represent the prior and posterior mean vector and covariance matrix of the model hyperparameters. First fit the model (with multivariate normal prior distribution on the hyperparameters) for any convenient choice of (μ₀, Δ₀), then use the equations above to estimate the location and spread of the likelihood for these parameters. This approximation relies on approximate normality of the hyperparameters. In large samples this should be true; in smaller samples transformations of parameters can make the normal approximation more accurate.

It’s also possible to check the accuracy of the likelihood approximation: Fit the model again using the estimated likelihood mean and variance as the prior mean and variance. The resulting posterior mean should approximately equal the prior mean, and the posterior variance should be about half the prior variance, if the likelihood approximation is good. If not, the process can be iterated: fit the model, estimate the likelihood mean and variance, use these as the prior mean and variance, fit the model again and compare the prior and posterior means and variances. Repeat until the prior and posterior means are approximately equal and the posterior variance is about half the prior variance. From here, a vague prior can be obtained by setting the prior mean to the estimated likelihood mean the prior variance equal to the estimated likelihood variance scaled by some large constant.

I’ve tried this method in some simulations and it seems to work, in the sense that after iterating the above procedure a few times you do obtain an estimated likelihood mean and variance that, when used as the prior mean and variance, lead to a posterior distribution with the same mean and half the variance. With simple or well-understood models, there are surely better ways than this to come up with a vague prior distribution, but in complex models this method could be a helpful last resort.

A wise statistician once told me that to succeed in statistics, one could either be really smart, or be Bayesian. He was joking (I think), but maybe an appropriate correlary to that sentiment is that to succeed in Bayesian statistics, one should either be really smart, or be a good programmer. There’s been an explosion in recent years in the number and type of algorithms Bayesian statisticians have available for fitting models (i.e., generating a sample from a posterior distribution), and it cycles: as computers get faster and more powerful, more complex model-fitting algorithms are developed, and we can then start thinking of more complicated models to fit, which may require even more advanced computational methods, which creates a demand for bigger better faster computers, and the process continues. As new computational methods are developed, there is rarely well-tested, publicly-available software for implementing the algorithms, and so statisticians spend a fair amount of time doing computer programming. Not that there’s anything wrong with that, but I know I at least have been guilty (once or twice, a long time ago) of being a little bit lax about making sure my programs actually work. It runs without crashing, it gives reasonable-looking results, it must be doing the right thing, right? Not necessarily, and this is why standard debugging methods from computer science and software engineering aren’t always helpful for testing statistical software. The point is that we don’t know exactly what the software is supposed output (if we knew what our parameter estimates were supposed to be, for example, we wouldn’t need to write the program in the first place), so if software has an error that doesn’t cause crashes or really crazy results, we might not notice.

So computing can be a problem. When fitting Bayesian models, however, it can also come to the rescue. (Like alcohol to Homer Simpson, computing power is both the cause of, and solution to, our problems. This particular problem, anyway.) The basic idea is that if we generate data from the model we want to fit and then analyze those data under the same model, we’ll know what the results should be (on average), and can therefore test that the software is written correctly. Consider a Bayesian model p(θ)p(y|θ), where p(y|θ) is the sampling distribution of the data and p(θ) is the prior distribution for the parameter vector θ. If you draw a “true” parameter value θ0 from p(θ), then draw data y from p(y|θ0), and then analyze the data (i.e., generate a sample from the posterior distribution, p(θ|y)) under this same model, θ0 and the posterior sample will both be drawn from the same distribution, p(θ|y), if the software works correctly. Testing that the software works then amounts to testing that θ0 looks like a draw from p(θ|y). There are various ways to do this. Our proposed method is based on the idea that if θ0 and the posterior sample are drawn from the same distribution, then the quantile of θ0 with respect to the posterior sample should follow a Uniform(0,1) distribution. Our method for testing software is as follows:

1. Generate θ0 from p(θ)
2. Generate y from p(y|θ0)
3. Generate a sample from p(θ|y) using the software to be tested.
4. Calculate the quantile of θ0 with respect to the posterior sample. (If θ is a vector, do this for each scalar component of θ.)

Steps 1-4 comprise one replication. Performing many replications gives a sample of quantiles that should be uniformly distributed if the software works. To test this, we recommend performing a z test (individually for each component of θ) on the following transformation of the quantiles: h(q) = (q-.5)2. If q is uniformly distributed, h(q) has mean 1/12 and variance 1/180. Click here for a draft of our paper on software validation [reference: www.techloris.com], which explains why we (we being me, Andrew Gelman, and Don Rubin) like this particular transformation of the quantiles, and also presents an omnibus test for all components of θ simultaneously. We also present examples and discuss design issues, the need for proper prior distributions, why you can’t really test software for implementing most frequentist methods this way, etc.

This may take a lot of computer time, but it doesn’t take much more programming time than that required to write the model-fitting program in the first place, and the payoff could be big.