<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Internet use and religion in Europe
This notebook presents a quick-and-dirty analysis of the association between Internet use and religion in Europe, using data from the European Social Survey (http
Step1:  The following function selects the columns I need.
Step2:  Read data from Cycle 1.
TODO
Step3:  Read data from Cycle 2.
Step4:  Read data from Cycle 3.
Step5:  Read data from Cycle 4.
Step6:  Read data from Cycle 5.
Step7:  Concatenate the cycles.
TODO
Step8:  TV watching time on average weekday
Step9:  Radio listening, total time on average weekday.
Step10:  Newspaper reading, total time on average weekday.
Step11:  TV watching
Step12:  Radio listening
Step13:  Newspaper reading
Step14:  Personal use of Internet, email, www
Step15:  Belong to a particular religion or denomination
Step16:  How religious
Step17:  Total household net income, all sources
TODO
Step18:  Shift income to mean near 0.
Step19:  Year born
Step20:  Shifted to mean near 0
Step21:  Number of years of education
Step22:  There are a bunch of really high values for eduyrs, need to investigate.
Step23:  Shift to mean near 0
Step24:  Country codes
Step25:  Make a binary dependent variable
Step26:  Run the model
Step27:  Here's the model with all control variables and all media variables
Step28:  Most of the media variables are not statistically significant.  If we drop the politial media variables, we get a cleaner model
Step29:  And if we fill missing values for income, cleaner still.
Step30:  Now all variables have small p-values.  All parameters have the expected signs
Step31:  Running a similar model with degree of religiosity.
Step32:  Group by country
Step33:  Run a sample country
Step34:  Run all countries
Step35:  In more than half of the countries, the association between Internet use and religious affiliation is statistically significant.  In all except two, the association is negative.
In many countries we've lost a substantial number of observations due to missing data.  Really need to fill that in!
<ASSISTANT_TASK:>
Python Code:
from __future__ import print_function, division
import numpy as np
import pandas as pd
import statsmodels.formula.api as smf
%matplotlib inline
Explanation: Internet use and religion in Europe
This notebook presents a quick-and-dirty analysis of the association between Internet use and religion in Europe, using data from the European Social Survey (http://www.europeansocialsurvey.org).
Copyright 2015 Allen Downey
MIT License: http://opensource.org/licenses/MIT
End of explanation
def select_cols(df):
    cols = ['cntry', 'tvtot', 'tvpol', 'rdtot', 'rdpol', 'nwsptot', 'nwsppol', 'netuse', 
            'rlgblg', 'rlgdgr', 'eduyrs', 'hinctnta', 'yrbrn', 'eisced', 'pspwght', 'pweight']
    df = df[cols]
    return df
Explanation: The following function selects the columns I need.
End of explanation
df1 = pd.read_stata('ESS1e06_4.dta', convert_categoricals=False)
df1['hinctnta'] = df1.hinctnt
df1 = select_cols(df1)
df1.head()
Explanation: Read data from Cycle 1.
TODO: investigate the difference between hinctnt and hinctnta; is there a recode that reconciles them?
End of explanation
df2 = pd.read_stata('ESS2e03_4.dta', convert_categoricals=False)
df2['hinctnta'] = df2.hinctnt
df2 = select_cols(df2)
df2.head()
Explanation: Read data from Cycle 2.
End of explanation
df3 = pd.read_stata('ESS3e03_5.dta', convert_categoricals=False)
df3['hinctnta'] = df3.hinctnt
df3 = select_cols(df3)
df3.head()
Explanation: Read data from Cycle 3.
End of explanation
df4 = pd.read_stata('ESS4e04_3.dta', convert_categoricals=False)
df4 = select_cols(df4)
df4.head()
Explanation: Read data from Cycle 4.
End of explanation
df5 = pd.read_stata('ESS5e03_2.dta', convert_categoricals=False)
df5 = select_cols(df5)
df5.head()
Explanation: Read data from Cycle 5.
End of explanation
df = pd.concat([df1, df2, df3, df4, df5], ignore_index=True)
df.head()
Explanation: Concatenate the cycles.
TODO: Have to resample each cycle before concatenating.
End of explanation
df.tvtot.replace([77, 88, 99], np.nan, inplace=True)
df.tvtot.value_counts().sort_index()
Explanation: TV watching time on average weekday
End of explanation
df.rdtot.replace([77, 88, 99], np.nan, inplace=True)
df.rdtot.value_counts().sort_index()
Explanation: Radio listening, total time on average weekday.
End of explanation
df.nwsptot.replace([77, 88, 99], np.nan, inplace=True)
df.nwsptot.value_counts().sort_index()
Explanation: Newspaper reading, total time on average weekday.
End of explanation
df.tvpol.replace([66, 77, 88, 99], np.nan, inplace=True)
df.tvpol.value_counts().sort_index()
Explanation: TV watching: news, politics, current affairs
End of explanation
df.rdpol.replace([66, 77, 88, 99], np.nan, inplace=True)
df.rdpol.value_counts().sort_index()
Explanation: Radio listening: news, politics, current affairs
End of explanation
df.nwsppol.replace([66, 77, 88, 99], np.nan, inplace=True)
df.nwsppol.value_counts().sort_index()
Explanation: Newspaper reading: politics, current affairs
End of explanation
df.netuse.replace([77, 88, 99], np.nan, inplace=True)
df.netuse.value_counts().sort_index()
Explanation: Personal use of Internet, email, www
End of explanation
df.rlgblg.replace([7, 8, 9], np.nan, inplace=True)
df.rlgblg.value_counts().sort_index()
Explanation: Belong to a particular religion or denomination
End of explanation
df.rlgdgr.replace([77, 88, 99], np.nan, inplace=True)
df.rlgdgr.value_counts().sort_index()
Explanation: How religious
End of explanation
df.hinctnta.replace([77, 88, 99], np.nan, inplace=True)
df.hinctnta.value_counts().sort_index()
Explanation: Total household net income, all sources
TODO: It looks like one cycle measured HINCTNT on a 12 point scale.  Might need to reconcile
End of explanation
df['hinctnta5'] = df.hinctnta - 5
df.hinctnta5.describe()
Explanation: Shift income to mean near 0.
End of explanation
df.yrbrn.replace([7777, 8888, 9999], np.nan, inplace=True)
df.yrbrn.describe()
Explanation: Year born
End of explanation
df['yrbrn60'] = df.yrbrn - 1960
df.yrbrn60.describe()
Explanation: Shifted to mean near 0
End of explanation
df.eduyrs.replace([77, 88, 99], np.nan, inplace=True)
df.loc[df.eduyrs > 25, 'eduyrs'] = 25
df.eduyrs.value_counts().sort_index()
Explanation: Number of years of education
End of explanation
df.eduyrs.describe()
Explanation: There are a bunch of really high values for eduyrs, need to investigate.
End of explanation
df['eduyrs12'] = df.eduyrs - 12
df.eduyrs12.describe()
Explanation: Shift to mean near 0
End of explanation
df.cntry.value_counts().sort_index()
Explanation: Country codes
End of explanation
df['hasrelig'] = (df.rlgblg==1).astype(int)
Explanation: Make a binary dependent variable
End of explanation
def run_model(df, formula):
    model = smf.logit(formula, data=df)    
    results = model.fit(disp=False)
    return results
Explanation: Run the model
End of explanation
formula = ('hasrelig ~ yrbrn60 + eduyrs12 + hinctnta5 +'
           'tvtot + tvpol + rdtot + rdpol + nwsptot + nwsppol + netuse')
res = run_model(df, formula)
res.summary()
Explanation: Here's the model with all control variables and all media variables:
End of explanation
formula = ('hasrelig ~ yrbrn60 + eduyrs12 + hinctnta5 +'
           'tvtot + rdtot + nwsptot + netuse')
res = run_model(df, formula)
res.summary()
Explanation: Most of the media variables are not statistically significant.  If we drop the politial media variables, we get a cleaner model:
End of explanation
def fill_var(df, var):
    fill = df[var].dropna().sample(len(df), replace=True)
    fill.index = df.index
    df[var].fillna(fill, inplace=True)
fill_var(df, var='hinctnta5')
formula = ('hasrelig ~ yrbrn60 + eduyrs12 + hinctnta5 +'
           'tvtot + rdtot + nwsptot + netuse')
res = run_model(df, formula)
res.summary()
Explanation: And if we fill missing values for income, cleaner still.
End of explanation
def extract_res(res, var='netuse'):
    param = res.params[var]
    pvalue = res.pvalues[var]
    stars = '**' if pvalue < 0.01 else '*' if pvalue < 0.05 else ''
    return res.nobs, param, stars
extract_res(res)
Explanation: Now all variables have small p-values.  All parameters have the expected signs:
Younger people are less affiliated.
More educated people are less affiliated.
Higher income people are less affiliated (although this could go either way)
Consumers of all media are less affiliated.
The strength of the Internet effect is stronger than for other media.
These results are consistent in each cycle of the data, and across a few changes I've made in the cleaning process.
However, these results should be considered preliminary:
I have not dealt with the stratification weights.
I have not dealt with missing data (particularly important for education)
Nevertheless, I'll run a breakdown by country.
Here's a function to extract the parameter associated with netuse:
End of explanation
formula = ('rlgdgr ~ yrbrn60 + eduyrs12 + hinctnta5 +'
           'tvtot + rdtot + nwsptot + netuse')
model = smf.ols(formula, data=df)    
res = model.fit(disp=False)
res.summary()
Explanation: Running a similar model with degree of religiosity.
End of explanation
grouped = df.groupby('cntry')
for name, group in grouped:
    print(name, len(group))
Explanation: Group by country:
End of explanation
gb = grouped.get_group('DK')
run_model(gb, formula).summary()
Explanation: Run a sample country
End of explanation
for name, group in grouped:
    try:
        fill_var(group, var='hinctnta5')
        res = run_model(group, formula)
        nobs, param, stars = extract_res(res)
        arrow = '<--' if stars and param > 0 else ''
        print(name, len(group), nobs, '%0.3g'%param, stars, arrow, sep='\t')
    except:
        print(name, len(group), ' ', 'NA', sep='\t')
Explanation: Run all countries
End of explanation
group = grouped.get_group('FR')
len(group)
for col in group.columns:
    print(col, sum(group[col].isnull()))
fill_var(group, 'hinctnta5')
formula
res = run_model(group, formula)
res.summary()
Explanation: In more than half of the countries, the association between Internet use and religious affiliation is statistically significant.  In all except two, the association is negative.
In many countries we've lost a substantial number of observations due to missing data.  Really need to fill that in!
End of explanation
<END_TASK>