In Think Stats (Chapter 7), I showed that the mean birth weight for first babies is lower than the mean for others by about 2 ounces, and that difference is statistically significant. But there is also a relationship between birth weight and mother's age, and the mothers of first babies tend to be younger. So the question is: if we control for the age of the mother, are first babies still lighter?
Several students in my class are working on projects involving multiple regression, so I want to use this question as an example. I will follow the steps I recommend to my students:
1) Before looking at relationships between variables, look at each variable in isolation. In particular, characterize the distributions and identify issues like outliers or long tails.
2) Look at the variables pairwise. For each pair, look at CDFs and scatterplots, and compute correlations and/or least squares fits.
3) If there seem to be relationships among the variables, look for ways to separate the effects, either by breaking the data into subsets or doing multiple regression.
So let's proceed.
One variable at a time
The variables I'll use are
- Mother's age in years,
- Birth weight in ounces, and
- first, which is a dummy variable, 1 for first babies and 0 for others.
It looks like there is a little skew to the right, but other than that, nothing to worry about.
We have 9038 records with valid weights, with this distribution:
The middle of the distribution is approximately normal, with some jaggies due to round-off. In the tails there are some values that are are certainly errors, but it is hard to draw a clear line between exceptional cases and bad data. It might be a good idea to exclude some extreme values, but for the analysis below I did not.
There are only two values for first, so it's not much of a distribution, but with categorical data, we should check that we have enough values in each bin. As it turns out, there are 4413 first babies and 4735 others, so that's just fine.
One pair at a time
To characterize the relationship between weight and first, we compare the CDF of weights for first babies and others:
There is some space between the distributions, and the difference in means is 2 ounces. It's not obvious whether that difference is statistically significant, but it is, with p < 0.001.
Similarly we can compare the CDF of mother's age for first babies and others:
Here there is clearly space between the distributions. The difference in means is 3.6 years, which is significant (no surprise this time).
It is not as easy to see the relationship between age and birthweight. I often tell my students to start with a scatterplot:
But in this case it's not much help. If there is a relationship there, it is hard to see. An alternative is to break the age range into bins and compute the mean in each bin. Here's what that looks like with 2-year bins:
We can't take the first and last points too seriously; there are not many cases in those bins. And ideally I should represent the variability in each bin so we have a sense of whether the apparent differences are real. Based on this figure, it looks like there is a relationship, but it might be nonlinear.
Pearson's coefficient of correlation between age and birthweight is 0.07, which is small but statistically significant. Spearman's coefficient of correlation is 0.10; the difference between the two coefficients is another warning that the relationship is nonlinear.
If we ignore the non-linearity for now, we can compute a least squares fit for birthweight as a function of age. The slope is 0.28, which means that we expect an additional 0.28 ounces per year of age. Since first mothers are 3.6 years younger than others, we expect their babies to be 1.0 ounces lighter. In fact, they are 2.0 ounces lighter, so the linear model of weight vs age accounts for 50% of the observed difference.
Multiple regression
So far I have been using the thinkstats libraries to compute correlation coefficients and least squares fits. But for multiple regression I am going to break out rpy, which is the Python interface to R, a "language and environment for statistical computing and graphics. It is a GNU project which is similar to the S language and environment which was developed at Bell Laboratories." To see how it works, you can download the code I used in this section: age_lm.py
The first model I ran is the same linear model we were just looking at:
weight = intercept + slope * age
In R's shorthand, that's:
weights ~ ages
Here are the results:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 109.28635 1.08775 100.470 < 2e-16 ***
ages 0.27926 0.04258 6.559 5.72e-11 ***
Multiple R-squared: 0.004738, Adjusted R-squared: 0.004628
The intercept is 109 ounces; the slope is 0.28 ounces per year, which what we saw before---comparing my implementation to R makes me more confident that R is correct :)
The standard error provides a confidence interval for the estimates; plus or minus 2 standard errors is (roughly) a 95% confidence interval.
The last column is the p-value, which is very small, indicating that the slope and intercept are significantly different from 0. The t-value is the test statistic used to compute the p-value, but I don't know why it gets reported; it doesn't mean much.
Since the p-values are so small, it might be surprising that R2 is so low, only 0.0047. But there is no contradiction. We can say with high confidence that there is a relationship between these variables; nevertheless, birth weight is highly variable, and even if you know the age of the mother, that does not reduce the variability by much.
To understand adjusted R2, consider this: if you add more explanatory variables to a model, R2 usually goes up even if there is no real relationship. The adjusted R2 takes this into account, which makes it more meaningful to compare models with a different number of variables.
Now let's add first as an explanatory variable, so the model looks like this:
weights ~ first + ages
Here are the results:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 110.62791 1.24198 89.073 < 2e-16 ***
first -1.11688 0.49940 -2.236 0.0253 *
ages 0.24708 0.04493 5.499 3.93e-08 ***
Multiple R-squared: 0.005289, Adjusted R-squared: 0.005069
The coefficient for first is -1.1, which means that we expect first babies to be 1.1 ounces lighter, controlling for age. The p-value for this estimate is 2.5%, which I consider borderline significant.
The coefficient for ages is about the same as before, and significant. And R2 is a little higher, but still small.
But remember that the relationship between weight and age is non-linear. We can explore that by introducing a new variable:
ages2 = ages^2
Now if we run this model:
weights ~ ages + ages2
We are effectively fitting a parabola to the weight vs age curve.
Estimate Std. Error t value Pr(>|t|)
(Intercept) 89.151237 4.407677 20.226 < 2e-16 ***
ages 1.898151 0.346071 5.485 4.25e-08 ***
ages2 -0.031002 0.006577 -4.714 2.47e-06 ***
Multiple R-squared: 0.00718, Adjusted R-squared: 0.00696
Estimates for both variables are significant. The coefficient for ages2 is negative, which means that the parabola has downward curvature, as expected. And R2 is a little bigger.
Finally, we can bring it all together:
weights ~ first + ages + ages2
With this model we get:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 91.07627 4.56813 19.937 < 2e-16 ***
first -0.80708 0.50373 -1.602 0.109
ages 1.79807 0.35163 5.113 3.23e-07 ***
ages2 -0.02953 0.00664 -4.447 8.80e-06 ***
Multiple R-squared: 0.007462, Adjusted R-squared: 0.007132
When we include the parabolic model of weight and age, the coefficient for first gets smaller, and the p-value is 11%.
I conclude that the difference in weight for first babies is explained by the difference in mothers' ages. When we control for age, the difference between first babies and others is no longer statistically significant. It is still possible that there is a small difference in weight for first babies, but this dataset provides little evidence for it.
I conclude that the difference in weight for first babies is explained by the difference in mothers' ages. When we control for age, the difference between first babies and others is no longer statistically significant. It is still possible that there is a small difference in weight for first babies, but this dataset provides little evidence for it.