<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Will Millennials Ever Get Married?
This notebook contains analysis related to a paper on marriage patterns in the U.S., based on data from the National Survey of Family Growth (NSFG).
It is based on Chapter 13 of Think Stats, 2nd Edition, by Allen Downey, available from thinkstats2.com
Step1:  Make a table showing the number of respondents in each cycle
Step2:  Check for missing values in agemarry
Step3:  Estimate the hazard function for the 80s cohort (curious to see what's going on during the "marriage strike")
Step4:  Run the same analysis for the 70s cohort (to extract $\lambda(33)$).
Step5:  Use the 30s cohort to demonstrate the simple way to do survival analysis, by computing the survival function directly.
Step6:  Then use the SurvivalFunction to compute the HazardFunction
Step7:  Make the first figure, showing sf and hf for the 30s cohort
Step8:  Make some pivot tables, just to see where the data are
Step9:  The following pivot table is not as helpful as it could be, since it doesn't show the number at risk.
Step10:  Estimate the survival curve for each cohort
Step11:  Check a sample
Step13:  Make the figure showing estimated survival curves
Step15:  Make a table of marriage rates for each cohort at each age
Step16:  Generate projections
Step17:  And make the figure showing projections
Step18:  Make the table again with the projections filled in.
Step19:  Marriage patterns for men
<ASSISTANT_TASK:>
Python Code:
%matplotlib inline
import pandas as pd
import numpy as np
import seaborn as sns
import math
import matplotlib.pyplot as plt
from matplotlib import pylab
from scipy.interpolate import interp1d
from scipy.misc import derivative
import thinkstats2
import thinkplot
from thinkstats2 import Cdf
import survival
import marriage
%time df = pd.read_hdf('FemMarriageData.hdf', 'FemMarriageData')
df.shape
Explanation: Will Millennials Ever Get Married?
This notebook contains analysis related to a paper on marriage patterns in the U.S., based on data from the National Survey of Family Growth (NSFG).
It is based on Chapter 13 of Think Stats, 2nd Edition, by Allen Downey, available from thinkstats2.com
End of explanation
df.cycle.value_counts().sort_index()
def format_date_range(array):
    a, b = array.astype(int)
    return '%d--%d' % (a, b)
def SummarizeCycle(cycle, df):
    ages = df.age.min(), df.age.max()
    ages= np.array(ages)
    
    intvws = df.cmintvw.min(), df.cmintvw.max()
    intvws = np.array(intvws) / 12 + 1900
    
    births = df.cmbirth.min(), df.cmbirth.max()
    births = np.array(births) / 12 + 1900
    
    intvw_dates = format_date_range(intvws)
    birth_dates = format_date_range(births)
    
    print(cycle, ' & ', intvw_dates, '&', len(df), '&', birth_dates, r'\\')
for cycle, group in df.groupby('cycle'):
    SummarizeCycle(cycle, group)
Explanation: Make a table showing the number of respondents in each cycle:
End of explanation
marriage.DigitizeResp(df)
grouped = df.groupby('birth_index')
for name, group in iter(grouped):
    age_range = '%d--%d' % (int(group.age.min()), int(group.age_index.max()))
    print(name, '&', len(group), '&', age_range,
                '&', len(group[group.evrmarry]), '&', sum(group.missing), r'\\')
def ComputeCutoffs(df):
    grouped = df.groupby('birth_index')
    cutoffs = {}
    for name, group in sorted(grouped):
        cutoffs[name] = int(group.age.max())
    return cutoffs
cutoffs = ComputeCutoffs(df)
cutoffs
Explanation: Check for missing values in agemarry:
def CheckAgeVars(df):
    print(sum(df[df.evrmarry].agemarry.isnull()))
for cycle, group in df.groupby('cycle'):
    CheckAgeVars(group)
Generate a table with the number of respondents in each cohort:
End of explanation
cohort = grouped.get_group(80)
missing = (cohort.evrmarry & cohort.agemarry.isnull())
cohort = cohort[~missing]
complete = cohort[cohort.evrmarry].agemarry_index
ongoing = cohort[~cohort.evrmarry].age_index
hf = survival.EstimateHazardFunction(complete, ongoing, verbose=True)
Explanation: Estimate the hazard function for the 80s cohort (curious to see what's going on during the "marriage strike")
End of explanation
cohort = grouped.get_group(70)
missing = (cohort.evrmarry & cohort.agemarry.isnull())
cohort = cohort[~missing]
complete = cohort[cohort.evrmarry].agemarry_index
ongoing = cohort[~cohort.evrmarry].age_index
hf = survival.EstimateHazardFunction(complete, ongoing, verbose=True)
Explanation: Run the same analysis for the 70s cohort (to extract $\lambda(33)$).
End of explanation
cohort = grouped.get_group(30)
sf = survival.MakeSurvivalFromSeq(cohort.agemarry_index.fillna(np.inf))
ts, ss = sf.Render()
print(ss)
thinkplot.Plot(ts, ss)
thinkplot.Config(xlim=[12, 42])
Explanation: Use the 30s cohort to demonstrate the simple way to do survival analysis, by computing the survival function directly.
End of explanation
hf = sf.MakeHazardFunction()
ts, lams = hf.Render()
print(lams)
thinkplot.Plot(ts, lams)
thinkplot.Config(xlim=[12, 42])
Explanation: Then use the SurvivalFunction to compute the HazardFunction:
End of explanation
options = dict(formats=['pdf', 'png'], clf=False)
thinkplot.PrePlot(rows=2)
thinkplot.Plot(sf, label='survival')
thinkplot.Config(xlim=[13, 41], ylim=[0, 1.05])
plt.ylabel('Survival Function')
thinkplot.SubPlot(2)
thinkplot.Plot(hf, label='hazard')
thinkplot.Config(xlabel='age(years)', ylabel='Hazard function', xlim=[13, 41])
plt.ylabel('Hazard Function')
plt.xlabel('Age (years)')
thinkplot.Save(root='figs/marriage1', **options)
thinkplot.Plot(sf, label='30s')
thinkplot.Config(xlim=[13, 41], ylim=[0, 1.05])
plt.xlabel('Age (years)', fontsize=14)
plt.ylabel('Survival function', fontsize=14)
thinkplot.Save(root='figs/marriage2', **options)
thinkplot.Plot(hf, label='30s')
thinkplot.Config(xlim=[13, 41])
plt.xlabel('Age (years)', fontsize=14)
plt.ylabel('Hazard function', fontsize=14)
thinkplot.Save(root='figs/marriage3', **options)
Explanation: Make the first figure, showing sf and hf for the 30s cohort:
End of explanation
pt = df.pivot_table(index='birth_index', columns='age_index', values='age', aggfunc=len, fill_value=0)
pt
Explanation: Make some pivot tables, just to see where the data are:
End of explanation
df.pivot_table(index='birth_index', columns='agemarry_index', values='age', aggfunc=len, fill_value=0)
Explanation: The following pivot table is not as helpful as it could be, since it doesn't show the number at risk.
End of explanation
df['complete'] = df.evrmarry
df['complete_var'] = df.agemarry_index
df['ongoing_var'] = df.age_index
df['complete_missing'] = df.complete & df.complete_var.isnull()
df['ongoing_missing'] = ~df.complete & df.ongoing_var.isnull()
# for some marriages, we don't have the date of marriage
for cycle, group in df.groupby('cycle'):
    print(cycle, sum(group.complete_missing), sum(group.ongoing_missing))
resps = [group for cycle, group in df.groupby('cycle')]
iters = 101
%time sf_map = marriage.EstimateSurvivalByCohort(resps, iters=iters, cutoffs=cutoffs)
del sf_map[30]
try: 
    del sf_map[100]
except KeyError:
    pass
Explanation: Estimate the survival curve for each cohort:
End of explanation
for sf in sf_map[90]:
    print(sf.ss)
    print(sf.Prob(34))
    break
for sf in sf_map[80]:
    print(sf.ss)
    print(sf.Prob(34))
    break
Explanation: Check a sample:
End of explanation
def PlotSurvivalFunctions(root, sf_map, sf_map_pred=None, **options):
    if sf_map_pred:
        marriage.PlotSurvivalFunctions(sf_map_pred, predict_flag=True)
    marriage.PlotSurvivalFunctions(sf_map)
    thinkplot.config(xlabel='Age (years)',
                     ylabel='Fraction never married',
                     xlim=[13, 50], 
                     ylim=[0, 1.05], 
                     loc='upper right', 
                     frameon=False,
                     **options)
    plt.tight_layout()
    thinkplot.save(root=root, formats=['pdf', 'png'])
def set_palette(*args, **kwds):
    Set the matplotlib color cycler.
    
    args, kwds: same as for sns.color_palette
    
    Also takes a boolean kwd, `reverse`, to indicate
    whether the order of the palette should be reversed.
    
    returns: list of colors
    
    reverse = kwds.pop('reverse', False)
    palette = sns.color_palette(*args, **kwds)
    
    palette = list(palette)
    if reverse:
        palette.reverse()
        
    cycler = plt.cycler(color=palette)
    plt.gca().set_prop_cycle(cycler)
    return palette
def draw_age_lines(ages):
    for age in ages:
        plt.axvline(age, color='gray', linestyle='dotted', alpha=0.3)
palette = set_palette('hls', 6)
draw_age_lines(ages)
options_w = dict(title='Women in the U.S. by decade of birth')
PlotSurvivalFunctions('figs/marriage4', sf_map, None, **options_w)
Explanation: Make the figure showing estimated survival curves:
End of explanation
def MakeTable(sf_map, ages):
    t = []
    for name, sf_seq in sorted(sf_map.items()):
        ts, ss = marriage.MakeSurvivalCI(sf_seq, [50])
        ss = ss[0]
        vals = [np.interp(age, ts, ss, right=np.nan) for age in ages]
        t.append((name, vals))
    return t
def MakePercentageTable(sf_map, ages):
    Prints percentage unmarried for each cohort at each age.
    
    t = MakeTable(sf_map, ages)
    for name, sf_seq in sorted(sf_map.items()):
        ts, ss = marriage.MakeSurvivalCI(sf_seq, [50])
        ss = ss[0]
        vals = [np.interp(age, ts, ss, right=np.nan) for age in ages]
        print(name, '&', ' & '.join('%0.0f' % (val*100) for val in vals), r'\\')
        
ages = [24, 34, 44]
MakePercentageTable(sf_map, ages=ages)
Explanation: Make a table of marriage rates for each cohort at each age:
End of explanation
%time sf_map_pred = marriage.EstimateSurvivalByCohort(resps, iters=iters, cutoffs=cutoffs, predict_flag=True)
del sf_map_pred[30]
try: 
    del sf_map[100]
except KeyError:
    pass
for cohort, seq in sf_map_pred.items():
    if cohort > 90:
        break
    medians = [sf.MakeCdf().Value(0.5) for sf in seq]
    print(cohort, np.median(medians))
Explanation: Generate projections:
End of explanation
palette = set_palette('hls', 6)
draw_age_lines(ages)
PlotSurvivalFunctions('figs/marriage5', sf_map, sf_map_pred, **options_w)
Explanation: And make the figure showing projections:
End of explanation
MakePercentageTable(sf_map_pred, ages)
def PlotFractions(sf_map, ages, label_flag=False, **options):
    t = MakeTable(sf_map, ages)
    cohorts, cols = zip(*t)
    rows = zip(*cols)
    thinkplot.PrePlot(3)
    t = list(zip(ages, rows))
    for age, row in reversed(t):
        label = 'at age %d' % age if label_flag else ''
        thinkplot.Plot(cohorts, row, label=label, **options)
PlotFractions(sf_map_pred, ages, color='gray', linestyle='dashed', linewidth=2)
PlotFractions(sf_map, ages, label_flag=True, alpha=1)
fontsize=12
thinkplot.Text(36, 0.26, '24', fontsize=fontsize)
thinkplot.Text(37, 0.13, '9', fontsize=fontsize)
thinkplot.Text(37, 0.07, '7', fontsize=fontsize)
thinkplot.Text(90, 0.85, '80', fontsize=fontsize)
thinkplot.Text(90, 0.56, '51', fontsize=fontsize)
thinkplot.Text(89.5, 0.47, '42', fontsize=fontsize)
thinkplot.Text(80, 0.42, '35', fontsize=fontsize)
thinkplot.Text(70, 0.18, '18', fontsize=fontsize)
thinkplot.Config(xlim=[34, 97], ylim=[0, 1], legend=True, loc='upper left',
                 xlabel='cohort (decade)', ylabel='Fraction unmarried', 
                 title='Women in the U.S.')
thinkplot.Save(root='figs/marriage6', **options)
Explanation: Make the table again with the projections filled in.
End of explanation
%time df2 = pd.read_hdf('MaleMarriageData.hdf', 'MaleMarriageData')
df2.shape
for cycle, group in df2.groupby('cycle'):
    SummarizeCycle(cycle, group)
sum(df2.missing)
marriage.DigitizeResp(df2)
grouped = df2.groupby('birth_index')
for name, group in iter(grouped):
    age_range = '%d--%d' % (int(group.age.min()), int(group.age_index.max()))
    print(name, '&', len(group), '&', age_range,
                '&', len(group[group.evrmarry]), '&', sum(group.missing), r'\\')
cutoffs2 = ComputeCutoffs(df2)
cutoffs2
resps2 = [group for cycle, group in df2.groupby('cycle')]
%time sf_map_male = marriage.EstimateSurvivalByCohort(resps2, iters=iters, cutoffs=cutoffs2)
del sf_map_male[100]
palette = set_palette('hls', 6)
draw_age_lines(ages)
options_m = dict(title='Men in the U.S. by decade of birth')
PlotSurvivalFunctions('figs/marriage7', sf_map_male, None, **options_m)
%time sf_map_male_pred = marriage.EstimateSurvivalByCohort(resps2, iters=iters, cutoffs=cutoffs2, predict_flag=True)
del sf_map_male_pred[100]
for cohort, seq in sf_map_male_pred.items():
    if cohort > 90:
        break
    medians = [sf.MakeCdf().Value(0.5) for sf in seq]
    print(cohort, np.median(medians))
palette = set_palette('hls', 6)
draw_age_lines(ages)
PlotSurvivalFunctions('figs/marriage8', sf_map_male, sf_map_male_pred, **options_m)
MakePercentageTable(sf_map_male, ages)
MakePercentageTable(sf_map_male_pred, ages)
PlotFractions(sf_map_male_pred, ages, color='gray', linestyle='dashed', linewidth=2)
PlotFractions(sf_map_male, ages, label_flag=True, alpha=1)
fontsize=12
thinkplot.Text(46, 0.69, '68', fontsize=fontsize)
thinkplot.Text(46, 0.30, '26', fontsize=fontsize)
thinkplot.Text(46, 0.20, '18', fontsize=fontsize)
thinkplot.Text(70, 0.18, '19', fontsize=fontsize)
thinkplot.Text(80, 0.43, '43', fontsize=fontsize)
thinkplot.Text(90, 0.89, '86', fontsize=fontsize)
thinkplot.Text(90, 0.56, '52', fontsize=fontsize)
thinkplot.Text(90, 0.40, '38', fontsize=fontsize)
thinkplot.Config(xlim=[34, 97], ylim=[0, 1], legend=True, loc='upper left',
                 xlabel='cohort (decade)', ylabel='Fraction unmarried', 
                 title='Men in the U.S.')
thinkplot.Save(root='figs/marriage9', **options)
Explanation: Marriage patterns for men
End of explanation
<END_TASK>