<jupyter_start><jupyter_text>NP plots<jupyter_code>revdf = revdf.copy()

revdf['quantile_ex_2'] = pd.qcut(revdf['REVENUE_AMT'], q=10, precision=0)
# describe(include='category')
test = revdf['quantile_ex_2'].value_counts()
# test = revdf.quantile_ex_2.value_counts()
test.index[4]
test2 = pd.DataFrame(test)

plt.bar(test2.index, test2.values)
#establish df for plotting num of NPs per State > city > zip
#plot on folium density map
#number of NPs registered with IRS by state


# #dual axis bar plot (reformat values to show in millions)

count = revdf.groupby('STATE_x').count()[['REVENUE_AMT']]
revs = revdf.groupby('STATE_x').sum()[['REVENUE_AMT']]

combined = pd.merge(count,revs, left_index=True,right_index=True).sort_values(by='REVENUE_AMT_x', ascending=False)

combined.REVENUE_AMT_y = round(combined['REVENUE_AMT_y']/1000000,2)

def dual_axis_plot(df, y1, y2, title, y1label, y2label, xlabel,filename):
    """function to plot bar and line on dual y-axis."""
    fig, ax1 = plt.subplots(figsize=(12,6))

    ax2 = ax1.twinx()
    
    ax1.bar(df.index, y1,label='Count')
    ax1.set_xticklabels(df.index, rotation=90)
    ax1.set_xlabel(f'{xlabel}', size=12)
    ax1.set_ylabel(f'{y1label}', size=12)
    ax1.set_title(f'{title}', size=15)
    ax1.grid()
    
    ax2.bar(df.index, y2, color='gold', alpha=.3, width=.5, align='center', label='Revenues')
    ax2.set_ylabel(f'{y2label}')
    ax2.grid()
    fig.legend(bbox_to_anchor=(.5, .85))
    plt.savefig(f'../images/{filename}.png')

    
dual_axis_plot(combined, 
               combined['REVENUE_AMT_x'], 
               combined['REVENUE_AMT_y'],
              'STATES - NONPROFITS REGISTERED (2019)',
              'COUNT',
              'REVENUES ($ IN MILLIONS)',
              'STATES',
              'num_nps2')
##### establish df for plotting num of NPs per State > city > zip
#plot on folium density map
#number of NPs registered with IRS by city

count2 = revdf.groupby('CITY').count()[['REVENUE_AMT']]
revs2 = revdf.groupby('CITY').sum()[['REVENUE_AMT']]

combined2 = pd.merge(count2,revs2, left_index=True,right_index=True).nlargest(20, columns='REVENUE_AMT_x').sort_values(by='REVENUE_AMT_x', ascending=False)

combined2.REVENUE_AMT_y = round(combined2['REVENUE_AMT_y']/1000000,2)


dual_axis_plot(combined2, 
               combined2['REVENUE_AMT_x'], 
               combined2['REVENUE_AMT_y'],
              'TOP 20 CITIES - NONPROFITS REGISTERED (2019)',
              'COUNT',
              'REVENUES ($ IN MILLIONS)',
              'CITY',
              'num_nps_city2')

new_index = [ 3,  9,  4, 12,  6,  5, 14,  8,  7, 27, 13,  2, 19, 25, 29, 10,  1,
            17, 11, 18, 23, 26, 15, 91, 16, 82,  0, 50, 20]


index_list = []

for i in new_index:
    index_list.append(f'501(c){i}')

#nps by subsection, relevant because this will tell me 
#type of organization and if contributions are deductible
#deductible is not beind analyzed at this moment

#which subsection have the most NPs and where are they in the US


count3 = revdf.groupby('SUBSECTION').count()[['REVENUE_AMT']]
revs3 = revdf.groupby('SUBSECTION').sum()[['REVENUE_AMT']]

combined3 = pd.merge(count3,revs3, left_index=True,right_index=True).nlargest(10,columns='REVENUE_AMT_x').sort_values(by='REVENUE_AMT_y', ascending=False)

combined3.REVENUE_AMT_y = round(combined3['REVENUE_AMT_y']/10000000,2)

# combined3.index.astype('str', copy=False)
combined3.index = combined3.index.astype(str, copy=False)
dual_axis_plot(combined3, 
               combined3['REVENUE_AMT_x'], 
               combined3['REVENUE_AMT_y'],
              'TOP 10 SUBSECTIONS 501(c)X - NUMBER OF NONPROFITS REGISTERED(2019)',
              'COUNT',
              'REVENUES ($ IN 10s MILLIONS)',
              'SUBSECTION (501(c)X)',
              'num_nps_sub2')
#nps by foundation, where they receive most of their funding

#which foundations have the most NPs and where are they in the US

count4 = revdf.groupby('FOUNDATION').count()[['REVENUE_AMT']]
revs4 = revdf.groupby('FOUNDATION').sum()[['REVENUE_AMT']]

combined4 = pd.merge(count4,revs4, left_index=True,right_index=True).nlargest(10,columns='REVENUE_AMT_y').sort_values(by='REVENUE_AMT_y', ascending=False)

combined4.REVENUE_AMT_y = round(combined4['REVENUE_AMT_y']/1000000,2)


combined4.index = combined4.index.astype(str, copy=False)

dual_axis_plot(combined4, 
               combined4['REVENUE_AMT_x'], 
               combined4['REVENUE_AMT_y'],
              'TOP 5 FOUNDATION CATEGORIES - NUMBER OF NONPROFITS REGISTERED (2019)',
              'COUNT',
              'REVENUES ($ IN MILLIONS)',
              'FOUNDATION CATEGORY',
              'num_nps_found2')
sub501c3 = revdf[revdf['SUBSECTION'] == 3]


count5 = sub501c3.groupby('STATE_x').count()[['REVENUE_AMT']]
revs5 = sub501c3.groupby('STATE_x').sum()[['REVENUE_AMT']]

combined5 = pd.merge(count5,revs5, left_index=True,right_index=True).sort_values(by='REVENUE_AMT_y', ascending=False)
#.nlargest(10,columns='REVENUE_AMT_y')
combined5.REVENUE_AMT_y = round(combined5['REVENUE_AMT_y']/1000000,2)


combined5.index = combined5.index.astype(str, copy=False)

dual_axis_plot(combined5, 
               combined5['REVENUE_AMT_x'], 
               combined5['REVENUE_AMT_y'],
              '501(c)3 NONPROFITS REGISTERED (2019) by STATE',
              'COUNT',
              'REVENUES ($ IN MILLIONS)',
              'STATE',
              'num_nps_sub2')
sub501c3 = subdf[subdf['FOUNDATION'] == 15]


count5 = sub501c3.groupby('STATE_x').count()[['REVENUE_AMT']]
revs5 = sub501c3.groupby('STATE_x').sum()[['REVENUE_AMT']]

combined5 = pd.merge(count5,revs5, left_index=True,right_index=True).sort_values(by='REVENUE_AMT_y', ascending=False)
#.nlargest(10,columns='REVENUE_AMT_y')
combined5.REVENUE_AMT_y = round(combined5['REVENUE_AMT_y']/1000000,2)


combined5.index = combined5.index.astype(str, copy=False)

dual_axis_plot(combined5, 
               combined5['REVENUE_AMT_x'], 
               combined5['REVENUE_AMT_y'],
              '501(c)3-15 NONPROFITS REGISTERED (2019) by STATE',
              'COUNT',
              'REVENUES ($ IN MILLIONS)',
              'STATE',
              'num_nps_subfound2')<jupyter_output><empty_output><jupyter_text>Income Plots<jupyter_code>x = foudf['avg_income']


plt.figure(figsize=(12,4))
plt.xlim(10, 250)
plt.axvline(x.median(), color='gold', label='median income')
plt.hist(x, bins=1000)
plt.xlabel('Income', size=12)
plt.ylabel('Count', size=12)
plt.title('Individual Income Reported to IRS 2019', size=15)
plt.grid()
plt.legend()
plt.savefig('../images/hist_ind_income3.png')

#pareto distribution which is very typical for income 80% pop only has 20% of wealth/income
fig, axs = plt.subplots(4,3, figsize=(14,14))

test_state_lst = ['CA','NY','FL', 'PA', 'DC','TX','OH', 'IL','VA','AZ', 'MD','WA']

for i, axs in zip(test_state_lst, axs.flatten()):

    state = revdf[revdf['STATE_x']==f'{i}']
    
    x = state['avg_income']
    
    axs.hist(x,bins=20)
    axs.axvline(x.median(), color='gold', label='median income')
    axs.set_title(f'{i}', size=12)
    axs.grid()
    axs.legend()
    fig.tight_layout
    fig.suptitle('Individual Income Distribution by Top 12 States', y=0.93,
                 size='xx-large', weight='bold',
                 stretch='expanded')
    fig.text(0.5, 0.1, 'Income (thousands)', ha='center', va='center',
             size='xx-large', stretch='semi-expanded')
    fig.text(0.06, 0.5, 'Count', ha='center', va='center',
             rotation='vertical', size='xx-large',
             stretch='semi-expanded')
    plt.savefig('../images/state_inc_hist3.png')<jupyter_output><empty_output><jupyter_text>NP vs. Income PlotsThere we begin to only focus on the nonprofits who have reported income to the IRS, i.e. including only revenues greater than 0. Also, we onyl analyze the Subsection 501(c)3 which are the majority of the nonprofits and revenues. Outliers, such as subsection 12, hospitals and medical research organizations were excluded due to having the amount of revenues they generate vs number of organizations. They would heavily skew our metrics.<jupyter_code>test_state_lst = ['CA','NY','FL', 'PA', 'DC','TX','OH', 'IL','VA','AZ', 'MD','WA']

fig, axs = plt.subplots(2,6, figsize=(22,10))



for i, axs in zip(test_state_lst, axs.flatten()):

    state = foudf[foudf['STATE_x']==f'{i}']
    
    x = state['avg_income']
    y = state['REVENUE_AMT']/1000000
    
    axs.scatter(x,y, alpha=.3)
    axs.set_title(f'{i}', size=12)
    axs.grid()
    fig.tight_layout
    fig.suptitle('NonProfit Revenue on Individual Income by Top 12 States', y=0.93,
                 size='xx-large', weight='bold',
                 stretch='expanded')
    fig.text(0.5, 0.08, 'Income (thousands)', ha='center', va='center',
             size='xx-large', stretch='semi-expanded')
    fig.text(0.08, 0.5, 'Revenues (millions)', ha='center', va='center',
             rotation='vertical', size='xx-large',
             stretch='expanded')
    plt.savefig('../images/rev_on_inc_scatter4.png')<jupyter_output><empty_output><jupyter_text>Scatter of revenue on individual incomeThese plots shows a very dense concentration on lower level revenue NPs. Since we are analyzing 501c3-15 only, these organizations only receive a substantial amount of support from governmental units or general public. This group is very interesting to me since, technically government funded organizations are general public funds since it's tax collected from the American citizen. The sole difference is the discretion lies on the government to donate those funds vs. the individual donating thier discretionary income to NPs. This is worth a deeper look. Also, it might be worth looking into if there is some relationship between income and lower revenue generating NPs. Large NPs understanably receive substantial government funds which can skew the story.<jupyter_code>revs = foudf['REVENUE_AMT'].sort_values(ascending=False)


x = np.logspace(1,10)

plt.figure(figsize=(12,4))
plt.hist(revs, bins=x, density=True)
plt.axvline(revs.mean(), color='gold', label='median revenues')
plt.xscale('log')
plt.title('501(c)3-15 NONPROFIT REVENUE DISTRIBUTION in 2019')
plt.ylabel('COUNT')
plt.xlabel('REVENUES (log scale)')
plt.legend()
plt.grid()
plt.show()
# plt.savefig('../images/np_rev6.png')
fig, axs = plt.subplots(12,1, figsize=(14,30))


for i, axs in zip(test_state_lst, axs.flatten()):
    
    scale = np.logspace(1,10)
    
    state = foudf[foudf['STATE_x']==f'{i}']
    
    x = state[state['REVENUE_AMT'] > 0]['REVENUE_AMT'].sort_values(ascending=False)
    
    axs.hist(x,bins=scale)
    axs.axvline(x.median(), color='gold', label='median revenue')
    axs.set_title(f'{i}', size=12)
    axs.set_xscale('log')
    axs.grid()
    axs.legend()
    fig.tight_layout
    fig.suptitle('501(c)3-15 Revenue Distribution by Top 12 States', y=0.93,
                 size='xx-large', weight='bold',
                 stretch='expanded')
    fig.text(0.5, 0.1, 'Revenue (log scale)', ha='center', va='center',
             size='xx-large', stretch='semi-expanded')
    fig.text(0.06, 0.5, 'Count', ha='center', va='center',
             rotation='vertical', size='xx-large',
             stretch='semi-expanded')
    plt.savefig('../images/state_rev_hist2.png')

fig, axs = plt.subplots(4,3, figsize=(14,14))


for i, axs in zip(test_state_lst, axs.flatten()):
    
    scale = np.logspace(1,10)
    
    state = foudf[foudf['STATE_x']==f'{i}']
    
    x = state[state['REVENUE_AMT'] > 0]['REVENUE_AMT'].sort_values(ascending=False)
    
    axs.hist(x,bins=scale)
    axs.axvline(x.median(), color='gold', label='median revenue')
    axs.set_title(f'{i}', size=12)
    axs.set_xscale('log')
    axs.grid()
    axs.legend()
    fig.tight_layout
    fig.suptitle('501(c)3-15 Revenue Distribution by Top 12 States', y=0.93,
                 size='xx-large', weight='bold',
                 stretch='expanded')
    fig.text(0.5, 0.1, 'Revenue (log scale)', ha='center', va='center',
             size='xx-large', stretch='semi-expanded')
    fig.text(0.06, 0.5, 'Count', ha='center', va='center',
             rotation='vertical', size='xx-large',
             stretch='semi-expanded')
    plt.savefig('../images/state_rev_hist4.png')<jupyter_output><empty_output><jupyter_text>stats analysis<jupyter_code># eda_df.info()
a = revdf[revdf['avg_income'] < revdf['avg_income'].median()]['REVENUE_AMT']
b = revdf[revdf['avg_income'] >= revdf['avg_income'].median()]['REVENUE_AMT']

# ttest, pvalue = stats.ttest_ind(a,b)

ttest, pvalue = stats.ttest_ind(a,b)
pvalue
t = foudf[foudf.REVENUE_AMT > 100000000]['REVENUE_AMT'].sum()

t2 = foudf[foudf.REVENUE_AMT <= 100000000]['REVENUE_AMT'].sum()

t, t2
x = np.logspace(1,10)

plt.figure(figsize=(12,4))
plt.xscale('log')
plt.hist(a, bins=x, alpha=.3)
plt.hist(b, bins=x, alpha=.3, color='gold')
plt.grid()
plt.show()
plt.savefig('../images/np_rev_dist.png');

#the probability of observing revenues above the median based on income is nearly 0. 

#the probability of observing revenues above the mean based on income is nearly 0 as well. 
#as the threshold 
# fig, axs = plt.subplots(4,3, figsize=(14,14))

# # for i, axs in zip(test_state_lst, axs.flatten()):
    
# #     state = foudf[foudf['STATE_x']==f'{i}']

# #     a = foudf[foudf['avg_income'] < foudf['avg_income'].mean()]['REVENUE_AMT']
# #     b = foudf[foudf['avg_income'] >= foudf['avg_income'].mean()]['REVENUE_AMT']


# #     min_x = a2.mean() - (a2.std()/(len(a2)**.5)*4)
# #     max_x = a2.mean() + (a2.std()/(len(a2)**.5)*10)

# #     min_y = b2.mean() - (b2.std()/(len(b2)**5)*4)
# #     max_y = b2.mean() + (b2.std()/(len(b2)**5)*10)


# #     z = np.linspace(min(min_x, min_y), max(max_x, max_y))

# #     one = stats.norm(loc=np.mean(a2), scale=np.std(a2)/(len(a2)**.5))
# #     two = stats.norm(loc=np.mean(b2), scale=np.std(b2)/(len(b2)**.5))


# #     axs.plot(z, one.pdf(z))
# #     axs.plot(z, two.pdf(z), color='gold')
# #     axs.grid()
    
# #     fig.tight_layout
# Plot the distributions of the samples and compare
a = foudf[foudf['avg_income'] < foudf['avg_income'].mean()]['REVENUE_AMT']
b = foudf[foudf['avg_income'] >= foudf['avg_income'].mean()]['REVENUE_AMT']


min_x = a2.mean() - (a2.std()/(len(a2)**.5)*4)
max_x = a2.mean() + (a2.std()/(len(a2)**.5)*10)

min_y = b2.mean() - (b2.std()/(len(b2)**5)*4)
max_y = b2.mean() + (b2.std()/(len(b2)**5)*10)


z = np.linspace(min(min_x, min_y), max(max_x, max_y))

one = stats.norm(loc=np.mean(a2), scale=np.std(a2)/(len(a2)**.5))
two = stats.norm(loc=np.mean(b2), scale=np.std(b2)/(len(b2)**.5))

fig, ax = plt.subplots(figsize=(6,6))

ax.plot(z, one.pdf(z))
ax.plot(z, two.pdf(z), color='gold')
ax.grid()

def dual_axis_plot(df, y1, y2, title, y1label, y2label, xlabel,filename):
    """function to plot bar and line on dual y-axis."""
    fig, ax1 = plt.subplots(figsize=(12,6))

    ax2 = ax1.twinx()
    
    ax1.bar(df.index, y1,label='Count')
    ax1.set_xticklabels(df.index, rotation=90)
    ax1.set_xlabel(f'{xlabel}', size=12)
    ax1.set_ylabel(f'{y1label}', size=12)
    ax1.set_title(f'{title}', size=15)
    ax1.grid()
    
    ax2.bar(df.index, y2, color='gold', alpha=.3, width=.5, align='center', label='Revenues')
    ax2.set_ylabel(f'{y2label}')
    ax2.grid()
    fig.legend(bbox_to_anchor=(.5, .85))
    plt.savefig(f'../images/{filename}.png')


fig, ax1 = plt.subplots(12,1, figsize=(14,45))


for i, ax1 in zip(test_state_lst, ax1.flatten()):
    
    state = foudf[foudf['STATE_x']==f'{i}']

    a = state[state['avg_income'] < state['avg_income'].median()]['REVENUE_AMT']
    b = state[state['avg_income'] >= state['avg_income'].median()]['REVENUE_AMT']
    
    ttest, pvalue = stats.ttest_ind(a,b)
    x = np.logspace(-1,10)
    
    ax2 = ax1.twinx()
    

    ax1.hist(a, bins=x, alpha=.3)
    ax1.set_title(f'{i} Non-Profit Revenue Distribution: T-Test={round(ttest,2)}, p-value={round(pvalue,4)}')
    ax1.set_xscale('log')
    ax1.grid()
    
    ax2.hist(b, bins=x, alpha=.3, color='gold')
    ax2.grid()
    ax2.set_xscale('log')
    
    fig.tight_layout

state = df[df['STATE_x']=='DC']

a = state[state['avg_income'] < state['avg_income'].median()]['REVENUE_AMT']
b = state[state['avg_income'] >= state['avg_income'].median()]['REVENUE_AMT']

a2 = a[a>0]
b2 = b[b>0]

ttest, pvalue = stats.ttest_ind(a2,b2)

ttest, pvalue


for i in test_state_lst:
    
    state = df[df['STATE_x']==f'{i}']

    a = state[state['avg_income'] < state['avg_income'].median()]['REVENUE_AMT']
    b = state[state['avg_income'] >= state['avg_income'].median()]['REVENUE_AMT']
    
    a2 = a[a>0]
    b2 = b[b>0]
    
    med = state['REVENUE_AMT'].mean()
    stda2 = a2.std()
    meda2 = a2.mean()
    ttest, pvalue = stats.ttest_ind(a2,b2)
subsection = df[df.SUBSECTION == 3]
fig, axs = plt.subplots(11,1, figsize=(14,40))

test_sub_lst = [3,4,9,10,11,12,13,14,15,17,18]

for i, axs in zip(test_sub_lst, axs.flatten()):
    
    scale = np.logspace(1,10)
    
    foundation = subsection[subsection['FOUNDATION']==i]
    
    x = foundation[foundation['REVENUE_AMT'] > 0]['REVENUE_AMT'].sort_values(ascending=False)
    
    axs.hist(x,bins=scale)
    axs.axvline(x.mean(), color='gold', label='mean revenue')
    axs.set_title(f'Category 501(c){i}', size=12)
    axs.set_xscale('log')
    axs.grid()
    axs.legend()
    fig.tight_layout
    fig.suptitle('Revenue Distribution by IRS Subsection Category', y=0.93,
                 size='xx-large', weight='bold',
                 stretch='expanded')
    fig.text(0.5, 0.1, 'Revenue (log scale)', ha='center', va='center',
             size='xx-large', stretch='semi-expanded')
    fig.text(0.06, 0.5, 'Count', ha='center', va='center',
             rotation='vertical', size='xx-large',
             stretch='semi-expanded')
    plt.savefig('../images/sub_rev2.png')

# Run a Welch's T-test on all samples to evaluate differences.
def ttest_for_rap(rev_list, alpha=0.05):
    '''
    Given a list of name-array tuples (created in df_to_scores), this 
    function runs a Welch's t-test between rap and all other genres with
    a Bonferroni correction.

    Parameters:
    score_list (list of tuples): list of name-series tuples
    alpha (float): desired alpha for each individual t-test

    Returns:
    diff_list (list): list of genres for which we reject the null
    same_list (list): list of genres for which we fail to reject the null
    '''
    diff_list = []
    same_list = []
    bonf = 24
    for tup in score_list:
        if tup[0] == 'rap':
            rap_scores = tup[1]
    for tup in score_list:
        if tup[0] == 'rap':
            continue
        else:
            tstat, pvalue = stats.ttest_ind(rap_scores, tup[1],
                                            equal_var=False)
            if pvalue < alpha/2/bonf:
                diff_list.append(tup[0])
            else:
                same_list.append(tup[0])
    return diff_list, same_list<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1689776981072389, 'fraction_numerical': 0.03759160647440291, 'mean_word_length': 3.6535819430814525, 'pattern_counts': {'":': 0, '<': 29, '<?xml version=': 0, '>': 39, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 3, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '27722', 'license': 'Naumen', 'n_tokens_mistral': 7803, 'n_tokens_neox': 7214, 'n_words': 1357, 'path': 'notebooks/.ipynb_checkpoints/day3_eda_-checkpoint.ipynb', 'repo_name': 'lesro/change_unlimited', 'size': '18967'}