A researcher who wishes to remain anonymous writes:
I am writing with a question about your article with Imbens, Why High-Order Polynomials Should Not Be Used in Regression Discontinuity Designs. In it, you discourage the use of high-order polynomials of the forcing variable when fitting models. I have a few questions about this:
(1) What are your thoughts about the use of restricted cubic splines (RCS) that are linear in both tails?
(2) What are your thoughts on the use of a generalized additive model with local regression (rather than with splines)?
(3) What are your thoughts on the use of loess to fit the regression models?
I wonder if the use of restricted cubic splines would be less susceptible to the difficulties that you describe given that it is linear in the tails.
My quick reply is that I wouldn’t really trust any estimate that jumps around a lot. I’ve seen too many regression discontinuity analyses that give implausible answers because the jump at the discontinuity cancels a sharp jump in the other direction in the fitted curve. When you look at the regression discontinuity analyses that work (in the sense of giving answers that make sense), the fitted curve is smooth.
The first question above is addressing the tail-wagging-the-dog issue, and that’s a concern as well. I guess I’d like to see models where the underlying curve is smooth, and if that doesn’t fit the data, then I think the solution is to restrict the range of the data where the model is fit, not to try to solve the problem by fitting a curve that gets all jiggy.
My other general advice, really more important than what I just wrote above, is to think of regression discontinuity as a special case of an observational study. You have a treatment or exposure z, an outcome y, and pre-treatment variables x. In a discontinuity design, one of the x’s is a “forcing variable,” for which z_i = 1 for cases where x_i exceeds some threshold, and z_i = 0 for cases where x_i is lower than the threshold. This is a design with known treatment assignment and zero overlap, and, yeah, you’ll definitely want to adjust for imbalance in that x-variable. My inclination would be to fit a linear model for this adjustment, but sometimes a nonlinear model will make sense, as long as you keep it smooth.
But . . . the forcing variable is, in general, just one of your pre-treatment variables. What you have is an observational study! And you can have imbalance on other pre-treatment variable also. So my main recommendation is to adjust for other important pre-treatment variables as well.
For an example, see here, where I discuss a regression discontinuity analysis where the outcome variable was length of life remaining, and the published analysis did not include age as a predictor. You gotta adjust for age! The message is: a discontinuity analysis is an observational study. The forcing variable is important, but it’s not the only thing in town. The big mistakes seem to come from: (a) unregularized regression on the forcing variable which randomly give you wild jumpy curves that pollute the estimate of the discontinuity, (b) not adjusting for other important pre-treatment predictors, and (c) taking statistically significant estimates and treating them as meaningful, without looking at the model that’s been fit.
We discuss some of this in Section 21.3 of Regression and Other Stories.
Personally I think the issue with high order polynomials is that people mainly use least squares and this is far too unregularized. Fortunately Bayesian methods give us the ability to regularize high order polynomials. Unfortunately this seems to be a secret that no one teaches. I’m putting it prominently in my book. My advice is this:
1) Use Chebyshev polynomial regression
2) Pick a finite length domain [a,b] that extends beyond the domain of the data sufficiently for any extrapolation you might want to do
3) Pick a high enough order that the order isn’t restricting the behavior of the function (I’d guess for a lot of problems this is 5-6th order or something)
4) Decide which representation you want: a) parameters are the coefficients or b) parameters are the value of the function at the chebyshev points (they are equivalent mathematically)
5) write a prior for the function as p(a_1)p(a_2)p(a_3)… Q(…)
where p(a_1) etc are individual simple factors for the coefficients (or for the chebyshev points) such as normal(a,b) type factors, and Q(…) is a factor that describes the *behavior of the entire function*. For example Q() could say that the function exceeds 1 throughout the region between x1 and x2, or that the mean squared second derivative is most likely less than some value. Or that the integral of the function between some values x1 and x2 is likely in some range… etc.
Why Chebyshev? Because it’s a complete orthogonal basis for functions on a closed interval [-1,1] which is easily transformed to any arbitrary [a,b], the closed interval [a,b] is all that is needed for the vast majority of applied problems, the coefficients will be generally fairly independent of each other (ie. low correlations in sampling) and it’s known that Chebyshev polynomials converge exponentially pointwise to any smooth function so you don’t need very many terms to express pretty much anything.
If you have a discontinuity in your RD design, add it in as a separate discontinuous basis function, basically f(x) = 1 if x > x* 0 otherwise.
If you’re putting a discontinuity in the function at point x* then you should probably have some information in your prior saying that the smooth polynomial isn’t doing anything crazy in that region. For example you might express information that the derivative at that point is near 0.
I still don’t understand why people thinking making all these arbitrary assumptions can replace a good theory and carefully collected data.
Look at the smoking and Parkinson’s thread. That correlation was apparently discovered around ~60 years ago and still remains unexplained. Eg, a pubmed search gives 1,323 results for “smoking parkinsons” going back to 1965.
If such methods were fruitful, wouldn’t one of the 1000+ papers have figured it out by now?
I think there’s two separate questions here. The one I’m really addressing is “how to represent a smooth function”. You can want smooth functions for lots of stuff other than the typical stuff you see here on the blog. For example, perhaps you want to calibrate a temperature vs voltage curve for a thermocouple to measure flame temperature in a boiler… or you want to estimate percentage body fat from some height and weight and other body size measurements, or you want to estimate the required calorie intake to maintain body weight for a pet dog, or estimate the trend in car accidents in the month surrounding the daylight saving time transition, plus any discontinuity on the day of the transition… etc.
You can use many models to make useful predictions. But afaict getting interpretable coefficients requires an entirely different approach.
These two different goals seems to have been conflated at some point, which has lead to people trying to do the impossible.
In analogy, I was playing 301 (darts) the other day, and someone who hadn’t played it before kept trying to double out when they had an odd number of points to get.* So their strategy never worked no matter how many rounds they had. After I realized what was going on I pointed it out and they were like oh yea, makes sense.
These efforts to interpret arbitrary coefficients/parameters seem like that.
*In the game you need to get to zero by hitting a double 1, 2, 3, etc (which would give 2, 4, 6, etc points). So if you have 15 points left, there is no double to hit. You first need to hit an odd number (less than 15) to get something divisible by 2.