This HTML version of is provided for convenience, but it is not the best format for the book. In particular, some of the symbols are not rendered correctly. You might prefer to read the PDF version, or you can buy a hardcopy from Amazon. Chapter 11 RegressionThe linear least squares fit in the previous chapter is an example of regression, which is the more general problem of fitting any kind of model to any kind of data. This use of the term “regression” is a historical accident; it is only indirectly related to the original meaning of the word. The goal of regression analysis is to describe the relationship between one set of variables, called the dependent variables, and another set of variables, called independent or explanatory variables. In the previous chapter we used mother’s age as an explanatory variable to predict birth weight as a dependent variable. When there is only one dependent and one explanatory variable, that’s simple regression. In this chapter, we move on to multiple regression, with more than one explanatory variable. If there is more than one dependent variable, that’s multivariate regression. If the relationship between the dependent and explanatory variable is linear, that’s linear regression. For example, if the dependent variable is y and the explanatory variables are x_{1} and x_{2}, we would write the following linear regression model:
where β_{0} is the intercept, β_{1} is the parameter associated with x_{1}, β_{2} is the parameter associated with x_{2}, and ε is the residual due to random variation or other unknown factors. Given a sequence of values for y and sequences for x_{1} and x_{2}, we can find the parameters, β_{0}, β_{1}, and β_{2}, that minimize the sum of ε^{2}. This process is called ordinary least squares. The computation is similar to thinkstats2.LeastSquare, but generalized to deal with more than one explanatory variable. You can find the details at https://en.wikipedia.org/wiki/Ordinary_least_squares The code for this chapter is in regression.py. For information about downloading and working with this code, see Section 0.2. 11.1 StatsModelsIn the previous chapter I presented thinkstats2.LeastSquares, an implementation of simple linear regression intended to be easy to read. For multiple regression we’ll switch to StatsModels, a Python package that provides several forms of regression and other analyses. If you are using Anaconda, you already have StatsModels; otherwise you might have to install it. As an example, I’ll run the model from the previous chapter with StatsModels: import statsmodels.formula.api as smf
live, firsts, others = first.MakeFrames()
formula = 'totalwgt_lb ~ agepreg'
model = smf.ols(formula, data=live)
results = model.fit()
statsmodels provides two interfaces (APIs); the “formula”
API uses strings to identify the dependent and explanatory variables.
It uses a syntax called patsy; in this example, the smf.ols takes the formula string and the DataFrame, live, and returns an OLS object that represents the model. The name ols stands for “ordinary least squares.” The fit method fits the model to the data and returns a RegressionResults object that contains the results. The results are also available as attributes. params is a Series that maps from variable names to their parameters, so we can get the intercept and slope like this: inter = results.params['Intercept']
slope = results.params['agepreg']
The estimated parameters are 6.83 and 0.0175, the same as from LeastSquares. pvalues is a Series that maps from variable names to the associated pvalues, so we can check whether the estimated slope is statistically significant: slope_pvalue = results.pvalues['agepreg']
The pvalue associated with agepreg is 5.7e11, which is less than 0.001, as expected. results.rsquared contains R^{2}, which is 0.0047. results also provides And results provides resid, a sequence of residuals, and fittedvalues, a sequence of fitted values corresponding to agepreg. The results object provides summary(), which represents the results in a readable format. print(results.summary())
But it prints a lot of information that is not relevant (yet), so I use a simpler function called SummarizeResults. Here are the results of this model: Intercept 6.83 (0)
agepreg 0.0175 (5.72e11)
R^2 0.004738
Std(ys) 1.408
Std(res) 1.405
Std(ys) is the standard deviation of the dependent variable, which is the RMSE if you have to guess birth weights without the benefit of any explanatory variables. Std(res) is the standard deviation of the residuals, which is the RMSE if your guesses are informed by the mother’s age. As we have already seen, knowing the mother’s age provides no substantial improvement to the predictions. 11.2 Multiple regressionIn Section 4.5 we saw that first babies tend to be lighter than others, and this effect is statistically significant. But it is a strange result because there is no obvious mechanism that would cause first babies to be lighter. So we might wonder whether this relationship is spurious. In fact, there is a possible explanation for this effect. We have seen that birth weight depends on mother’s age, and we might expect that mothers of first babies are younger than others. With a few calculations we can check whether this explanation is plausible. Then we’ll use multiple regression to investigate more carefully. First, let’s see how big the difference in weight is: diff_weight = firsts.totalwgt_lb.mean()  others.totalwgt_lb.mean()
First babies are 0.125 lbs lighter, or 2 ounces. And the difference in ages: diff_age = firsts.agepreg.mean()  others.agepreg.mean()
The mothers of first babies are 3.59 years younger. Running the linear model again, we get the change in birth weight as a function of age: results = smf.ols('totalwgt_lb ~ agepreg', data=live).fit()
slope = results.params['agepreg']
The slope is 0.0175 pounds per year. If we multiply the slope by the difference in ages, we get the expected difference in birth weight for first babies and others, due to mother’s age: slope * diff_age
The result is 0.063, just about half of the observed difference. So we conclude, tentatively, that the observed difference in birth weight can be partly explained by the difference in mother’s age. Using multiple regression, we can explore these relationships more systematically. live['isfirst'] = live.birthord == 1
formula = 'totalwgt_lb ~ isfirst'
results = smf.ols(formula, data=live).fit()
The first line creates a new column named isfirst that is True for first babies and false otherwise. Then we fit a model using isfirst as an explanatory variable. Here are the results: Intercept 7.33 (0)
isfirst[T.True] 0.125 (2.55e05)
R^2 0.00196
Because isfirst is a boolean, ols treats it as a categorical variable, which means that the values fall into categories, like True and False, and should not be treated as numbers. The estimated parameter is the effect on birth weight when isfirst is true, so the result, 0.125 lbs, is the difference in birth weight between first babies and others. The slope and the intercept are statistically significant, which means that they were unlikely to occur by chance, but the the R^{2} value for this model is small, which means that isfirst doesn’t account for a substantial part of the variation in birth weight. The results are similar with agepreg: Intercept 6.83 (0)
agepreg 0.0175 (5.72e11)
R^2 0.004738
Again, the parameters are statistically significant, but R^{2} is low. These models confirm results we have already seen. But now we
can fit a single model that includes both variables. With the
formula Intercept 6.91 (0)
isfirst[T.True] 0.0698 (0.0253)
agepreg 0.0154 (3.93e08)
R^2 0.005289
In the combined model, the parameter for isfirst is smaller by about half, which means that part of the apparent effect of isfirst is actually accounted for by agepreg. And the pvalue for isfirst is about 2.5%, which is on the border of statistical significance. R^{2} for this model is a little higher, which indicates that the two variables together account for more variation in birth weight than either alone (but not by much). 11.3 Nonlinear relationshipsRemembering that the contribution of agepreg might be nonlinear, we might consider adding a variable to capture more of this relationship. One option is to create a column, agepreg2, that contains the squares of the ages: live['agepreg2'] = live.agepreg**2
formula = 'totalwgt_lb ~ isfirst + agepreg + agepreg2'
Now by estimating parameters for agepreg and agepreg2, we are effectively fitting a parabola: Intercept 5.69 (1.38e86)
isfirst[T.True] 0.0504 (0.109)
agepreg 0.112 (3.23e07)
agepreg2 0.00185 (8.8e06)
R^2 0.007462
The parameter of agepreg2 is negative, so the parabola curves downward, which is consistent with the shape of the lines in Figure 10.2. The quadratic model of agepreg accounts for more of the variability in birth weight; the parameter for isfirst is smaller in this model, and no longer statistically significant. Using computed variables like agepreg2 is a common way to fit polynomials and other functions to data. This process is still considered linear regression, because the dependent variable is a linear function of the explanatory variables, regardless of whether some variables are nonlinear functions of others. The following table summarizes the results of these regressions:
The columns in this table are the explanatory variables and the coefficient of determination, R^{2}. Each entry is an estimated parameter and either a pvalue in parentheses or an asterisk to indicate a pvalue less that 0.001. We conclude that the apparent difference in birth weight is explained, at least in part, by the difference in mother’s age. When we include mother’s age in the model, the effect of isfirst gets smaller, and the remaining effect might be due to chance. In this example, mother’s age acts as a control variable; including agepreg in the model “controls for” the difference in age between firsttime mothers and others, making it possible to isolate the effect (if any) of isfirst. 11.4 Data miningSo far we have used regression models for explanation; for example, in the previous section we discovered that an apparent difference in birth weight is actually due to a difference in mother’s age. But the R^{2} values of those models is very low, which means that they have little predictive power. In this section we’ll try to do better. Suppose one of your coworkers is expecting a baby and there is an office pool to guess the baby’s birth weight (if you are not familiar with betting pools, see https://en.wikipedia.org/wiki/Betting_pool). Now suppose that you really want to win the pool. What could you do to improve your chances? Well, the NSFG dataset includes 244 variables about each pregnancy and another 3087 variables about each respondent. Maybe some of those variables have predictive power. To find out which ones are most useful, why not try them all? Testing the variables in the pregnancy table is easy, but in order to use the variables in the respondent table, we have to match up each pregnancy with a respondent. In theory we could iterate through the rows of the pregnancy table, use the caseid to find the corresponding respondent, and copy the values from the correspondent table into the pregnancy table. But that would be slow. A better option is to recognize this process as a join operation as defined in SQL and other relational database languages (see https://en.wikipedia.org/wiki/Join_(SQL)). Join is implemented as a DataFrame method, so we can perform the operation like this: live = live[live.prglngth>30]
resp = chap01soln.ReadFemResp()
resp.index = resp.caseid
join = live.join(resp, on='caseid', rsuffix='_r')
The first line selects records for pregnancies longer than 30 weeks, assuming that the office pool is formed several weeks before the due date. The next line reads the respondent file. The result is a DataFrame with integer indices; in order to look up respondents efficiently, I replace resp.index with resp.caseid. The join method is invoked on live, which is considered the “left” table, and passed resp, which is the “right” table. The keyword argument on indicates the variable used to match up rows from the two tables. In this example some column names appear in both tables,
so we have to provide rsuffix, which is a string that will be
appended to the names of overlapping columns from the right table.
For example, both tables have a column named race that encodes
the race of the respondent. The result of the join contains two
columns named race and The pandas implementation is fast. Joining the NSFG tables takes less than a second on an ordinary desktop computer. Now we can start testing variables. t = []
for name in join.columns:
try:
if join[name].var() < 1e7:
continue
formula = 'totalwgt_lb ~ agepreg + ' + name
model = smf.ols(formula, data=join)
if model.nobs < len(join)/2:
continue
results = model.fit()
except (ValueError, TypeError):
continue
t.append((results.rsquared, name))
For each variable we construct a model, compute R^{2}, and append the results to a list. The models all include agepreg, since we already know that it has some predictive power. I check that each explanatory variable has some variability; otherwise the results of the regression are unreliable. I also check the number of observations for each model. Variables that contain a large number of nans are not good candidates for prediction. For most of these variables, we haven’t done any cleaning. Some of them are encoded in ways that don’t work very well for linear regression. As a result, we might overlook some variables that would be useful if they were cleaned properly. But maybe we will find some good candidates. 11.5 PredictionThe next step is to sort the results and select the variables that yield the highest values of R^{2}. t.sort(reverse=True)
for mse, name in t[:30]:
print(name, mse)
The first variable on the list is Similarly prglngth has useful predictive power, but for the office pool we assume pregnancy length (and the related variables) are not known yet. The first useful predictive variable is babysex which indicates whether the baby is male or female. In the NSFG dataset, boys are about 0.3 lbs heavier. So, assuming that the sex of the baby is known, we can use it for prediction. Next is race, which indicates whether the respondent is white, black, or other. As an explanatory variable, race can be problematic. In datasets like the NSFG, race is correlated with many other variables, including income and other socioeconomic factors. In a regression model, race acts as a proxy variable, so apparent correlations with race are often caused, at least in part, by other factors. The next variable on the list is nbrnaliv, which indicates whether the pregnancy yielded multiple births. Twins and triplets tend to be smaller than other babies, so if we know whether our hypothetical coworker is expecting twins, that would help. Next on the list is paydu, which indicates whether the respondent owns her home. It is one of several incomerelated variables that turn out to be predictive. In datasets like the NSFG, income and wealth are correlated with just about everything. In this example, income is related to diet, health, health care, and other factors likely to affect birth weight. Some of the other variables on the list are things that would not be known until later, like bfeedwks, the number of weeks the baby was breast fed. We can’t use these variables for prediction, but you might want to speculate on reasons bfeedwks might be correlated with birth weight. Sometimes you start with a theory and use data to test it. Other times you start with data and go looking for possible theories. The second approach, which this section demonstrates, is called data mining. An advantage of data mining is that it can discover unexpected patterns. A hazard is that many of the patterns it discovers are either random or spurious. Having identified potential explanatory variables, I tested a few models and settled on this one: formula = ('totalwgt_lb ~ agepreg + C(race) + babysex==1 + '
'nbrnaliv>1 + paydu==1 + totincr')
results = smf.ols(formula, data=join).fit()
This formula uses some syntax we have not seen yet: C(race) tells the formula parser (Patsy) to treat race as a categorical variable, even though it is encoded numerically. The encoding for babysex is 1 for male, 2 for female; writing babysex==1 converts it to boolean, True for male and false for female. Similarly nbrnaliv>1 is True for multiple births and paydu==1 is True for respondents who own their houses. totincr is encoded numerically from 114, with each increment representing about $5000 in annual income. So we can treat these values as numerical, expressed in units of $5000. Here are the results of the model: Intercept 6.63 (0)
C(race)[T.2] 0.357 (5.43e29)
C(race)[T.3] 0.266 (2.33e07)
babysex == 1[T.True] 0.295 (5.39e29)
nbrnaliv > 1[T.True] 1.38 (5.1e37)
paydu == 1[T.True] 0.12 (0.000114)
agepreg 0.00741 (0.0035)
totincr 0.0122 (0.00188)
The estimated parameters for race are larger than I expected, especially since we control for income. The encoding is 1 for black, 2 for white, and 3 for other. Babies of black mothers are lighter than babies of other races by 0.27–0.36 lbs. As we’ve already seen, boys are heavier by about 0.3 lbs; twins and other multiplets are lighter by 1.4 lbs. People who own their homes have heavier babies by about 0.12 lbs, even when we control for income. The parameter for mother’s age is smaller than what we saw in Section 11.2, which suggests that some of the other variables are correlated with age, probably including paydu and totincr. All of these variables are statistically significant, some with very low pvalues, but R^{2} is only 0.06, still quite small. RMSE without using the model is 1.27 lbs; with the model it drops to 1.23. So your chance of winning the pool is not substantially improved. Sorry! 11.6 Logistic regressionIn the previous examples, some of the explanatory variables were numerical and some categorical (including boolean). But the dependent variable was always numerical. Linear regression can be generalized to handle other kinds of dependent variables. If the dependent variable is boolean, the generalized model is called logistic regression. If the dependent variable is an integer count, it’s called Poisson regression. As an example of logistic regression, let’s consider a variation on the office pool scenario. Suppose a friend of yours is pregnant and you want to predict whether the baby is a boy or a girl. You could use data from the NSFG to find factors that affect the “sex ratio”, which is conventionally defined to be the probability of having a boy. If you encode the dependent variable numerically, for example 0 for a girl and 1 for a boy, you could apply ordinary least squares, but there would be problems. The linear model might be something like this:
Where y is the dependent variable, and x_{1} and x_{2} are explanatory variables. Then we could find the parameters that minimize the residuals. The problem with this approach is that it produces predictions that are hard to interpret. Given estimated parameters and values for x_{1} and x_{2}, the model might predict y=0.5, but the only meaningful values of y are 0 and 1. It is tempting to interpret a result like that as a probability; for example, we might say that a respondent with particular values of x_{1} and x_{2} has a 50% chance of having a boy. But it is also possible for this model to predict y=1.1 or y=−0.1, and those are not valid probabilities. Logistic regression avoids this problem by expressing predictions in terms of odds rather than probabilities. If you are not familiar with odds, “odds in favor” of an event is the ratio of the probability it will occur to the probability that it will not. So if I think my team has a 75% chance of winning, I would say that the odds in their favor are three to one, because the chance of winning is three times the chance of losing. Odds and probabilities are different representations of the same information. Given a probability, you can compute the odds like this: o = p / (1p)
Given odds in favor, you can convert to probability like this: p = o / (o+1)
Logistic regression is based on the following model:
Where o is the odds in favor of a particular outcome; in the example, o would be the odds of having a boy. Suppose we have estimated the parameters β_{0}, β_{1}, and β_{2} (I’ll explain how in a minute). And suppose we are given values for x_{1} and x_{2}. We can compute the predicted value of logo, and then convert to a probability: o = np.exp(log_o)
p = o / (o+1)
So in the office pool scenario we could compute the predictive probability of having a boy. But how do we estimate the parameters? 11.7 Estimating parametersUnlike linear regression, logistic regression does not have a closed form solution, so it is solved by guessing an initial solution and improving it iteratively. The usual goal is to find the maximumlikelihood estimate (MLE), which is the set of parameters that maximizes the likelihood of the data. For example, suppose we have the following data: >>> y = np.array([0, 1, 0, 1])
>>> x1 = np.array([0, 0, 0, 1])
>>> x2 = np.array([0, 1, 1, 1])
And we start with the initial guesses β_{0}=−1.5, β_{1}=2.8, and β_{2}=1.1: >>> beta = [1.5, 2.8, 1.1]
Then for each row we can compute >>> log_o = beta[0] + beta[1] * x1 + beta[2] * x2
[1.5 0.4 0.4 2.4]
And convert from log odds to probabilities: >>> o = np.exp(log_o)
[ 0.223 0.670 0.670 11.02 ]
>>> p = o / (o+1)
[ 0.182 0.401 0.401 0.916 ]
Notice that when The likelihood of an outcome is p when y==1 and 1p when y==0. For example, if we think the probability of a boy is 0.8 and the outcome is a boy, the likelihood is 0.8; if the outcome is a girl, the likelihood is 0.2. We can compute that like this: >>> likes = y * p + (1y) * (1p)
[ 0.817 0.401 0.598 0.916 ]
The overall likelihood of the data is the product of likes: >>> like = np.prod(likes)
0.18
For these values of beta, the likelihood of the data is 0.18. The goal of logistic regression is to find parameters that maximize this likelihood. To do that, most statistics packages use an iterative solver like Newton’s method (see https://en.wikipedia.org/wiki/Logistic_regression#Model_fitting). 11.8 ImplementationStatsModels provides an implementation of logistic regression called logit, named for the function that converts from probability to log odds. To demonstrate its use, I’ll look for variables that affect the sex ratio. Again, I load the NSFG data and select pregnancies longer than 30 weeks: live, firsts, others = first.MakeFrames()
df = live[live.prglngth>30]
logit requires the dependent variable to be binary (rather than boolean), so I create a new column named boy, using astype(int) to convert to binary integers: df['boy'] = (df.babysex==1).astype(int)
Factors that have been found to affect sex ratio include parents’ age, birth order, race, and social status. We can use logistic regression to see if these effects appear in the NSFG data. I’ll start with the mother’s age: import statsmodels.formula.api as smf
model = smf.logit('boy ~ agepreg', data=df)
results = model.fit()
SummarizeResults(results)
logit takes the same arguments as ols, a formula in Patsy syntax and a DataFrame. The result is a Logit object that represents the model. It contains attributes called endog and exog that contain the endogenous variable, another name for the dependent variable, and the exogenous variables, another name for the explanatory variables. Since they are NumPy arrays, it is sometimes convenient to convert them to DataFrames: endog = pandas.DataFrame(model.endog, columns=[model.endog_names])
exog = pandas.DataFrame(model.exog, columns=model.exog_names)
The result of model.fit is a BinaryResults object, which is similar to the RegressionResults object we got from ols. Here is a summary of the results: Intercept 0.00579 (0.953)
agepreg 0.00105 (0.783)
R^2 6.144e06
The parameter of agepreg is positive, which suggests that older mothers are more likely to have boys, but the pvalue is 0.783, which means that the apparent effect could easily be due to chance. The coefficient of determination, R^{2}, does not apply to logistic regression, but there are several alternatives that are used as “pseudo R^{2} values.” These values can be useful for comparing models. For example, here’s a model that includes several factors believed to be associated with sex ratio: formula = 'boy ~ agepreg + hpagelb + birthord + C(race)'
model = smf.logit(formula, data=df)
results = model.fit()
Along with mother’s age, this model includes father’s age at birth (hpagelb), birth order (birthord), and race as a categorical variable. Here are the results: Intercept 0.0301 (0.772)
C(race)[T.2] 0.0224 (0.66)
C(race)[T.3] 0.000457 (0.996)
agepreg 0.00267 (0.629)
hpagelb 0.0047 (0.266)
birthord 0.00501 (0.821)
R^2 0.000144
None of the estimated parameters are statistically significant. The pseudoR^{2} value is a little higher, but that could be due to chance. 11.9 AccuracyIn the office pool scenario, we are most interested in the accuracy of the model: the number of successful predictions, compared with what we would expect by chance. In the NSFG data, there are more boys than girls, so the baseline strategy is to guess “boy” every time. The accuracy of this strategy is just the fraction of boys: actual = endog['boy']
baseline = actual.mean()
Since actual is encoded in binary integers, the mean is the fraction of boys, which is 0.507. Here’s how we compute the accuracy of the model: predict = (results.predict() >= 0.5)
true_pos = predict * actual
true_neg = (1  predict) * (1  actual)
results.predict returns a NumPy array of probabilities, which we
round off to 0 or 1. Multiplying by actual
yields 1 if we predict a boy and get it right, 0 otherwise. So,
Similarly, acc = (sum(true_pos) + sum(true_neg)) / len(actual)
The result is 0.512, slightly better than the baseline, 0.507. But, you should not take this result too seriously. We used the same data to build and test the model, so the model may not have predictive power on new data. Nevertheless, let’s use the model to make a prediction for the office pool. Suppose your friend is 35 years old and white, her husband is 39, and they are expecting their third child: columns = ['agepreg', 'hpagelb', 'birthord', 'race']
new = pandas.DataFrame([[35, 39, 3, 2]], columns=columns)
y = results.predict(new)
To invoke results.predict for a new case, you have to construct a DataFrame with a column for each variable in the model. The result in this case is 0.52, so you should guess “boy.” But if the model improves your chances of winning, the difference is very small. 11.10 ExercisesMy solution to these exercises is in Exercise 1
Suppose one of your coworkers is expecting a baby and you are
participating in an office pool to predict the date of birth.
Assuming that bets are placed during the 30th week of pregnancy, what
variables could you use to make the best prediction? You should limit
yourself to variables that are known before the birth, and likely to
be available to the people in the pool.
Exercise 2
The TriversWillard hypothesis suggests that for many mammals the
sex ratio depends on “maternal condition”; that is,
factors like the mother’s age, size, health, and social status.
See https://en.wikipedia.org/wiki/TriversWillard_hypothesis
Some studies have shown this effect among humans, but results are mixed. In this chapter we tested some variables related to these factors, but didn’t find any with a statistically significant effect on sex ratio. As an exercise, use a data mining approach to test the other variables in the pregnancy and respondent files. Can you find any factors with a substantial effect? Exercise 3
If the quantity you want to predict is a count, you can use Poisson
regression, which is implemented in StatsModels with a function called
poisson. It works the same way as ols and logit.
As an exercise, let’s use it to predict how many children a woman
has born; in the NSFG dataset, this variable is called numbabes.
Suppose you meet a woman who is 35 years old, black, and a college graduate whose annual household income exceeds $75,000. How many children would you predict she has born? Exercise 4
If the quantity you want to predict is categorical, you can use
multinomial logistic regression, which is implemented in StatsModels
with a function called mnlogit. As an exercise, let’s use it to
guess whether a woman is married, cohabitating, widowed, divorced,
separated, or never married; in the NSFG dataset, marital status is
encoded in a variable called rmarital.
Suppose you meet a woman who is 25 years old, white, and a high school graduate whose annual household income is about $45,000. What is the probability that she is married, cohabitating, etc? 11.11 Glossary

Are you using one of our books in a class?We'd like to know about it. Please consider filling out this short survey.
