Congressional Elections, Logistic Regression, and Feature Selection

In the spirit of our recent election season (hah), I’ve decided to do some exploration in using boosted trees to identify features that have the most predictive power. Using those selected features and some propreitary intuition, I’ll then use a logistic model to predict whether a congressional district will be democrat or republican.

For comparison’s sake, I’ll begin by specifying a logistic model with features that I pick without any empirical evidence to back up my decision.

The data describe all 435 districts in the 105th congress from 1997-1998. It contains demographic and employment data for each district as well as a variable indicating whether the district’srepresentative in the House was a Republican or not. Republican districts were coded as 1 while Democratic/independent districts were coded as 0. There are a total of 31 variables, as described below:

import pandas as pd
import numpy as np
import statsmodels.api as sm
from statsmodels.discrete.discrete_model import Logit
import xgboost as xgb
from xgboost.sklearn import XGBClassifier
from __future__ import division
import seaborn as sns
from sklearn.metrics import confusion_matrix, classification_report

%matplotlib inline
import matplotlib.pylab as plt
from matplotlib.pylab import rcParams
rcParams['figure.figsize'] = 18.5, 10.5

congress = pd.read_stata("http://rlhick.people.wm.edu/econ407/data/congressional_105.dta" )

For the naive logit model, I will estimate $repub = \beta_0 + \beta_1 per\_age65 + \beta_2 per\_black + \beta_3 per\_ bluecllr + \beta_4 city + \beta_5 mdnincm + \beta_6 per\_unemployed + \beta_7 union$

where per_ indicates that variable is a percentage of the district. I use percentages so that (presumably) the results are agnostic of the size of a district and can be generalized to any region in the US.

#making percentages
variables = ['unemplyd', 'age65', 'black', 'blucllr']

for v in variables:
    congress['per_' + v] = congress[[v, 'populatn']].apply(lambda row: (row[0] / row[1])*100, axis = 1)

Before I do anything, I check for null values and get a sense of the data

congress.isnull().sum()

state       0
fipstate    0
sc          0
cd          0
repub       1
age65       0
black       0
blucllr     0
city        0
coast       0
construc    0
cvllbrfr    0
enroll      0
farmer      0
finance     0
forborn     0
gvtwrkr     0
intrland    0
landsqmi    0
mdnincm     0
miltinst    0
miltmajr    0
miltpop     0
nucplant    0
popsqmi     0
populatn    0
rurlfarm    0
transprt    0
unemplyd    0
union       0
urban       0
whlretl     0
dtype: int64

After doing some research on Oklahoma’s first congressional district, they had republican congressman Steve Largent from 1994 to 2002. I will code that row as republican.

congress.repub.fillna(value = 1, inplace = True)

Lets see how balanced my data set is

sns.countplot(x = 'repub', data = congress, palette = 'muted');

print congress.repub.value_counts()

1.0    228
0.0    207
Name: repub, dtype: int64

So our data is fairly balanced, although republicans have the slight majority. Lets check out the relationship between the blue collar workers and black population:

sns.regplot(congress.per_black, congress.per_blucllr);

There doesn’t seem to be a definite relationship between the percentage of blue collar workers and percentage of black population. What about unemployment and the black population?

sns.regplot(congress.per_black, congress.per_unemplyd);

Those two variables are positively correlated.

Lets take a look at differences in race that elected a republican congressman:

sns.boxplot(x = 'repub', y = 'per_black', data = congress );
sns.swarmplot(x = 'repub', y = 'per_black', data = congress, color = ".15" )

<matplotlib.axes._subplots.AxesSubplot at 0x119de4a10>

facetgrid = sns.FacetGrid(congress,hue='repub',size=6)
facetgrid.map(sns.kdeplot,'per_black',shade=True)
facetgrid.set(xlim=(0, congress.per_black.max()))
facetgrid.add_legend()

<seaborn.axisgrid.FacetGrid at 0x119e327d0>

print "Districts that went Republican: \n",congress[congress.repub == 1].per_black.describe()

print "\nDistricts that went Democrat: \n",congress[congress.repub == 0].per_black.describe()

Districts that went Republican: 
count    228.000000
mean       7.139820
std        7.773183
min        0.178235
25%        2.148288
50%        4.678620
75%        8.204570
max       40.724916
Name: per_black, dtype: float64

Districts that went Democrat: 
count    207.000000
mean      17.020763
std       20.083805
min        0.131355
25%        2.776919
50%        8.122120
75%       21.369574
max       73.948710
Name: per_black, dtype: float64

The median value for districts that went Democrat is 4 % higher than those that went republican. In addition to the districts that went Republican are more closely centered around the mean and median, where as in the districts that went Democrat, there is a significant discrepancy between the mean and median. In fact, 25 % of the data falls between 8.12 % and 21.36 % for Democratic districts, but in Republican districts that number is 4.6 % and 8.2 %.

I expect that the percentage of the population that is black will have significant predicitive power in our model.

Examine the correlation between our variables of interest:

ind_vars = [ 'per_age65', 'per_black', 'per_blucllr', 'city','mdnincm', 'per_unemplyd', 'union']
congress[ind_vars].corr()

Nothing seriously outrageous to see here. Percent black and unemployment have the highest correlation at 54.9%, followed by city and unemployment at 29.9%. Something to keep in mind when we are doing feature selections.

quick aside about MLE and Logit

I’m not going to go into a ton detail about logistic regression, because it has been done ad nauseam. If you would like more info about it, the wiki is a great place to start.

I’d like to detail a bit about MLE or maximum likelihood estimation in the context of the logit model. Stats models uses maximum likelihood to estimate the coefficients. In order to use MLE, we assume that the errors in our model are independent and identically Logit distributed. The cumulative logit function is given by:

The log-likelihood is is then:

When I run a logitstic model in statsmodels, it estimates parameters $b$ such that the above sum is maximized.

So now that we know what’s going on behind the scenes, lets heckin do it :)

dep_var = 'repub'

x_const = sm.add_constant(congress[ind_vars])
y = congress[dep_var]

#Logit and probit models

logit_results = Logit(y, x_const).fit()
logit_results.summary()

Optimization terminated successfully.
         Current function value: 0.542084
         Iterations 7

logit_results.get_margeff(dummy = True).summary()

Before I look at how well the model is classifying, I’ll interpret the marginal effects:

For each variable that is a percentage, the marginal effect is the amount that the probability drops for each percent increase in that variable. Those include per_age65, per_black, per_bluecllr, per_unemplyd, and union. The most notable of these is that for every percentage increase in unemployment, the probability that the district elects a republican congressman or woman decreases by almost 27%!! Oh how the times have changed.

Another interesting result is that for every 10,000 $ increase in median income, the probability of a republican goes down by about 10.7 %. This goes against conventional wisdom that says that richer folks prefer republican canidates. Again, I think that if region was controlled for, the effect may be different.

Something that I think would be interesting to explore is how religious preferences within each district affect republican chances. Anyway, on to the classification reports

predictions = np.array([1 if x >= 0.5 else 0 for x in logit_results.predict()])

print classification_report(congress.repub, predictions)
print "Confusion matrix:\n"
print confusion_matrix(congress.repub, predictions)

             precision    recall  f1-score   support

        0.0       0.74      0.67      0.70       207
        1.0       0.72      0.78      0.75       228

avg / total       0.73      0.73      0.73       435

Confusion matrix:

[[139  68]
 [ 50 178]]

Okay, so the model isn’t awful, although it could use some serious improvement. Democratic districts are only classified correctly (139 / (139 + 68)) = 67% of the time. Okay so that’s not very good at all. The pseudo R2 is 0.2166, which is okay. A note about that:

McFadden’s (pseudo) R2 = (1 - (log likelihood full model / log likelihood only intercept)). The log likelihood of the intercept model is treated as a total sum of squares, and the log likelihood of the full model is treated as the sum of squared errors (like in approach 1). The ratio of the likelihoods suggests the level of improvement over the intercept model offered by the full model. A higher R2 suggests a better model that explains more of the variance.

We proceed boldy forward by diving in to Feature Selection

Using a Boosted classifier to do Feature Selection

A cool feature of the xgboost library is that after it estimates all the trees, it calculates for each tree the amount that each attribute split point improves the performance measure, weighted by the number of observations the node is responsible for. The performance measure may be the purity (Gini index) used to select the split points or another more specific error function. For more information, check out this stack overflow question.

I start with using all variables except those that are identifying information like state, FIPS code etc. I also make the rest of the population count variables in to percentages so that they are comparable to other districts.

needs_convert = ['construc', 'cvllbrfr', 'enroll', 'farmer', 'finance', 'forborn', 'gvtwrkr', 'miltpop',
                 'rurlfarm', 'transprt', 'urban', 'whlretl']
for v in needs_convert:
    congress['per_' + v] = congress[[v, 'populatn']].apply(lambda row: (row[0] / row[1])*100, axis = 1)

all_ind_vars = ['per_age65', 'per_black', 'per_blucllr', 'city','coast', 'per_construc', 'per_cvllbrfr', 
                'per_enroll','per_farmer', 'per_finance', 'per_forborn', 'per_gvtwrkr', 'intrland', 'landsqmi', 
                'mdnincm', 'miltinst', 'miltmajr', 'per_miltpop', 'nucplant', 'popsqmi','populatn', 'per_rurlfarm', 
                'per_transprt', 'per_unemplyd', 'union', 'per_urban','per_whlretl']

boosted = XGBClassifier(max_depth= 4, objective= 'binary:logistic')
boosted.fit(congress[all_ind_vars], congress.repub)

XGBClassifier(base_score=0.5, colsample_bylevel=1, colsample_bytree=1,
       gamma=0, learning_rate=0.1, max_delta_step=0, max_depth=4,
       min_child_weight=1, missing=None, n_estimators=100, nthread=-1,
       objective='binary:logistic', reg_alpha=0, reg_lambda=1,
       scale_pos_weight=1, seed=0, silent=True, subsample=1)

xgb.plot_importance(boosted.booster());

Some surprising results, some expected results.

Unexpected:

The percentage of the population that works construction is the most important feature… interesting. Personally, I would consider construction work to be blue collar work, however, they weren’t counted as part of the population that constituted blue collar workers. Maybe instead of campaigning with the steel workers, congressmen should pander to construction unions (do those even exist?). The size of the district also plays a significant role. It will be interesting to see whether size contributes positively or negatively to the probability of a republican.

Government, transport and utility workers also contribute significantly. A common theme here is features that are most important are either directly related to government work or somehow often implicated with government. Indeed, whether a person is employed or not seriously colors their view of how the government should be run.

Expected

I was right about unemployment and the black population playing a large role in elections! Sweet validation.

New Model

The new model I will estimate is:

I include all the top four variables that were selected for by the boosted trees. All of those features, with the exception of black, are indicators of employment. To control for age, race, and economic status of the district I include percentage black, over age 65, and median income. I include the amount of National Park land because that land is owned by the government, and provides employment and tourist attractions for citizens. City and coast are to control for urban populations and location respectively. Coast especially should account for the democratic “coastal elite”. Enroll is to control for those voters whose children are enrolled in public schools, also a factor that would color their view of the government.

Lets let er rip

new_vars = ['per_construc', 'per_black', 'per_unemplyd', 'per_transprt', 'per_gvtwrkr', 'per_enroll', 'per_age65',
            'mdnincm','populatn',  'intrland', 'city', 'coast' ]
new_x_const = sm.add_constant(congress[new_vars])

new_model = Logit(congress.repub, new_x_const).fit()
new_model.summary()

Optimization terminated successfully.
         Current function value: 0.504768
         Iterations 7

new_predictions = np.array([1 if x >= 0.5 else 0 for x in new_model.predict()])

print classification_report(congress.repub, new_predictions)
print "Confusion matrix:\n"
print confusion_matrix(congress.repub, new_predictions)

             precision    recall  f1-score   support

        0.0       0.76      0.68      0.72       207
        1.0       0.73      0.81      0.77       228

avg / total       0.75      0.74      0.74       435

Confusion matrix:

[[140  67]
 [ 44 184]]

Our pseudo R2 improves significantly from 0.2166 to 0.2705. The precision and recall also improve slightly, although they are still not what I would like them to be. Lets look at the marginal effects:

new_model.get_margeff(dummy=True).summary()

Percent unemployment remains the biggest factor for probability of republican, followed by population and then percentage of the population that is black. Interestingly, a larger population benefits republican canidates. So does having children enrolled in public schools.

What would happen if we controlled for region of the country with a series of dummy variables?

# these are regions as defined by the bureau of economic analysis
regions = {'new_england':['CT', 'ME', 'MA', 'NH', 'RI', 'VT'] , 
           'mid_east': ['DE', 'DC', 'MD', 'NJ', 'NY', 'PA'],
           'great_lakes': ['IL', 'IN', 'MI', 'OH', 'WI'], 
           'plains': ['IA', 'KS', 'MN', 'MO', 'NE', 'ND', 'SD'], 
           'south_east':['AL', 'AR', 'FL', 'GA', 'KY', 'LA', 'MS', 'NC', 'SC', 'TN', 'VA', 'WV'], 
           'south_west': ['AZ', 'NM', 'OK', 'TX'], 
           'rocky_mountain': ['CO', 'ID', 'MT', 'UT', 'WY'], 
           'far_west': ['CA', 'AK', 'HI', 'NV', 'OR', 'WA']}


#assigns 1 to that column if state is in that region
for key,val in regions.iteritems():
    congress['in_' + key] = congress.state.apply(lambda state: 1 if state in regions[key] else 0)

controlled_region = ['per_construc', 'per_black', 'per_unemplyd', 'per_transprt', 'per_gvtwrkr', 'per_enroll', 
                     'per_age65', 'mdnincm','populatn',  'intrland', 'city', 'coast', 'in_rocky_mountain', 
                     'in_plains', 'in_new_england', 'in_great_lakes', 'in_mid_east', 'in_south_west', 'in_south_east',
                     'in_far_west']

controlled_model = Logit(congress.repub, sm.add_constant(congress[controlled_region])).fit()
controlled_model.summary()

Optimization terminated successfully.
         Current function value: 0.482721
         Iterations 7


/Users/benjamindykstra/anaconda/lib/python2.7/site-packages/statsmodels/base/model.py:971: RuntimeWarning: invalid value encountered in sqrt
  return np.sqrt(np.diag(self.cov_params()))
/Users/benjamindykstra/anaconda/lib/python2.7/site-packages/scipy/stats/_distn_infrastructure.py:875: RuntimeWarning: invalid value encountered in greater
  return (self.a < x) & (x < self.b)
/Users/benjamindykstra/anaconda/lib/python2.7/site-packages/scipy/stats/_distn_infrastructure.py:875: RuntimeWarning: invalid value encountered in less
  return (self.a < x) & (x < self.b)
/Users/benjamindykstra/anaconda/lib/python2.7/site-packages/scipy/stats/_distn_infrastructure.py:1814: RuntimeWarning: invalid value encountered in less_equal
  cond2 = cond0 & (x <= self.a)

There seems to be a problem with the parameter estimates on the regions. Lets check the rank of our matrix and covariance matrix

np.linalg.matrix_rank(sm.add_constant(congress[controlled_region])), sm.add_constant(congress[controlled_region]).shape[1]

(20, 21)

controlled_model.cov_params()[['in_rocky_mountain','in_plains', 'in_new_england', 'in_great_lakes', 'in_mid_east', 
                               'in_south_west', 'in_south_east', 'in_far_west']]

So the matrix is not quite full rank and there is very little covariance between some of the dummy variable parameters. This is probably the cause of the problems.

Marginal effects:

controlled_model.get_margeff(dummy = True).summary()

/Users/benjamindykstra/anaconda/lib/python2.7/site-packages/statsmodels/discrete/discrete_margins.py:345: RuntimeWarning: invalid value encountered in sqrt
  return cov_me, np.sqrt(np.diag(cov_me))

new_predictions = np.array([1 if x >= 0.5 else 0 for x in controlled_model.predict()])

print classification_report(congress.repub, new_predictions)
print "Confusion matrix:\n"
print confusion_matrix(congress.repub, new_predictions)

             precision    recall  f1-score   support

        0.0       0.80      0.71      0.75       207
        1.0       0.76      0.83      0.80       228

avg / total       0.78      0.78      0.78       435

Confusion matrix:

[[148  59]
 [ 38 190]]

The model improves both true negatives and true positives. The R2 is also at a high of 0.3024 because of the better log likelihood value. It seems controlling for region of the country really helped out the model with its predictive accuracy. According to the marginal effects, a district in the south east has a 8 % chance higher of being republican, all else equal. Employment is still the most significant factore in the model and has the highest effect on probability as in the other models.

Using Sklearn’s Recursive Feature Elimination for Feature Selection

For comparison’s sake, I will use recursive feature elimination to select features from the dataset. From sklearn’s documentation:

Features whose absolute weights are the smallest are pruned from the current set features. That procedure is recursively repeated on the pruned set until the desired number of features to select is eventually reached.

In order to use this method, I will need to use sklearn’s logistic regression instead of statsmodels. Something to be aware of is the regularization term ‘C’ in sklearn - it needs to be set to an extremely high value so that the parameter estimates are the same as statsmodels.

from sklearn.feature_selection import RFECV
from sklearn.linear_model import LogisticRegression

new_ind_vars = ['per_age65', 'per_black', 'per_blucllr', 'city', 'coast', 'per_construc','per_cvllbrfr',
                'per_enroll','per_farmer','per_finance','per_forborn','per_gvtwrkr','intrland','landsqmi',
                'mdnincm','miltinst','miltmajr','per_miltpop','nucplant','popsqmi','populatn','per_rurlfarm',
                'per_transprt','per_unemplyd','union','per_urban','per_whlretl','in_rocky_mountain', 'in_plains', 
                'in_new_england','in_great_lakes', 'in_mid_east', 'in_south_west', 'in_south_east','in_far_west']

x = sm.add_constant(congress[new_ind_vars])
y = congress.repub

#cross validated recursive feature elimination, with auc scoring
rfecv = RFECV(estimator= LogisticRegression(C = 1e10),scoring = 'roc_auc' )
x_new = rfecv.fit_transform(x, y)

print x_new.shape, x.shape


#get the selected features
good_features = []
for i in range(len(x.columns.values)):
    if rfecv.support_[i]:
        good_features.append(x.columns.values[i])

print "The selected features: \n", good_features

(435, 14) (435, 36)
The selected features: 
['const', 'per_age65', 'city', 'coast', 'per_construc', 'per_cvllbrfr', 'per_finance', 'per_gvtwrkr', 'miltinst', 'per_miltpop', 'per_transprt', 'per_unemplyd', 'in_rocky_mountain', 'in_new_england']

A lot of the same features that were previously included by xgboost were selected for. Thats rather encouraging. I’ll fit the statsmodels logit to the data using the new features now.

selected_model = Logit(congress.repub, x[good_features]).fit()
selected_model.summary()

Optimization terminated successfully.
         Current function value: 0.497313
         Iterations 7

selected_model.get_margeff(dummy = True).summary()

Wow, so in this model percent unemployment doesn’t have nearly the same maginitude effect that it does in the other models. However, I like what I see with regions - New England has always been very liberal, and that is reflected in the dummy variables ‘in_new_england’ and ‘coast’.

Interestingly, percentage black was not selected for.. I still believe that that is a very important feature and should be included. Lets add it in and reexamine the results.

good_features.append('per_black')
new_selected_model = Logit(congress.repub, x[good_features]).fit()
new_selected_model.summary()

Optimization terminated successfully.
         Current function value: 0.488322
         Iterations 7

new_selected_model.get_margeff().summary()

new_selected_model.bic, selected_model.bic, controlled_model.bic

(515.97027448372125, 517.71751551732223, 541.47425095794426)

The percentage of unemployment has been restored to it’s former prominence. Other important variables (statistically and magnitude) include coast, construction workers, and New England. New England has the largest magnitude with a 30% reduction in probability of being republican if the district is in the region. The new model also has the lowest BIC.

Lets look at the predictive power:

new_selected_predictions = np.array([1 if x >= 0.5 else 0 for x in new_selected_model.predict()])

print classification_report(congress.repub, new_selected_predictions)
print "Confusion matrix:\n"
print confusion_matrix(congress.repub, new_selected_predictions)

             precision    recall  f1-score   support

        0.0       0.75      0.70      0.72       207
        1.0       0.74      0.79      0.76       228

avg / total       0.74      0.74      0.74       435

Confusion matrix:

[[144  63]
 [ 49 179]]

These are actually worse predictions that the previous controlled model. We choose that model instead.

Conclusion

While there are certainly more powerful (prediction wise) algorithms out there, logistic regression is nice because it can be interpreted meaningfully and actionably. With the help of xgboost, I picked out some features that I wouldn’t have otherwise and was able to improve the logit model. After controlling for region, I achieved the best results, with in class predictions at 71% and 83% for democrats and republicans, respectively. With more detailed information on things like district religious preferences, age distributions, previous voting preferences etc, I’m sure that the model would make much better predictions.

Thanks for following along! :)