23 thoughts on “Bias remaining after adjusting for pre-treatment variables. Also the challenges of learning through experimentation.”

1. I think that students learn more from a single class like this than a semester of theorem-proof style statistics.

   For me, the first big revelation in statistics was that you CAN simulate stuff: it is pretty easy and very instructive. The second, which came years later, is that you SHOULD simulate almost always. At the very least, test inference with simulated data.

   Reply ↓

2. A good way to get lots of bias simply is with a step function, eg:
   
z = rbinom(n, 1, .1 + .8 * (x > .4))
y = 0.2 + 0.1*(x > .4) + 0.5*z + rnorm(n, 0, 0.4)


   Reply ↓

   • I think that’s just more the extreme dependence of z on x and less the step function in y. Andrews example only slowly goes from 25% to 75% Z=1 with x. Using a more comparable setup I find the result to not be biased in b_z at all but only for b_x and the intercept. I’m not sure why.
     n <- 1000000
     x <- runif(n, -1, 1)
     z .4) + 0.5*z + rnorm(n, 0, 0.4)

     fake <- data.frame(x, y, z)

     fit <- lm(y ~ x + z, data=fake)
     summary(fit)
     estimate |t|)
     (Intercept) 0.3507742 0.0006784 517.097 <2e-16 ***
     x 0.3149222 0.0007625 412.990 <2e-16 ***
     z -0.0013109 0.0008955 -1.464 0.143

     Reply ↓

     • Well that didn’t look like I thought it would, but I hope it’s still readable.

       Reply ↓

3. I think the issue is that x ranges over -1,1

   If the true function is f(x) then you can replace it with a 1st order Taylor series and as long as the Taylor series does an ok job over the range [-1,1] your fit will be ok. Sort of.

   So you had to construct an example model where over This range stuff didn’t look “fairly linear”, or have some symmetry property that resulted in automatically averaging to 0 or similar…

   Now I’m going to go to Desmos and see if I can construct some graphs that illustrate this and see if I’m right.

   Reply ↓

   • https://www.desmos.com/calculator/wenvfufnmf

     That shows all the different x functions you used, plus has a line which is parameterized by two parameters a, b that are associated with sliders you can fiddle with.

     over the range [-1,1] only the 10x^3 function deviates a LOT from a straight line. The other functions could be approximated as a line segment with some small deviations through the region [-1,1] by fiddling with the sliders.

     Reply ↓

     • Daniel:

       Indeed, it’s kind of like this amusing example from a few years back.

       Reply ↓

4. Doesn’t stan_glm have default priors on the coefficients that actually make this a weakly regularized linear model? If so I would have expected that there is a small but desirable bias towards 0 for both. Am I wrong about stan_glm, am I misunderstanding what the default priors do, or is it that this particular bias is just too small to be visible in that example?

   Now that I’m writing this I see that you’re switching to `glm` in the simulations though `stan_glm` is used in the question. Is it intended? Would the results change with stan_glm? I can’t run a test right now.

   Reply ↓

   • Étienne:

     Oh, good catch. I meant to switch entirely to “lm” in the above code to avoid having to think about the Bayesian angle (indeed, the “stan_glm” estimate is not, strictly, unbiased). In practice the prior won’t matter much here, but it’s a cleaner problem to just stick with least squares.

     Reply ↓

5. >140 standard errors away from zero . . . yeah, that’s a real bias!

   Yeah but it’s only a 0.4749/0.5 = ~5% bias away from the true effect. Seems like it would be good to crank the bias so that it becomes large in magnitude or flips sign etc.

   Reply ↓

   • Anon:

     Agreed. That’s why I wrote above, “I’m sure that with a little bit of thought I could tweak the example to make the bias much bigger.”

     Reply ↓

6. So how does this square with Pearl’s framework of “closing the back-door path” by conditioning on x? I figured out the original question by sketching out the DAG, and thought you simply have to condition on x. But are you saying you actually have to condition on the “correct function of x on y ” to completely close that back-door?

   And how would you go about doing that in an lm formulation?

   I tried increasing the n of the dataset to even out of some of the randomness, increasing the range of x to allow for more effect of the non-linearity following Daniel’s comment and setting a seed. With an x^2 generating function, the lm(y ~ x + z) seems to be doing just as well as an lm(y ~ x^2 + z) formulation. But with an x^3 generating function the z estimation is way off, even if you specify lm(y ~ x^3 + z).

   #Test with lm(y ~ x^2…) vs. as here lm(y ~ x…)
   library(“arm”)
   n_loop <- 10^4
   set.seed(1234)
   estimate <- rep(NA, n_loop)
   for (loop in 1:n_loop) {
   n <- 10^4
   x <- runif(n, -3, 3)
   z <- rbinom(n, 1, invlogit(x))
   y <- 0.2 + 10*x^2 + 0.5*z + rnorm(n, 0, 0.4)
   fake <- data.frame(x, y, z)
   fit <- lm(y ~ x + z, data=fake)
   estimate[loop] <- coef(fit)["z"]
   }
   print(c(mean(estimate), sd(estimate)/sqrt(n_loop)), digits=2)

   Lastly, could you explain the formula behind the statement "140 standard errors away from zero"?

   Reply ↓

   • Jens:

     Pearl’s framework, or what I’ve seen of it, does not address misalignment of the data generated process and the fitted model. The same is true of how we present causal inference in our BDA book: in the example of the above post, the treatment assignment is ignorable given x, and so the model for y needs to condition on x in order for it to be the appropriate likelihood. But these general principles don’t dictate what particular family of parametric models should be used—or, in the case of general nonparametric models, they don’t dictate what regularizations or constraints should be used.

     In short, yes, you will have bias if you adjust for the wrong model. In real life, we will always be adjusting for the wrong model, so we always have bias, as is generally the case in regression modeling. In some settings, the bias from fitting the wrong model can be larger than various other biases of design and measurement.

     The other relevant factor in this problem is the quantity of interest. When the goal is to estimate an average treatment effect, the bias correction comes from adjusting for average values of x in the two groups. As Daniel Lakeland points out in comments, if you’re performing a linear adjustment, this bias will be most problematic when the underlying relationship between y^0 and x is far from linear in the range of the data at hand, and that depends on all sorts of things in the model. If you fit a regression model that’s of the same family as the data generating process, it should work fine.

     Finally, my calculation is (0.5 – 0.4749) / 0.0032, which is 7.8, rather than 140. I don’t know where that 140 came from! I fixed it in the above post.

     Reply ↓

     • When I complain about economists not being interested in dynamics, or people doing regression discontinuity by fitting a different line on each side of some time point, or fitting a model to an outcome in a particular single time point some time after application of a treatment (5 yr survival rates, or income 20 years after kindergarten intervention or whatever) it’s all in some ways similar to this issue.

       Social scientists tend to view the whole world as a bunch of lines going up or lines going down. This is all that the world is, so the key to doing an analysis is to merely include the right dimensions (the right variables). They seem to act as if this gives them carte Blanche to claim findings on the basis of the coefficients in their models…. It drives me bananas.

       This concept that “everything is linear” makes sense **at the differential level** that is, when independent variables can change only over very small ranges. Hence in a differential equation… You can say we know the value of these variables today, and we know the rate of change today, and so we can predict the values tomorrow. But the rates of change themselves change through time , often nonlinearly, and so the resulting curves do whatever… Oscillate, bend, become attracted nonlinearly to certain values, chaotically bounce between attractors…

       Just take the growth in height of a male child between birth and three years. https://www.cdc.gov/growthcharts/data/set1clinical/cj41c017.pdf

       Growth is extremely fast initially, and the rate decays through time down to a kind of constant. Of course linear length growth still means nonlinear weight growth which we just begin to see at the end of this timeframe but is more prominent in later time frames

       https://www.cdc.gov/growthcharts/data/set1clinical/cj41c021.pdf

       When outcomes depend on multiple factors in nonlinear dynamic ways even simple ones like child height, you need a model that takes this into account, because malnutrition in the early phase of growth can lead to decades of poor health later. That’s why public health people push for breastfeeding for 6 months.

       It’s not just about height, and weight, it’s about brain development, immune system health, allergies, autoimmune, asthma, etc etc.

       Anyway it’s great that we have this simple example but really this sort of thing should be the MAIN content of stats classes. The probability and soforth is kinda ancillary, functional forms, dynamics, approximation theory, control theory, chaos, and similar should be major components of our education.

       Reply ↓

       • With all these factors (which I agree matter), most machine learning models make no assumptions regarding functional form, so should be considered superior to fitting regression models, linear or otherwise – yes? I know your preference for dynamic models of generative processes, but they will be limited in trying to explain things like school attendance or earnings at age 30. Anatomy and biology might provide some structure for how growth varies with age, but we are looking at an example where the dependent variable seems to elude specifying a generative model. You could make some assumptions and develop one, but why do that if you could just use predictive models that are totally inductive rather than deductive?

       • Dale:

         You write, “most machine learning models make no assumptions regarding functional form, so should be considered superior to fitting regression models, linear or otherwise – yes?”

         The answer is . . . Maybe. As I wrote elsewhere in this thread, these general principles don’t dictate what particular family of parametric models should be used—or, in the case of general nonparametric models, they don’t dictate what regularizations or constraints should be used.

         For a nonparametric method, the analogy to “function form of the curve” is “regularizer.” A nonparametric method will be regularized, and you need to choose the regularizer (or estimate the regularizer from some family of regularizers, etc.). It depends on context what will work better. If you only have 100 data points, the best regularizer might be something that pretty much smooths your fit all the way down to a low-dimensional model.

       • Dale, the problem with machine learning is that it’s kind of taking a general purpose function approximator and throwing data at it and asking it to follow the patterns in the data. In general to be effective you need a large number of parameters (ie the fancy LLMs these days have tens of billions or a hundred billion parameters) and a large quantity of data. If you don’t have that, you need information about behavior that comes from somewhere other than the data. What Andrew calls regularization. Some of that regularization can come in the form of providing a prior on a function’s behavior, or an alternative is to provide a more limited functional form based on knowledge. The most limited functional form that isn’t a constant is a line, so in some sense doing a line fit to data is regularization. But sometimes fitting a line is just vastly too regularized, particularly if the shape of the curve is what’s interesting.

         If a boy child starts at one place in the growth chart and doesn’t track along approximately with his cohort of similar starting point children, we are going to start asking questions. Does the child gain height too quickly? Perhaps there’s a hormone issue. Or maybe the child gains height too slowly and there’s a digestive issue. What happens if the child starts out too quick and then becomes too slow? That’s a different issue. Or starts out too slow and becomes too quick?

         Suppose we have 4 children, they all start out same height, they all wind up in the same height at age 1.5, but one is doing “what’s normal”, one grew extra fast and then slow (r shaped), one grew slow and then fast (S shaped), and one grew along a straight line (which is abnormal).

         A linear regression of height at age 1.5 years vs some differences in the kids blood markers or nutrition or whatever would show that the differences don’t matter, the kid winds up the same height no matter what. If however you look at the shape of the curve, A doctor on the other hand would immediately say 3 of these children have serious health issues that need investigation.

         Too much of social science is about regression of outcome at one single time point vs some stuff that is thought to affect that outcome and not plausibly related to a *trajectory* whose shape might tell us a lot more.

7. I really like this example, as it also nicely teaches the difference between the extra-statistical model that a causal graph would encode, showing that z can logically be identified without bias given x, and the statistical model that leads to issues around specifying a functional form, such as residual bias, but also bias/variance trade-offs in the real-world exploration to find a “good” model.

   Reply ↓

8. In practice, what is the solution to this? Propensity score models? including a bunch of splines and interactions that you *might* not need?

   Reply ↓

   • The solution is model building. Scientific models not just statistical fitting of some restricted class. And not just some default widening the restrictions to use splines or whatever.

     What’s the functional form you expect for y dependence on x? Why? Why is it independent of time, or is it? Is Y composed of many different pieces of y? Is this overall average from multiple subtypes of y which have different forms for each type?

     Assuming you’ve made the decision that y = f(x) + g(z(x)) + noise as in this model, then, yes you can now say something like x is in a compact support and so chebyshev polynomials, or splines, or radial basis functions all provide a complete basis for continuous functions on this space, so we can try them out… But in general measurement noise will cause higher order terms to explode, so we should provide a prior that describes how much wiggling the function should do.

     At each step we use domain knowledge and encode it in some form.

     Reply ↓

     • Sorry, I was actually thinking this was the thread from a few days ago, but I’ll use it as an example. The subject was how changes in maternity leave influence subsequent school attendance and earnings at age 30. I don’t see a meaningful model being developed along the lines you suggest (although the effort would prove interesting anyway). There is domain knowledge, but it doesn’t go much beyond specifying a large number of factors that would be relevant. That is why I was asking about data dredging versus trying to specify the generative process.

       Reply ↓

       • Dale, weather models are pretty good, they’re based on real physics of convection and soforth, but they can’t predict the weather 3 weeks from now much less a year or a decade. Heck, you can write down the full equations for the motion of a two bar physical pendulum and yet you can’t start with a 1 second video and then predict its future trajectory more than a few seconds out. Here’s a LEGO version of that experiment with 6 different realizations all synchronized to the “same” starting point.

         https://en.wikipedia.org/wiki/File:Double_pendulum_simultaneous_realisations.ogv

         Sometimes what we discover is that measuring what happens to a cohort of people 20-30 years from intervention and trying to relate it to a particular 6 week window at the beginning of their life is just hopeless and you shouldn’t be doing it and especially shouldn’t be looking for when the noise looks the way you like then claiming a strong effect that should be informing public policy.

         On the other hand, the weather model can predict for 2 weeks-ish, and maybe we can predict some stuff about what happens to kids between birth and age say 3 via this information about their mother staying home for some number of weeks. Can we justify the cost of the intervention based on this more reasonably predictable trajectory for the first few years of life? Maybe, or maybe not.

         It is totally plausible to me that we could study maternal leave and quantify very clear differences in things like average time the baby nurses on breastmilk, the average level of stress hormones in the baby and mother’s body, the relationship between those hormones and illnesses the child has between age 0 and 3, the relationship between illnesses in childhood and rates of allergy, blood pressures, or autoimmune diseases as an adult…. If each step in this causal chain produces a clear picture, then we can be more confident that the information we get about adult outcomes from increased maternal leave time represent real effects rather than merely noise trying to connect directly a 6 week treatment to a 25 year later income increment… :-|

         Also we should probably include into the calculus the fact that there may be some different intervention that has clear strong benefits in the short term and maybe we shouldn’t even be asking the question about a much less clear less good intervention.

9. Thank you!

   Reply ↓