<jupyter_start><jupyter_text>L&S 88 - Final Project [Sample] - Insurance Charges and Smoking Habits_<NAME>, [other teammates]_In this notebook, we will attempt to reproduce the result of the analysis by <NAME> on the relationship between insurance charges and smoking habits. The original analysis is hosted [on Github](https://github.com/chrispyles/jupyter/tree/master/content/notebooks/insurance).**Table of Contents**1. [Exploratory Data Analysis](Exploratory-Data-Analysis)2. [Reproducing the Original A/B Test](Reproducing-the-Original-A/B-Test)3. [A/B Testing Without Outliers](A/B-Testing-Without-Outliers)4. [Reliability of the Data](Reliability-of-the-Data)5. [Conclusions](Conclusions)6. [Sources](Sources)<jupyter_code>from datascience import *
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
plt.style.use('fivethirtyeight')<jupyter_output><empty_output><jupyter_text>Exploratory Data AnalysisWe begin by loading in the dataset which was downloaded from Github.<jupyter_code>data = Table.read_table('../data/insurance.csv')
data.show(5)<jupyter_output>/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/datascience/tables.py:132: FutureWarning: read_table is deprecated, use read_csv instead.
  df = pandas.read_table(filepath_or_buffer, *args, **vargs)<jupyter_text>In order to better understand the types of variables that are recording in the `data` table, we undertake some EDA below. We first take a look at the unique values in the `region` column:<jupyter_code>pd.Series(data.column('region')).unique()<jupyter_output><empty_output><jupyter_text>Based on this, `region` appears to be a categorical variable. According to the documentation for the original data set on Kaggle, this column describes the region of the US in which the individual lives. The `sex` and `smoker` columns are self-explanatory, and contains only the expected values:<jupyter_code>print('Unique values for sex:', np.unique(data['sex']))
print('Unique values for smoker:', np.unique(data['smoker']))<jupyter_output>Unique values for sex: ['female' 'male']
Unique values for smoker: ['no' 'yes']<jupyter_text>To better understand the numerical values in the table, we look at some important statistics for each:<jupyter_code>data.drop('sex', 'smoker', 'region').stats(ops = [min, max, np.mean, np.std])<jupyter_output><empty_output><jupyter_text>Because we are interested in the variables `charges` and `smoking`, we plot box plots below for the insurance charges separated into smoker and non-smoker groups.<jupyter_code>smoker = data.where('smoker', 'yes').column('charges')
nonsmoker = data.where('smoker', 'no').column('charges')
plt.boxplot([smoker, nonsmoker])
plt.title('Boxplots of Smoker and Nonsmoker Charges')
plt.ylabel('Charges ($)')
plt.xlabel('Smoker; 1 = Y, 2 = N')
plt.savefig("../output/smoker-boxplots-datascience.png"); # export to output folder<jupyter_output><empty_output><jupyter_text>Reproducing the Original A/B TestIn this section of the notebook, we will attempt to reproduce the original A/B test that was run on the data. The test statistic that was used was the absolute difference between the mean charges for the smoker and non-smoker groups. That is, for two vectors of $n$ smoker charges $\mathbf{s}$ and $m$ non-smoker charges $\mathbf{\bar{s}}$, the test statistic is computed as$$ \left | \frac{1}{n} \sum_{i=1}^n s_i - \frac{1}{m} \sum_{i=1}^m \bar{s}_i \right | $$To accomplish this, we first define a function `test_statistic` that takes in a table with columns `smoker` and `charges` and computes the value of the test statistic.<jupyter_code>def test_statistic(tbl):
    grouped = (tbl
               .group('smoker', np.mean)             # group by smoker values and take mean of cols
               .select('smoker', 'charges mean'))    # select smoker and charges mean cols
    smoker = (grouped
              .where('smoker', 'yes')                # select the row where smoker is 'yes'
              .column('charges mean'))[0]            # take that value of charges mean col
    nonsmoker = (grouped
                 .where('smoker', 'no')              # select the row where smoker is 'yes'
                 .column('charges mean'))[0]         # take that value of charges mean col
    
    return abs(smoker - nonsmoker)                   # return absolute difference

test_statistic(data)    # running the function on the original data<jupyter_output><empty_output><jupyter_text>We now turn our attention to the A/B test. The purpose of the A/B test is to determine whether or not the charges for smokers and non-smokers come from the same underlying distribution. For this A/B test the null and alternative hypotheses are as follows:* $H_0$: There is not difference between the distributions of charges for smokers and non-smokers. Any difference observed is due to random chance.* $H_A$: Being a smoker affects what you are charged in insurance premia; that is, the smoker and non-smoker charges are not from the same distribution.To check the intuition of $H_A$, we first plot histograms of the data on the same axes to see if the distributions match visually.<jupyter_code>plt.hist(smoker, density=True, alpha=0.5)       # reusing the smoker array created earlier
plt.hist(nonsmoker, density=True, alpha=0.5)    # reusing the nonsmoker array created earlier
plt.title('Histograms of Smoker and Non-Smoker Charges')
plt.xlabel('Charges ($)')
plt.legend(['Smoker', 'Non-Smoker'])
plt.savefig("../output/charges-hist-dstascience.png"); # export to output folder<jupyter_output><empty_output><jupyter_text>Based on the above visualization, it appears our intuition is correct, but we should nevertheless verify this statistically. The A/B test will occur in steps that will culminate in calculating a p-value to determine whether we accept $H_0$ or $H_A$. The steps we will follow are enumerated below:1. Define a function `shuffle` to shuffle the table's smoker classifications2. Repeatedly shuffle the table and apply the function `test_statistic`3. Calculate the proportion of values of the test statistic that are **greater than or equal to** the value for the original tableTo begin, we define the function `shuffle` below which shuffles the table.<jupyter_code>def shuffle(tbl):
    new_tbl = tbl.copy()    # copy the table so we don't mutate the original
    smoker_status = (new_tbl
                     .sample(with_replacement=False)    # shuffle the table by sampling without replacement
                     .column('smoker'))                 # select smoker column from shuffle table
    new_tbl['smoker'] = smoker_status
    return new_tbl

shuffle(data).show(5)    # test the function on the original data table<jupyter_output><empty_output><jupyter_text>For the second step of the process, we use a for loop below to run through the repetitions. A list is used instead of an array below to make the code run faster.<jupyter_code>np.random.seed(13)       # seeding np.random to make results reproducible
reps = 10000             # number of repetitions to run through
test_statistics = []     # initialize array to collect values
for _ in np.arange(reps):
    shuffled_data = shuffle(data)                 # shuffle data table
    test_stat = test_statistic(shuffled_data)     # calculate test statistic on shuffled table
    test_statistics += [test_stat]                # add value to test_statistics list<jupyter_output><empty_output><jupyter_text>Finally, we calculate the p-value below.<jupyter_code># calculate value of test statistic for original data
original = test_statistic(data)

# calculate p-value
p_value = sum(test_statistics >= original) / len(test_statistics)     

print("The p-value for this A/B test is {}.".format(p_value))<jupyter_output>The p-value for this A/B test is 0.0.<jupyter_text>The p-value is apparently 0, which is significant by any convention. So, in this case, we accept $H_A$, that the smoker and non-smoker charges do not come from the same underlying distribution.  A/B Testing Without OutliersThe original analysis fails to consider the effect of outliers on the results that it derives. In order to see if this has an effect, we now run the same test on a data set where we have taken out the outliers _within each smoking group_. That is, we first separate the individuals into smokers and non-smokers, and then remove outliers from each subset.For the purposes of this analysis, we define outliers to be values that are more than 2.5 standard deviations ($\sigma$) away from the mean ($\mu$). That is, the set of outliers $O$ is$$O = \{ c_i \mid c_i  \mu + 2.5 \sigma \}$$In the cell below, we remove these data from the data table.<jupyter_code>all_smokers = data.where('smoker', 'yes')
mu_smokers, sigma_smokers = np.mean(all_smokers['charges']), np.std(all_smokers['charges'])
smokers = all_smokers.where('charges', are.between_or_equal_to(
    mu_smokers - 2.5 * sigma_smokers,
    mu_smokers + 2.5 * sigma_smokers
))

print('Number of removed smokers: {}'.format(all_smokers.num_rows - smokers.num_rows))
all_non_smokers = data.where('smoker', 'no')
mu_non_smokers, sigma_non_smokers = np.mean(all_non_smokers['charges']), np.std(all_non_smokers['charges'])
non_smokers = all_non_smokers.where('charges', are.between_or_equal_to(
    mu_non_smokers - 2.5 * sigma_non_smokers,
    mu_non_smokers + 2.5 * sigma_non_smokers
))

print('Number of removed non-smokers: {}'.format(all_non_smokers.num_rows - non_smokers.num_rows))<jupyter_output>Number of removed non-smokers: 41<jupyter_text>Because we removed outliers in the data set, we now recreate the data table without outliers from the `smokers` and `non_smokers` dataframes and run the A/B test.<jupyter_code>no_outliers = smokers.with_rows(non_smokers.rows)
print(no_outliers.num_rows == data.num_rows - 43)
np.random.seed(13)       # seeding np.random to make results reproducible
reps = 10000             # number of repetitions to run through
no_outliers_test_statistics = []     # initialize array to collect values
for _ in np.arange(reps):
    shuffled_data = shuffle(no_outliers)          # shuffle data table
    test_stat = test_statistic(shuffled_data)     # calculate test statistic on shuffled table
    no_outliers_test_statistics += [test_stat]    # add value to test_statistics list<jupyter_output><empty_output><jupyter_text>Finally, we calculate the p-value below.<jupyter_code># calculate value of test statistic for original data
no_outliers_original = test_statistic(no_outliers)

# calculate p-value
no_outliers_p_value = sum(no_outliers_test_statistics >= no_outliers_original) \
        / len(no_outliers_test_statistics)     

print("The p-value for this A/B test is {}.".format(no_outliers_p_value))<jupyter_output>The p-value for this A/B test is 0.0.
{'fraction_non_alphanumeric': 0.09297807382736609, 'fraction_numerical': 0.006383569247849014, 'mean_word_length': 4.529411764705882, 'pattern_counts': {'":': 0, '<': 54, '<?xml version=': 0, '>': 56, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '14', 'dir': 'jupyter-structured-clean-dedup', 'id': '585099', 'license': 'CC0-1.0', 'n_tokens_mistral': 3381, 'n_tokens_neox': 3170, 'n_words': 1153, 'path': 'notebooks/insurance-datascience.ipynb', 'repo_name': 'zehao-zhao/datah195a', 'size': '10819'}