<jupyter_start><jupyter_text>**Data Science in Python workshop** <NAME>, Applied Science @ Discovery Channel___ ** teaching materials in these notebooks are solely the property of <NAME>, and cannot be used for other workshops or courses other than by the author or by license. Students are open to review these notebooks at any time, but with attribution of credit back to the author.  Learning Objectives- Review the relation between correlation and causation- Understand what a hypothesis test is, and why it is important- Identify what missing data is and how to handle it- Test a hypothesis using a sample case study--- Data SourceToday, we’ll use advertising data from an example in the book [An Introduction to Statistical Learning](http://www-bcf.usc.edu/~gareth/ISL/).- This is a well-known, standard introduction to machine learning.- The book has a more advanced version — [Elements of Statistical Learning](http://web.stanford.edu/~hastie/ElemStatLearn/) — if you are comfortable with linear algebra and statistics at the graduate level.<jupyter_code>#import pandas, numpy, and matplotlib packages, 
#include jupyter magic
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

#this allows us to view visualizations in our notebook
%matplotlib inline
#For google colab, we need to upload files - there are other methods worth exploring on your own. 
from google.colab import files
uploaded = files.upload()
#Here I'm importing that data table
df = pd.read_csv('TV-Radio-Newspaper-Sales.txt',  index_col = 0)
df.info()
#Show the head
df.head()
#describe your data using the pandas method for summary statistics
df.describe()<jupyter_output><empty_output><jupyter_text>The Normal distribution---The Normal distribution is arguably the most commonly used distribution in all of statistics. **Normality** is an assumption that underlies many statistical tests and serves as a convenient model for the distribution of many (but not all!) variables.The Normal distribution is parameterized by two parameters: - the population mean- the population standard deviation. If a variable follows a Normal distribution exactly, its mean, median, and mode will all be equal.<jupyter_code>#We can approximate whether we have a normal distribution using .describe() by comparing the mean and the median
df.sales.describe()
#We can also visualize it!
plt.hist(df['sales'], bins = 5);<jupyter_output><empty_output><jupyter_text>Here we are starting to do some statistical computing. We might think this is very abnormally distributed, but it's because we don't have the right bins factor. [Read this post if you want to learn all of the different rules for binning](https://medium.com/analytics-vidhya/data-scientists-stop-randomly-binning-histograms-1069d7380c3a)Typically, though, you should start with scott's rule. Let's look at the pieces of statistical knowledge we need to get there!  Outliers and Standard DeviationAn outlier in a data series is a signficant deviation from some reasonable value. Let's start with a 'small data' example. In this dataframe, the 'count' is the number of 'values' in the dataset -- this is essentially a visualization of a `value_count()` applied to sales.<jupyter_code>count = [2, 3, 5, 2, 3, 12, 5, 3, 4]
value = [1, 2, 3, 4, 5, 6, 7, 8, 9]

dataframe = pd.DataFrame(list(zip(value, count)), columns = ['value', 'count'])
plt.scatter(x = dataframe['value'], y = dataframe['count'])
plt.title('Number of sales (y) worth specific values (x)');<jupyter_output><empty_output><jupyter_text>It's pretty intuitive for us humans to call out which of these does not belong. But computers don't have intuition, so we need to help it out. To find the anomaly in a series of numbers, we need to first define what a reasonable value is. A good place to start looking is the mean!<jupyter_code>dataframe['count'].mean()<jupyter_output><empty_output><jupyter_text>... and then, we need to define the deviation. The standard deviation is an easy measure for this. Let's implement this by hand! First, we calculate variance - the average squared distance from the mean. Then, to get to the standard deviation, we simply take the square root of the variance.<jupyter_code>dataframe.head()
#start with an empty list to store our distance calcuations
distance_calc = []

#We start with the top right of the equation: 
## squared value of xi (x at index i) minus the mean of x (x-bar)
for cnt in dataframe['count']:
  distance_calc.append((cnt - dataframe['count'].mean()) ** 2)

#take the sum of that list (Sigma = sum of...), and divide by the length of the results. The is stored as s^2
s2 = round(sum(distance_calc) / (len(distance_calc)-1), 2)

#Take the square root of s2 to get to the standard deviation
standard_deviation = np.sqrt(s2)

standard_deviation<jupyter_output><empty_output><jupyter_text>Now that we've defined a reasonable value and a deviation, we can define a range of acceptable, expected values. as something that is within 1 standard deviation: mean +- 1*std. Any value outside of that is considered an outlier.<jupyter_code>lower_bound = dataframe['count'].mean() - standard_deviation
upper_bound = dataframe['count'].mean() + standard_deviation

print('lowest expected value in a single deviation is {}'.format(lower_bound))
print('highest expected value in a single deviation is {}'.format(upper_bound))

#Create a feature that flags data that is outside of 1 standard deviation in dataframe
dataframe['outlier'] = ((dataframe['count'] < lower_bound) | (dataframe['count'] > upper_bound))
dataframe<jupyter_output><empty_output><jupyter_text>Using this approach, you should have found that value 6 (count 12) is outside the range of acceptable values. This is fairly imprecise though. What if we want to experiment with what the acceptable number of deviations is? We can instead return a Z-score. z-score, or **Standard Score**, is the number of standard deviations from the mean a value is. That means, it's obviously derivative from what we've just done! The Z score for a value is the value-mean, divided by the standard deviation - this returns how many standard deviations outside of the norm the value is.<jupyter_code>#lets implement this as a function that you'll apply to each cell of the dataframe!
def zScore(cell):
  '''
  Calculates the zscore, returning the zscore per instance. Accepts a single cell at a time. 
  '''
  return (cell - dataframe['count'].mean()) / standard_deviation
#this is the first time you're learning about applying functions - check out this syntax, and let's talk about it's power. 
dataframe['zScore'] = dataframe['count'].apply(zScore)
dataframe<jupyter_output><empty_output><jupyter_text>Let's say we wanted to flag any instances outside of 1.5 standard deviations as an outlier, we could create a boolean feature.<jupyter_code>#create a feature that says whether something is an outlier
dataframe['outlier'] = ((dataframe['zScore'] > 1.5 )| (dataframe['zScore'] < -1.5))<jupyter_output><empty_output><jupyter_text>If the business asked you "what percentage of our ideal dataframe is an outlier? Could you answer it now?<jupyter_code>#Check out all of your options in the .value_counts() method
dataframe['outlier'].value_counts(normalize = True)<jupyter_output><empty_output><jupyter_text>This is a great exercise, and having this skill will make you an easy top-hire for any data scientist. New algorithms are constantly being created, and if you can implement and make significant money for the business, you'll be well ahead of teams that have to wait for someone else to implement an algorithm in a standard package that passes security review. Let's practice using the sales data!<jupyter_code>#get a list of distances
distance_cal = []

#calculate the variance
for val in df.sales:
  distance_cal.append((val - df.sales.mean())**2)

#calculate the standard deviation
s2 = round(sum(distance_cal)/(len(df)-1), 2)

standard_deviation = np.sqrt(s2)

#define the upper and lower bounds for 1 standard deviation
upper_bound = df.sales.mean() + standard_deviation
lower_bound = df.sales.mean() - standard_deviation

#flag any values outside 1 standard deviation, name the feature stdFlag
df['stdFlag'] = (df.sales > upper_bound) | (df.sales < lower_bound)

df.head()
standard_deviation
df.stdFlag.value_counts()
#create a standard score function
def zScore(cell):
    return (cell - df['sales'].mean()) / standard_deviation
#apply standard score function to sales
df['zscore'] = df.sales.apply(zScore)
#flag any sales that are outside of TWO standard deviations
df['2_stdFlag'] = (df.zscore > 2) | (df.zscore < -2)
df['2_stdFlag'].value_counts()<jupyter_output><empty_output><jupyter_text>While useful as an exercise and skill, sometimes, we would just rather get through this work and onto actually developing insights out of the data. That's why we can use built in functions from pandas to help us with this work!<jupyter_code>#here, let's get the standard deviation. Does it match your answer above? Beware of rounding :)
df.sales.std()<jupyter_output><empty_output><jupyter_text>Now back to our histogram. The best way to bin a histogram so that you can evaluate whether the distribution is normal is by using scott rule. Scott rule is defined as a way to define the width of the bins in terms of values in the series.   Where  - R is`max(value) - min(value)` - n is the number of samples - σ is the standard deviationTo get to bin count, you then muliply by 1/Range.<jupyter_code>plt.hist(df.sales, bins = 'scott')<jupyter_output><empty_output><jupyter_text>Now you understand the math! Why do you think this presents the most empirical view of whether the distribution is normal? How does this visualization capture what you just learned? Below is a table comparing some example sample statistics and population parameters:Metric  | Statistic  | Parameter ------------- | --------------- | -------------mean   | $$\bar{x} = \frac{\sum x}{n}$$ | $$ \mu = \frac{\sum x}{N} $$ variance | $$ s^2 = \frac{\sum_i (x_i - \bar{x})^2}{n-1} $$ | $$ \sigma^2 = \frac{\sum_i (x_i - \mu)^2}{N}  $$standard deviation   | $$ s = \sqrt{\frac{\sum_i (x_i - \bar{x})^2}{n-1}} $$ | $$ \sigma = \sqrt{\frac{\sum_i (x_i - \mu)^2}{N} } $$correlation | $$ r = \frac{\sum(x_i - \bar{x})(y_i - \bar{y})}{(n-1)s_xs_y }$$ | $$r = \frac{\sum(x_i - \mu_x)(y_i - \mu_y)}{(n)\sigma_x\sigma_y }$$  Introduction to Hypothesis Testing---**Objective**: Test a hypothesis within a sample case study.You'll remember that we've worked previously on descriptive statistics such as mean and variance. How would we tell if there is a difference between our groups? How would we know if this difference was real or if our finding is simply the result of chance?Hypothesis Testing is a family of statistical methods used to identify whether a sample of observed data can be used to accept or reject a predefined hypothesis.  --- Imagine you’re checking whether an apple is good or bad, using a random selection of apples from a very large box of apples.In the below image, we see a clear effect of sample size on the measurements: the effect of one different apple is very large for small samples and the effect of 1 apple becomes smaller and smaller the larger our sample becomes.Understanding the effect of sample size is the first basis towards an understanding of Hypothesis Testing. We could start arguing that 0.5 on 2 apples, would be a 1 apple difference: very likely to happen. But on 100 apples, 0.5 would represent a 50 apple difference: an extremely strong difference!? Generally speaking, we start with a **null hypothesis** and an **alternative hypothesis**, which is the opposite of the null. Then, you check whether the data support rejecting your null hypothesis or fail to reject the null hypothesis.> - The __null hypothesis__ is, in general, the boring stuff i.e. it assumes that nothing interesting happens/happened.- The __alternative hypothesis__ is, where the action is i.e. some observation/ phenomenon is real (i.e. not a fluke) and statistical analysis will give us more insights on that.For example:> - **Null hypothesis**: There is no relationship between gender and sales.- **Alternative hypothesis**: There is a relationship between gender and sales.Note that "failing to reject" the null hypothesis is not the same as "accepting" it. Your alternative hypothesis may indeed be true, but you don't necessarily have enough data to show that yet.This distinction is important for helping you avoid overstating your findings. You should only state what your data and analysis can truly represent.  Example* The **null hypothesis**, denoted $H_0$, is that there's no difference in two **means** ($\mu_1 = \mu_2$) of a sample* The **alternative hypothesis**, denoted $H_1$, is that that there's a difference between the means of a sample. When doing calculations, we want to be very clear what we're testing. So, we need to define what we mean by 'difference' in our hypothesis test. Some differences can be:  * low-tailed: $H_1: \mu_1 > \mu_2$  * high-tailed: $H_1: \mu_1 < \mu_2$  * two-tailed: $H_1: \mu_1 \neq \mu_2$---Let's illustrate this with an example. * Hypothesis: There is a statistically significant difference the mean height of men and women. * Data: our data comprises a sample of men's heights and a sample of women's heights.* $H_0: \mu_{men} = \mu_{women}$* Alternative:  * low-tailed: $H_1: \mu_{men} > \mu_{women}$  * high-tailed: $H_1: \mu_{men} < \mu_{women}$  * two-tailed: $H_1: \mu_{men} \neq \mu_{women}$Wait... what do you mean statistically significant?  The p-ValueWe are trying to answer a crucial question: “What is the probability that observing a test-statistic, this extreme, for this sample (considering its size and a suitable probabilistic dynamics governing it), is the product of pure chance?"The simple description of the **p-value** is the probability that, **GIVEN THE NULL HYPOTHESIS IS TRUE**, we would observe the current set of data.  The actual definition: > For a given statistical model, **GIVEN THE NULL HYPOTHESIS IS TRUE**, the _p-value_ is the probability of observing a statistical summary at least as extreme as the actual observed results.If your p-value is less than a chosen significance level $\alpha$, then you may choose to reject the null hypothesis and accept there is reasonable evidence to support the alternative hypothesis.How we check that p-value depends on our research question and the relationship between our values, population, and sample! Let's walk through a few examples.<jupyter_code>#import helper functions
import statsmodels.api as sm

#set a random seed so we are all getting the same answers: 
np.random.seed = 42<jupyter_output>/usr/local/lib/python3.6/dist-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<jupyter_text>Hypothesis Tests Guide-Alternative hypothesis | test-- Mean of a variable is different from a reference value | One Sample t-test- Means of two groups are different | Two Sample t-test- Means of >2 gropus are different | ANOVA- % of success in a binary variable is different from a refernece value | Proportions z-test- % of success in a binary variable is different between two groups | two-propritions z-test  One Population Proportion__Research Question__In previous years, 52% of parents believed that electronics and social media was the cause of their teenager’s lack of sleep. Do more parents today believe that their teenager’s lack of sleep is caused due to electronics and social media?- __Population__: Parents with a teenager (age 13-18)- __Parameter of Interest__: p- __Null Hypothesis__: p = 0.52- __Alternative Hypthosis__: p > 0.52 (note that this is a one-sided test)- __Data__: 1018 people were surveyed. 56% of those who were surveyed believe that their teenager’s lack of sleep is caused due to electronics and social media.<jupyter_code>n = 1018 #number of surveys
pnull = .52 #previous years' percentage of parents. 
phat = .56 #new group of parents surveyed<jupyter_output><empty_output><jupyter_text>Was the increase in responses statistically significant? To test this, we use a one-sided proportions ztest. It is one-sided because we only want to know if the INCREASE is statistically significant - we already know there was an increase. __Use of `proportions_ztest()` from `statsmodels`__:Note the argument `alternative="larger"` indicating a one-sided test. The function returns two values - the z-statistic and the corresponding p-value.<jupyter_code>h_test = sm.stats.proportions_ztest(phat * n, n, pnull, alternative='larger')
print('z-stat: {}'.format(h_test[0]))
print('p-value: {}'.format(h_test[1]))<jupyter_output>z-stat: 2.571067795759113
p-value: 0.005069273865860533<jupyter_text>**Conclusion of the hypothesis test**:  Since the calculated p-value of the z-test is <.05, there is a low probability that, assuming the null hypothesis, the results of this survey increase are randomly by chance. We can thus reject the Null hypothesis that the percentage of parents is the same as previous years' estimate, i.e. 52%.Although, we do not accept the alternate hypothesis, this informally means that there is a good chance of this proportion being more than 52%. ---  Difference in Population Proportions__Research Question__Is there a significant difference between the population proportions of parents of African American children and parents of Hispanic children who report that their child has had some swimming lessons?- __Populations__: All parents of African American children age 6-18 and all parents of Hispanic children age 6-18- __Parameter of Interest__: p1 - p2, where p1 = African American and p2 = hispanic- __Null Hypothesis__: p1 - p2 = 0- __Alternative Hypthesis__: p1 - p2 $\neq$ = 0- __Data__: 247 Parents of African American Children. 36.8% of parents report that their child has had some swimming lessons. 308 Parents of Hispanic Children. 38.9% of parents report that their child has had some swimming lessons.<jupyter_code>n1 = 247 #parents of African American Children
p1 = .36 #percent reporting that their children had swimming lessons

n2 = 308 #parents of hispanic children
p2 = .38 #percentage reporting swimming lessons

#simulate survey data that adheres to the report above
## We do not care about individuals or orders, so this will allow us to cmopare samples
population1 = np.random.binomial(1, p1, n1) 
population2 = np.random.binomial(1, p2, n2) 
population1<jupyter_output><empty_output><jupyter_text>__Use of `ttest_ind()` from `statsmodels`__Difference in population proportion needs t-test. Also, the population follows a binomial distribution here. We can just pass on the two population quantities with the appropriate binomial distribution parameters to the t-test function.The function returns three values: (a) test statisic, (b) p-value of the t-test, and (c) degrees of freedom used in the t-test.<jupyter_code>sm.stats.ttest_ind(population1, population2)[1]<jupyter_output><empty_output><jupyter_text>**Conclusion of the hypothesis test**: Since the p-value is quite high, we cannot reject the Null hypothesis in this case i.e. the difference in the population proportions are not statistically significant. ---  One Population Mean__Research Question__Let's say a cartwheeling competition was organized for some adults. The data looks like following,(80.57, 98.96, 85.28, 83.83, 69.94, 89.59, 91.09, 66.25, 91.21, 82.7 , 73.54, 81.99, 54.01, 82.89, 75.88, 98.32, 107.2 , 85.53, 79.08, 84.3 , 89.32, 86.35, 78.98, 92.26, 87.01)Is there a statistically significant chance that the average cartwheel distance (in inches) for adults more than 80 inches?- __Population__: All adults- __Parameter of Interest__: $\mu$, population mean cartwheel distance.- __Null Hypothesis__: $\mu$ = 80 - __Alternative Hypothesis__: $\mu$ > 80- __Data__:     - 25 adult participants.     - $\mu = 83.84$     - $\sigma = 10.72$<jupyter_code>cwdata = np.array([80.57, 98.96, 85.28, 83.83, 69.94, 89.59, 91.09, 66.25, 91.21, 82.7 , 73.54, 81.99, 54.01, 
                 82.89, 75.88, 98.32, 107.2 , 85.53, 79.08, 84.3 , 89.32, 86.35, 78.98, 92.26, 87.01])
n = len(cwdata) #number of samples
mean = cwdata.mean() #mean
sd = cwdata.std() #standard deviation
(n, mean, sd)
sm.stats.ztest(cwdata, value = 80, alternative = "larger")<jupyter_output><empty_output><jupyter_text>__Conclusion of the Hypothesis Test__Since the p-value (0.0394) is lower than the standard confidence level 0.05, we can reject the Null hypothesis that the mean cartwheel distance for adults (a population quantity) is equal to 80 inches. There is strong evidence in support for the alternatine hypothesis that the mean cartwheel distance is, in fact, higher than 80 inches. Note, we used alternative="larger" in the z-test.We can also plot the histogram of the data to check if it approximately follows a Normal distribution.<jupyter_code>#let's verify with a visualization!
plt.hist(cwdata,bins='scott',edgecolor='k')
plt.show()<jupyter_output><empty_output><jupyter_text>What can you draw from this?  --- Difference in Population MeansResearch Question Considering adults in the NHANES data, do males have a significantly higher mean Body Mass Index than females? Population: Adults in the NHANES data. Parameter of Interest:  𝜇1−𝜇2 , Body Mass Index.Null Hypothesis:  𝜇1=𝜇2 Alternative Hypothesis:  𝜇1≠𝜇2 Data:2976 Females2759 Male Adults<jupyter_code>#import the BMI_NHANES.csv data
from google.colab import files
uploaded = files.upload()
#read into pandas
health = pd.read_csv('BMI_NHANES.csv')
health.info()
#filter the male and female dataframes into two separate objects. Females RIAGENDR = 2, males = 1
females = health[health.RIAGENDR == 2].copy()
males = health[health.RIAGENDR == 1].copy()
#get the number of females
n1 = males.shape[0]
#calculate the mean of BMXBMI
mu1 = males.BMXBMI.mean()
#calculate the standard deviation of BMXBMI
sd1 = males.BMXBMI.std()
#repeat the above cell calculations for men
n2 = females.shape[0]
mu2 = females.BMXBMI.mean()
sd2 = females.BMXBMI.std()
#calcualte if the averages are the same using a comparitor
print(mu1 == mu2)
#show the difference in the means
print(mu1 - mu2)
#now we can see if these are statistically representative 
pval = sm.stats.ztest(females["BMXBMI"].dropna(), males["BMXBMI"].dropna(),alternative='smaller')[1]
round(pval, 2)<jupyter_output><empty_output><jupyter_text>**Conclusion of the hypothesis test**: Since the p-value (6.59e-10) is extremely small, we can reject the Null hypothesis that the mean BMI of males is same as that of females. Note, we used alternative="two-sided" in the z-test because here we are checking for inequality.Create a histogram of the data to check if it approximately follows a Normal distribution.<jupyter_code>#matplotlib code
plt.hist(females["BMXBMI"].dropna(),edgecolor='k',color='pink',bins='scott')
plt.hist(males["BMXBMI"].dropna(),edgecolor='k',color='blue',bins='scott')
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.08855992728217114, 'fraction_numerical': 0.021295935592780157, 'mean_word_length': 4.895381474866038, 'pattern_counts': {'":': 0, '<': 92, '<?xml version=': 0, '>': 97, 'https://': 1, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '21', 'dir': 'jupyter-structured-clean-dedup', 'id': '357674', 'license': 'MIT', 'n_tokens_mistral': 6962, 'n_tokens_neox': 6363, 'n_words': 3012, 'path': '3_statistics.ipynb', 'repo_name': 'burgamacha/DS-Masterclass', 'size': '23115'}