import pandas as pd
import numpy as np
import seaborn as sns
import matplotlib as mat
import matplotlib.pyplot as plt

#1b
refpop=np.array(["UN"]*19726+["AN"]*467) #creates a reference population with the unvaccinated proportions
simpop=[] #empty list to store values
vaxprop=153/6731 #finds the proportion of the vaccinated population
for i in range(10000): #runs simulation 10000 times
    p_repolist=np.random.choice(refpop, 6731) #creates new population out of the original data
    simpop.append(np.sum(p_repolist=="AN")/6731) #finds the proportion of babies born to vaccinated mothers in the sim that had anomolies and adds it to the simpop list
#creates histogram of the simpop list and adds titles
p=sns.displot(data=simpop)
p.set(xlabel="No of babies with anomalies")
p.set(title="proportions of babies with anomalies")
plt.axvline(vaxprop, color="red") #plots the observed value from the raw data

simpop1=np.array(simpop) #creates array out of population proportions
count=np.sum(simpop1<=vaxprop) #divides the amount of times the simulated population produced a result as or more extreme than the observed value, for the lower end
print(2*count) #using 2 tailed p value because we want the values that are as or more extreme than the observed value
p=(2*count)/10000
print("The p value of this data is:", p)

#1c
simpop1.sort() #orders simpop1 array from smallest to largest
lower=simpop1[49] #sets the 49th value in list equal to lower
upper=simpop1[9950] #sets the 950th value in list equal to upper
obs=467/20193 #observed proportion of babies born to unvaccinated mothers with anomalies
lowerCI=(2*obs)-upper #finds cutoff for lower confidence interval
upperCI=(2*obs)-lower #finds cutoff for upper confidence interval
print("lower confidence interval:",lowerCI)#prints the lower confidence interval
print("upper confidence interval:",upperCI) #prints the upper confidence interval
p=sns.displot(data=simpop) #creates histogram of data
p.set(xlabel="No of babies with anomalies") #sets x axis label
p.set(title="proportions of babies with anomalies") #sets title
plt.axvline(vaxprop, color="red") #observed measure
plt.axvline(upperCI, color="green") #upper confidence interval cutoff
plt.axvline(lowerCI, color="green") #lower confidence interval cutoff

#2a
Oatdata=pd.read_csv("cholesterol levels - Sheet1 (1).csv") #imports cholesterol data measurements
Oatdata #prints the table to ensure correctness

def slopePlot(Oatdata, labels=["", ""], line_color="gray", point_color="black"): #"def"- defines function
    import matplotlib.pyplot as plt
    from numpy import array

    OatArr = array(Oatdata) #turn data into an array in order to use plt and array
    fig, ax = plt.subplots(figsize=(4, 3)) #figure size

    x1=0.8
    x2=1.2
    n = OatArr.shape[0] #gives the size of the column
    for i in range(n): #runs a loop
        ax.plot([x1, x2], [OatArr[i,0], OatArr[i,1]], color=line_color) 

        # Plot the points
        ax.scatter(n*[x1-0.01], OatArr[:,0], color=point_color, s=25, label=labels[0])
        ax.scatter(n*[x2+0.01], OatArr[:,1], color=point_color, s=25, label=labels[1])

    # Fix the axes and labels
    ax.set_xticks([x1, x2])
    _ = ax.set_xticklabels(labels, fontsize='x-large')

    return ax

slopePlot(Oatdata, labels=["Corn Flakes", "Oat Bran"], line_color="blue", point_color="orange")
plt.xlabel("Type of Cereal") 
plt.ylabel("LDL Cholesterol Level") 
plt.title("paired cholesterol levels")

#b
p=sns.swarmplot(data=Oatdata)
p.set(ylabel="cholesterol levels")
p.set(title="cholesterol levels observed in different diets of cereal")

##partC Checking Work
#How to calculate the appropriate measure to compare the groups: Median Differences

CornFlakes=Oatdata["CornFlakes"] #create a list of only the cornflakes column from data
OatBran=Oatdata["OatBran"] #create a list of only the oat bran column from data
Corn_arr=np.array(CornFlakes) #turns list of cornflakes into an array for calculations
Oat_arr=np.array(OatBran) #turn list of oat bran into an array for calculations
diff_array=Corn_arr-Oat_arr #subtracts oat bran array from the cornflakes array to get an array of differences
Mobs_array=np.median(diff_array) #find the median difference of the sample
print("The difference in medians between the cornflakes diet effects and the oat bran cereal effects",Mobs_array) #print the difference in medians

#d

#Calculate the p-value for NHST

diff_array=Corn_arr-Oat_arr #subtract oat bran array from the cornflakes array to get an array of differences

Mobs_array=np.median(diff_array) #find the median difference of the sample

onesarray=np.random.choice([-1,1],14) #make an array of randomly chosen +1 and -1 values that will be the same length as the original data

mult_array=onesarray*diff_array #multiply the array of differences and the sign array

newmed_array=np.median(mult_array) #find the median difference for this new array

#do 10000 simulations for NHST
results=np.zeros(10000)
for i in range(10000):
    onesarray = np.random.choice([-1,1],14)
    mult_array= onesarray*diff_array
    newmed_array=np.median(mult_array)
    results[i]=newmed_array

#Making a Histogram   
m=sns.displot(results) #make a histogram of simulation results
m.set(xlabel="Median Differences in Pairs", ylabel="Count")
m.set(title="Difference of Medians of Cholesterol")
plt.axvline(Mobs_array, color="red") #plot where the observed effect size is (Mobs)
plt.axvline((-Mobs_array), color="red") #plot the the lower limit of the effect size

#p-value calculations
var=np.array(results)
high_extremes=np.sum(var>=Mobs_array) #find results that are as or more extreme for upper bound

low_extremes=np.sum(var<=(-Mobs_array))#find resullts that are as or more extreme for lower bound

pval=(high_extremes+low_extremes)/10000 #calculate p-value
print("The p-value we found from NHST is",pval) #print p-value

#3d
pvals=[0.004,0.017,0.021,0.04,0.23,0.32] #import p values
import statsmodels.stats.multitest as smm
smm.multipletests(pvals, alpha=0.05, method='fdr_bh') #preform benjamini hochberg analysis