# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# ## WGP Followup - Variability
# The purpose of this lab will be to continue on with the follow-up rounds of the Water Guard Promotion study.  We will also continue look at the Water Quality variables because they are continuous variables. 
#
# We should program the following
# 1. Graphing Data ?  Whisker Plot?
# 2. Variance - Standard Error - Z score
# 3. Confidence Interval
# 4. Hypothesis test ( and its ensuing CI etc) 
#
#
#

# +
# run this cell - obligatory imports in every file
from datascience import *
import pandas as pd
from pandas import read_stata
import numpy as np
# %matplotlib inline
import matplotlib.pyplot as plt

plt.style.use('fivethirtyeight')
from pylab import plot, show, savefig, xlim, figure, \
                hold, ylim, legend, boxplot, setp, axes
# -

# read in stata data set using pandas ( pd) 
#pd.read_stata('stata.dta')
#pd.read_stata('bwm_child_EVDvars.dta')
WGP3rds_df = pd.DataFrame(pd.read_stata('WGP_3waves_Data8.dta'))
WGP3rds_table= Table.from_df(WGP3rds_df)
#WGP_df.head(10)

# ##### This is a large dataset, basically three datasets merged together, one for baseline, one for short term follow up and one for long term followup
#  - Round = 1 : baseline
#  - Round = 2 : 3 week followup
#  - Round = 3 : 3 month followup
#  
# For many of the variables they are only asked in one of the three rounds

# **From Lab 4:**
#  - The variable for self_reported chlorine use was "c6n" in Round 2, and 'c5n' in Round 3
#  - The variable for chlorine use is "c12n21pnk" in Round 2 and 'c15npt2or1pnk' in Round3
#  
#  
# ** These variables have been combined across rounds for the ease of programming **
#  - 'Selfrptpct' is self reported chlorine use in both round 2 and round 3
#  - 'Vldclpct'  is validated chlorine use in both rounds
#  
#  

# +
# Here is a list of all of the possible categories / columns
#list(WGP3rds_df)
# -

#WGP_baseline.groupby('e1_num_kids_under_5').size()
WGP3rds_df.groupby('treatment_arm').size()

WGP3rds_df.groupby('round').size()

#  Here is a new way to look at the data - grouping by two variables in Pandas
WGP3rds_df.groupby([ 'round', 'treatment_arm',]).size()

# Which brings us to a new  way to look at the outcome data - grouping by two variables in Pandas an calling for means
# Make a smaller df with just the outcome variables
WGP_3rds_outcomesonly= pd.DataFrame(WGP3rds_df,columns=['round','treatment_arm','Selfrptpct','Vldclpct'])
# Then call for the means
WGP_3rds_outcomesonly.groupby([ 'round','treatment_arm']).mean()

WGPRd2 = WGP3rds_table.where("round",2).select("a1_cmpd_id","treatment_arm",
                                           "c6_current_water_treated_wg", 
                                           'c6_curr_water_treat_other_c',
                                           'c12_chlorine_meter_reading',
                                           'c11_chlorine_color','c12n21pnk', 'c6n'
                                          )
WGPRd2

# make a pandas dataframe df for Round 2
WGP_rnd2df= WGP3rds_df.where(WGP3rds_df['round'] ==2)

# we can start by looking at the estimate of WG use in Round 2 across all treatment arms ( for the entire sample) 
np.mean(WGPRd2.column('c12n21pnk'))

# this is the numpy stadard error
np.std(WGPRd2.column('c12n21pnk'))

# +

# - but we have to correct in this by the square root of the number of observations
# which happens in thecase of a proportion of a categorical variable 
np.std(WGPRd2.column('c12n21pnk'))/WGPRd2.num_rows**0.5
# -

# Here is the histogram to remind us that our categorical data doesnt look like a bell curve
WGPRd2.hist('c12n21pnk')


# Compute Standard Error of a proportion
# We can do this the first time for the entire sample
# Save the followng values
p_all=np.mean(WGPRd2.column('c12n21pnk'))
N_all=WGPRd2.num_rows

#  The formula for the SE of a proportion is sqrt((p1(1-p1))/N)
se_all=((p_all*(1-p_all))/N_all)**0.5
se_all

# Upper bound of CI
# Using Z = 1.96
Z=1.96
upper_all = p_all+Z*se_all
lower_all = p_all-Z*se_all
CI_all=(upper_all,lower_all)
CI_all

# ** Using this technique we can look at the Confidence Intervals for each of the treatment arms **
#



# Here is the code for just arm 1
p_1=np.mean(WGPRd2.where("treatment_arm",1).column('c12n21pnk'))
N_1=WGPRd2.where("treatment_arm",1).num_rows
N_1
print("In Round 2, the number of households in arm 1 Control is:", N_1)
print("The share of households with Validated WG use is: ", p_1)

se_1=((p_1*(1-p_1))/N_1)**0.5
print("The standard error of the Validated WG estimate is: ", se_1)

# Upper bound of CI
# Using Z = 1.96
Z=1.96
upper_1 = p_1+Z*se_1
lower_1 = p_1-Z*se_1
CI_1=(upper_1,lower_1)
print("The 95% Confidence Interval for Validated WG use is: ", CI_1)

# ** Repeating for Arm 2 **
#

# Here is the code for just arm 2
p_2=np.mean(WGPRd2.where("treatment_arm",2).column('c12n21pnk'))
N_2=WGPRd2.where("treatment_arm",2).num_rows
N_2
print("In Round 2, the number of households in arm 1 Control is:", N_2)
print("The share of households with Validated WG use is: ", p_2)
se_2=((p_2*(1-p_2))/N_2)**0.5
print("The standard error of the Validated WG estimate is: ", se_2)
# Upper bound of CI
# Using Z = 1.96
Z=1.96
upper_2 = p_2+Z*se_2
lower_2 = p_2-Z*se_2
CI_2=(upper_2,lower_2)
print("The 95% Confidence Interval for Validated WG use is: ", CI_2)


# ## Question 1.1 - What can we tell by comparing the confidence intervals for Arm 1 and Arm 2.  Do they overlap?  What does that mean?

# ***put your answer here***

# **Hypothesis test**

# +
#  Lets test the difference between arm1 and arm 2

# The Null Hypotetheis is that the means are the same
# So we can test if the difference in means is not =0
Difference_1to2 = p_1-p_2
print ( "The difference between Arm 1 and Arm 2  is",Difference_1to2)



# +
# The formula for this SE is Sqrt( SE_1 + SE_2)
se_diff_1to2 = (se_1+se_2)**0.5
CI_diff_1to2 = (Difference_1to2+1.96*se_diff_1to2, Difference_1to2-1.96*se_diff_1to2)
print("The 95% Confidence Interval for Validated WG use is: ", CI_diff_1to2)


# -

# Lets look at the test startistic for the Null Hypothesis
Z_diff_1to2 = Difference_1to2/se_diff_1to2
print("The Test Statistic for the difference between Validated WG use in Arm 1 and Arm2 is: ",Z_diff_1to2)

# ## Question 1.2 - What can we tell by looking at the confidence interval for the difference?  Do they overlap with zero?  Should we reject the null hypothesis? What does this mean?

# ***put your answer here***

# ## Question 1.3 - Propose an additional hypothesis to test with this data?  Program the hypothesis test, difference, and confidence interval

# ***written answer proposing a hypothesis test ***

# +
# coding section

# +
# coding section
# -

# ***written answer discussing outcome ***

# ## Working with Springs Water Quality Data again

# read in stata data set using pandas ( pd) to dataframe
# This is the dataset for water quality testing at springs
# From the BiWeekly Monitoring study 
# The most important variables will be trt_07, and bwm_round and ecmpn1 and tcmpn1
BWM_WQ_s = pd.DataFrame(pd.read_stata('WQ_sp_forData8.dta'))
# make a datascience table
BWM_WQ_spring = Table.from_df(BWM_WQ_s)
# make log10 transforms of E coli numbers ( adding columns)
BWM_WQ_spring = BWM_WQ_spring.with_column('log10ecmpn1',BWM_WQ_spring.apply(np.log10,'ecmpn1'))
BWM_WQ_spring = BWM_WQ_spring.with_column('log10tcmpn1',BWM_WQ_spring.apply(np.log10,'tcmpn1'))
BWM_WQ_spring = BWM_WQ_spring.drop("dd1","mm1","yy1","chlorine")
BWM_WQ_spring


# This is the mean of the Water Quality Before
np.mean(BWM_WQ_spring.where("bwm_round",are.below(6)).column('log10ecmpn1'))

# This is the mean of Water Quality After
np.mean(BWM_WQ_spring.where("bwm_round",are.above(5)).column('log10ecmpn1'))

# A test of whether water quality changed 
difference_allsprings = np.mean(BWM_WQ_spring.where("bwm_round",are.below(6)).column('log10ecmpn1'))-np.mean(BWM_WQ_spring.where("bwm_round",are.above(5)).column('log10ecmpn1'))
print ("the difference is", difference_allsprings)
SE_allsprings = np.std(BWM_WQ_spring.column('log10ecmpn1'))
print ("the SE is", SE_allsprings)
Z_allsprings= difference_allsprings /SE_allsprings
print ("the Z score is", Z_allsprings)

# ## But what kind of Test are we really interested in?
#

# Mean of WQ at treated springs after
np.mean(BWM_WQ_spring.where("bwm_round",are.above(5)).where('trt_07',1).column('log10ecmpn1'))

# Mean of of WQ at control  springs after
np.mean(BWM_WQ_spring.where("bwm_round",are.above(5)).where('trt_07',2).column('log10ecmpn1'))

np.std (BWM_WQ_spring.where("bwm_round",are.above(5)).column('log10ecmpn1'))

# ## Question 2 - can you construct a hypothesis test of whether the Water Quality is cleaner at protected springs after the spring protection?
#

# +
# Coding section

# +
# Coding section
# -



## lets look at the data graphically 
BWM_WQ_spring.hist("log10ecmpn1")

#  We can also look at Four relevant Histograms
BWM_WQ_spring.where('bwm_round',are.below(6)).where('trt_07',1).hist("log10ecmpn1")
BWM_WQ_spring.where('bwm_round',are.below(6)).where('trt_07',2).hist("log10ecmpn1")
BWM_WQ_spring.where('bwm_round',are.above(5)).where('trt_07',1).hist("log10ecmpn1")
BWM_WQ_spring.where('bwm_round',are.above(5)).where('trt_07',2).hist("log10ecmpn1")

# ## Question 3.1 - can you discuss how our hypothesis tests relate to these 4 histograms.  Which ones are we comparing and why?

# +
# Here we can do a boxplot graph to compare the Treatment Springs before and After
# A boxplot is a common graph that describes the data by quartiles - 25% of the data
# The box in the middle displays the middle two quartiles and the red line is the median

before_treatment = BWM_WQ_spring.where('bwm_round',are.below(6)).where('trt_07',1).column('log10ecmpn1')
after_treatment = BWM_WQ_spring.where('bwm_round',are.above(5)).where('trt_07',1).column('log10ecmpn1')

# -

data = [before_treatment, after_treatment]
fig, ax = plt.subplots()
ax.set_title('Comparing ECMPN Before and After ')
ax.boxplot(data)
plt.show()

# ## Question 3 - can you redo this boxplot graph to compare the Treatment Springs after Treatment to the Control Springs after Treatment

# +
#code here for subsetting data

# +
# code here for new boxplot

StarcoderdataJupyterScriptsDedupFiltered