Python 因果推断（下）

import numpy as np
import pandas as pd
pd.set_option('precision', 3)

# Data from Oreopoulos (2011)
path = "https://github.com/causal-methods/Data/raw/master/" 
df = pd.read_stata(path + "oreopoulos.dta")
df.head(5) 

# Transform the variable of interest in % 
df["callback"] = 100*df["callback"] 

# Restrict data to April and August 2008
df0 = df[(df.fall_data == 0)] 

mean = df0.groupby('name_ethnicity').agg([np.mean, np.size])
mean["callback"] 

prop = pd.crosstab(index= df0['name_ethnicity'], columns=df0['female'], 
            values=df0['callback'], aggfunc='mean')
prop 

# Restrict data to February and September 2009
df2 = df[(df.fall_data == 2)] 

mean = df2.groupby('name_ethnicity').agg([np.mean, np.size])
mean["callback"] 

import plotly.express as px

y = mean["callback"].values[:, 0]
x = mean["callback"].index

fig = px.bar(df2, x, y, color = x,
             title="Callback Rate for Interview by Name Ethnicity",
             labels={ "y": "Callback Rate (%)",
                      "x": "Name Ethnicity",                 
                      "color": ""} )

fig.update_layout(font_size = 17)
fig.show() 

prop = pd.crosstab(index= df2['name_ethnicity'], columns=df2['female'], 
            values=df2['callback'], aggfunc='mean')
prop 

import plotly.graph_objects as go

ethnicity = prop.index
male = prop.values[:,0]
female = prop.values[:,1]

fig = go.Figure(data=[
         go.Bar(name='Male', x = ethnicity, y = male),
         go.Bar(name='Female', x = ethnicity, y = female) ])

fig.update_layout(barmode='group', font_size = 17,
      title = "Callback Rate for Interview by Gender",
      yaxis = dict(title='Callback Rate (%)'),
      xaxis = dict(title='Name Ethnicity') )

fig.show() 

control = ['additional_credential', 'ba_quality',
           'extracurricular_skills', 'language_skills',
           'certificate', 'ma', 'same_exp', 'exp_highquality',
           'skillspeaking', 'skillsocialper', 'skillwriting']

df2.groupby('female').agg([np.mean])[control] 

df2.groupby(['female', 'name_ethnicity']).agg([np.mean])[control] 

import statsmodels.formula.api as smf

# Sample Restriction based on name ethnicity
Canada = df2.name_ethnicity == "Canada"
Indian = df2[(Canada) | (df2.name_ethnicity == "Indian")]
Chinese = df2[(Canada) | (df2.name_ethnicity == "Chinese")]
Chn_Cdn = df2[(Canada) | (df2.name_ethnicity == "Chn-Cdn")]
Greek = df2[(Canada) | (df2.name_ethnicity == "Greek")]

sample = [Indian, Chinese, Chn_Cdn, Greek]

#  Run the simple model for each ethnicity
# and save the results
model1 = "callback ~ female + C(type)"

result1 = []
for data in sample:
   ols = smf.ols(model1, data).fit(cov_type='HC1')
   result1.append(ols) 

C:\Anaconda\envs\textbook\lib\site-packages\statsmodels\tools\_testing.py:19: FutureWarning:

pandas.util.testing is deprecated. Use the functions in the public API at pandas.testing instead. 

#  Library to print professional publication
# tables in Latex, HTML, etc.
!pip install stargazer 

Requirement already satisfied: stargazer in c:\anaconda\envs\textbook\lib\site-packages (0.0.5) 

WARNING: Error parsing requirements for numpy: [Errno 2] No such file or directory: 'c:\\anaconda\\envs\\textbook\\lib\\site-packages\\numpy-1.19.2.dist-info\\METADATA' 

# Settings for a nice table
from stargazer.stargazer import Stargazer
stargazer = Stargazer(result1)

stargazer.title('Table 1 - Callback Rates by Resume Type')

names = ['Indian', 'Chinese', 'Chn_Cdn', 'Greek']
stargazer.custom_columns(names, [1, 1, 1, 1])

order = ['female', 'Intercept', 'C(type)[T.1.0]',
         'C(type)[T.2.0]', 'C(type)[T.3.0]',
         'C(type)[T.4.0]']     
stargazer.covariate_order(order)

dict1 = {'C(type)[T.1.0]': '1) Foreign Name, Cdn Educ and Exp',
         'C(type)[T.2.0]': '2) Foreign Name and Educ, Cdn exp',
         'C(type)[T.3.0]': '3) Foreign Name and Educ, Mixed Exp',
         'C(type)[T.4.0]': '4) All Foreign (Name, Educ, and Exp)',
              'Intercept': '0) English Name, Cdn Educ and Exp',
                 'female': 'Female'}
stargazer.rename_covariates(dict1)

stargazer 

model2 = "callback ~ female*C(type)"

result2 = []
for data in sample:
   ols = smf.ols(model2, data).fit(cov_type='HC1')
   result2.append(ols) 

stargazer = Stargazer(result2)

stargazer.title('Table 2 - Callback Rates by Resume Type and Gender')

stargazer.custom_columns(names, [1, 1, 1, 1])

dict2 = {'female:C(type)[T.1.0]':'[Female]x[1]',
         'female:C(type)[T.2.0]':'[Female]x[2]',
         'female:C(type)[T.3.0]':'[Female]x[3]',
         'female:C(type)[T.4.0]':'[Female]x[4]'}

list2 = list(dict2.keys())

dict2.update(dict1)
stargazer.rename_covariates(dict2)

list2 = order + list2
stargazer.covariate_order(list2)

stargazer 

control1 = "+ ba_quality + extracurricular_skills + language_skills"
control2 = "+ ma + exp_highquality"
model3 = "callback ~ female*C(type)" + control1 + control2

result3 = []
for data in sample:
   ols = smf.ols(model3, data).fit(cov_type='HC1')
   result3.append(ols) 

stargazer = Stargazer(result3)

stargazer.title('Table 3 - Callback Rates and Robustness Checks')

stargazer.custom_columns(names, [1, 1, 1, 1])

dict3 = {'ba_quality':'Top 200 world ranking university',
         'exp_highquality':'High quality work experience',
         'extracurricular_skills'	:'List extra-curricular activities',
         'language_skills':'Fluent in French and other languages',
         'ma':'Canadian master’s degree'}

list3 = list(dict3.keys())

dict3.update(dict2)
stargazer.rename_covariates(dict3)

list3 = list2 + list3
stargazer.covariate_order(list3)

stargazer 

pd.crosstab(index= [df2['name_ethnicity'], df2['female'],
                           df2['name']], columns=df2['type'], 
                         values=df2['callback'], aggfunc='mean') 

pd.crosstab(index= df2['occupation_type'],
                   columns=[df2['name_ethnicity'], df2['female']], 
                   values=df2['callback'], aggfunc='mean') 

FE = "+ C(occupation_type) + C(city) + C(name)"
model4 = "callback ~ female*C(type) " + control1 + control2 + FE

result4 = []
for data in sample:
   ols = smf.ols(model4, data).fit(cov_type='HC1')
   result4.append(ols) 

stargazer = Stargazer(result4)

stargazer.title('Table 4 - Callback Rates and Fixed Effects')
stargazer.custom_columns(names, [1, 1, 1, 1])
stargazer.rename_covariates(dict3)
stargazer.covariate_order(list3)

stargazer.add_line('Fixed Effects', ['', '', '', ''])
stargazer.add_line('Occupation', ['Yes', 'Yes', 'Yes', 'Yes'])
stargazer.add_line('Name', ['Yes', 'Yes', 'Yes', 'Yes'])
stargazer.add_line('City', ['Yes', 'Yes', 'Yes', 'Yes'])

stargazer 

C:\Anaconda\envs\textbook\lib\site-packages\statsmodels\base\model.py:1752: ValueWarning:

covariance of constraints does not have full rank. The number of constraints is 43, but rank is 1

C:\Anaconda\envs\textbook\lib\site-packages\statsmodels\base\model.py:1752: ValueWarning:

covariance of constraints does not have full rank. The number of constraints is 41, but rank is 39

C:\Anaconda\envs\textbook\lib\site-packages\statsmodels\base\model.py:1752: ValueWarning:

covariance of constraints does not have full rank. The number of constraints is 37, but rank is 35

C:\Anaconda\envs\textbook\lib\site-packages\statsmodels\base\model.py:1752: ValueWarning:

covariance of constraints does not have full rank. The number of constraints is 31, but rank is 29 

import numpy as np
import pandas as pd
pd.set_option('precision', 3)

# Data from Bisschop et al. (2017)
path = "https://github.com/causal-methods/Data/raw/master/" 
df = pd.read_stata(path + "CBSregist2015.dta")
df.head(5) 

big_cities = ["Amsterdam", "Rotterdam", "Den Haag"]
medium_cities = ["Utrecht", "Nijmegen", "Groningen",
                 "Heerlen", "Eindhoven", "Arnhem"]

# Classify cities
def classify(var):
    if var in big_cities:
        return "Big Cities"
    elif var in medium_cities:
        return "Medium Cities"    
    else:   
        return "No tippelzone"

df['group'] = df["city"].apply(classify) 

outcome = ['sexassaultpcN', 'rapepcN', 
	'drugspcN', 'maltreatpcN', 'weaponspcN']

df.groupby('group')[outcome].mean().T 

demographics = ['popul_100', 'pop_dens_100', 'popmale1565_100',
            'inkhh', 'educhpc', 'nondutchpc', 'insurWWAO_pc']

df.groupby('group')[demographics].mean().T 

political_party = ['mayorSoc', 'mayorLib', 'mayorChr']
df.groupby('group')[political_party].mean().T 

df['year'] = pd.DatetimeIndex(df['year']).year
df['Dyear'] = pd.Categorical(df.year)

# Set Panel Data
# Set city as the unit of analysis
df25 = df.set_index(['city1', 'year']) 

import statsmodels.formula.api as smf

Ys = ["lnsexassaultN", "lnrapeN", "lndrugsN",
      "lnweaponsN", "lnmaltreatN"]

base = "~ 1 + opening"
fe = "+ C(city) + C(Dyear)"

controls = ['logpopmale1565', 'logpopdens', 'inkhh', 
	'educhpc', 'nondutchpc', 'insurWWAO', 'mayorCDA',
  'mayorCU', 'mayorD66', 'mayorVVD']

Xs = ""
for var in controls:
    Xs = Xs + '+' + var

columns = []
for Y in Ys:
  result = smf.ols(Y + base + fe + Xs, df25).fit(cov_type='cluster',
                cov_kwds={'groups': df25['city']})
  columns.append(result) 

C:\Anaconda\envs\textbook\lib\site-packages\statsmodels\tools\_testing.py:19: FutureWarning: pandas.util.testing is deprecated. Use the functions in the public API at pandas.testing instead.
  import pandas.util.testing as tm 

#  Library to print professional publication
# tables in Latex, HTML, etc.
!pip install stargazer 

Requirement already satisfied: stargazer in c:\anaconda\envs\textbook\lib\site-packages (0.0.5) 

WARNING: Error parsing requirements for numpy: [Errno 2] No such file or directory: 'c:\\anaconda\\envs\\textbook\\lib\\site-packages\\numpy-1.19.2.dist-info\\METADATA' 

# Settings for a nice table
from stargazer.stargazer import Stargazer
stargazer = Stargazer(columns)

stargazer.title('The Impact of Tippelzone on Crime')

names = ['Sex Abuse', 'Rape', 'Drugs', 'Weapons', 'Assault']
stargazer.custom_columns(names, [1, 1, 1, 1, 1])

stargazer.covariate_order(['opening'])

stargazer.add_line('Covariates', ['Yes', 'Yes', 'Yes', 'Yes', 'Yes'])

stargazer.add_line('City Fixed Effects', ['Yes', 'Yes', 'Yes', 'Yes', 'Yes'])
stargazer.add_line('Year Fixed Effects', ['Yes', 'Yes', 'Yes', 'Yes', 'Yes'])

stargazer 

C:\Anaconda\envs\textbook\lib\site-packages\statsmodels\base\model.py:1752: ValueWarning: covariance of constraints does not have full rank. The number of constraints is 52, but rank is 24
  'rank is %d' % (J, J_), ValueWarning)
C:\Anaconda\envs\textbook\lib\site-packages\statsmodels\base\model.py:1752: ValueWarning: covariance of constraints does not have full rank. The number of constraints is 52, but rank is 24
  'rank is %d' % (J, J_), ValueWarning)
C:\Anaconda\envs\textbook\lib\site-packages\statsmodels\base\model.py:1752: ValueWarning: covariance of constraints does not have full rank. The number of constraints is 52, but rank is 24
  'rank is %d' % (J, J_), ValueWarning)
C:\Anaconda\envs\textbook\lib\site-packages\statsmodels\base\model.py:1752: ValueWarning: covariance of constraints does not have full rank. The number of constraints is 52, but rank is 24
  'rank is %d' % (J, J_), ValueWarning)
C:\Anaconda\envs\textbook\lib\site-packages\statsmodels\base\model.py:1752: ValueWarning: covariance of constraints does not have full rank. The number of constraints is 52, but rank is 24
  'rank is %d' % (J, J_), ValueWarning) 

import math
math.exp(-0.302) - 1 

-0.2606619351104681 

import numpy as np
import pandas as pd
pd.set_option('precision', 3)

# Data from Edelman et al. (2017)
path = "https://github.com/causal-methods/Data/raw/master/" 
df = pd.read_csv(path + "Airbnb.csv")
df.head(5) 

# Data for bar chart
count = pd.crosstab(df["graph_bins"], df["guest_black"])

import plotly.graph_objects as go

node = ['Conditional No', 'Conditional Yes', 'No',
        'No Response', 'Yes']
fig = go.Figure(data=[
    go.Bar(name='Guest is white', x=node, y=count[0]),
    go.Bar(name='Guest is African American', x=node, y=count[1]) ])

fig.update_layout(barmode='group',
  title_text = 'Host Responses by Race',
  font=dict(size=18) )

fig.show() 

import statsmodels.api as sm

df['const'] = 1 

# Column 1
#  The default missing ='drop' of statsmodels doesn't apply
# to the cluster variable. Therefore, it is necessary to drop
# the missing values like below to get the clustered standard 
# errors.
df1 = df.dropna(subset=['yes', 'guest_black', 'name_by_city'])
reg1 = sm.OLS(df1['yes'], df1[['const', 'guest_black']])
res1 = reg1.fit(cov_type='cluster',
                cov_kwds={'groups': df1['name_by_city']})

# Column 2
vars2 = ['yes', 'guest_black', 'name_by_city', 
        'host_race_black', 'host_gender_M']
df2 = df.dropna(subset = vars2)
reg2 = sm.OLS(df2['yes'], df2[['const', 'guest_black',
                    'host_race_black', 'host_gender_M']])
res2 = reg2.fit(cov_type='cluster',
                cov_kwds={'groups': df2['name_by_city']})

# Column 3
vars3 = ['yes', 'guest_black', 'name_by_city', 
         'host_race_black', 'host_gender_M',
         'multiple_listings', 'shared_property',
         'ten_reviews', 'log_price']
df3 = df.dropna(subset = vars3)
reg3 = sm.OLS(df3['yes'], df3[['const', 'guest_black',
                    'host_race_black', 'host_gender_M',
                    'multiple_listings', 'shared_property',
                    'ten_reviews', 'log_price']])
res3 = reg3.fit(cov_type='cluster',
                cov_kwds={'groups': df3['name_by_city']})

columns =[res1, res2, res3] 

C:\Anaconda\envs\textbook\lib\site-packages\statsmodels\tools\_testing.py:19: FutureWarning:

pandas.util.testing is deprecated. Use the functions in the public API at pandas.testing instead. 

#  Library to print professional publication
# tables in Latex, HTML, etc.
!pip install stargazer 

Requirement already satisfied: stargazer in c:\anaconda\envs\textbook\lib\site-packages (0.0.5) 

WARNING: Error parsing requirements for numpy: [Errno 2] No such file or directory: 'c:\\anaconda\\envs\\textbook\\lib\\site-packages\\numpy-1.19.2.dist-info\\METADATA' 

# Settings for a nice table
from stargazer.stargazer import Stargazer
stargazer = Stargazer(columns)
stargazer.title('The Impact of Race on Likelihood of Acceptance')
stargazer 

control = ['host_race_white', 'host_race_black', 'host_gender_F', 
	'host_gender_M', 'price', 'bedrooms', 'bathrooms', 'number_of_reviews', 
	'multiple_listings', 'any_black', 'tract_listings', 'black_proportion']

df.describe()[control].T 

result = []

for var in control:
    # Do the T-test and save the p-value
    pvalue = sm.OLS(df[var], df[['const', 'guest_black']],
               missing = 'drop').fit().pvalues[1]
    result.append(pvalue) 

ttest = df.groupby('guest_black').agg([np.mean])[control].T
ttest['p-value'] = result
ttest 

female = df['guest_gender']=='female'
df[female].groupby(['guest_race', 'guest_first_name'])['yes'].mean() 

guest_race  guest_first_name
black       Lakisha             0.433
            Latonya             0.370
            Latoya              0.442
            Tamika              0.482
            Tanisha             0.413
white       Allison             0.500
            Anne                0.567
            Kristen             0.486
            Laurie              0.508
            Meredith            0.498
Name: yes, dtype: float64 

male = df['guest_gender']=='male'
df[male].groupby(['guest_race', 'guest_first_name'])['yes'].mean() 

guest_race  guest_first_name
black       Darnell             0.412
            Jamal               0.354
            Jermaine            0.379
            Kareem              0.436
            Leroy               0.371
            Rasheed             0.409
            Tyrone              0.377
white       Brad                0.419
            Brent               0.494
            Brett               0.466
            Greg                0.467
            Jay                 0.581
            Todd                0.448
Name: yes, dtype: float64 

pd.crosstab(index= [df['host_gender_F'], df['host_race']],
            columns=[df['guest_gender'], df['guest_race']], 
            values=df['yes'], aggfunc='mean')