{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Descriptive Statistics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#importing the basic libraries\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#creating a dataframe\n",
    "df = pd.DataFrame([\n",
    " { 'name': 'Ahmed', 'salary': 1200, 'age': 22} ,\n",
    " { 'name': 'Saleh', 'salary': 22000, 'age': 33},\n",
    " { 'name': 'Hashir', 'salary': 25000, 'age': 44 },\n",
    " { 'name': 'Qasim', 'salary': 35000, 'age': 23 },\n",
    " { 'name': 'Junaid', 'salary': 25000, 'age': 42 },\n",
    " { 'name': 'Ammaar', 'salary': 15000, 'age': 32 },\n",
    " { 'name': 'Bill Gates', 'salary': 300000, 'age': 52 }    \n",
    "])\n",
    "\n",
    "#print it\n",
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#another way of creating a data frame\n",
    "df2 = pd.DataFrame({'Name':['Taha', 'Mahad'],\n",
    "                     'Salary':[12000,14000]})\n",
    "\n",
    "#print it\n",
    "df2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Central Tendency (Mean)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#calculate the mean\n",
    "print(\"Mean: \",df['salary'].mean())\n",
    "\n",
    "#convert to integer\n",
    "print(\"Integerized Mean: \", int(df['salary'].mean()))\n",
    "\n",
    "#in the standard format\n",
    "print(\"Standardized Mean: \",\"{:,}\".format(int(df['salary'].mean())))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Central Tendency (Median)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Mean is on the higher side due to Bill Gates Salary (an outlier)\n",
    "#Let's calculate the Median (central value)\n",
    "print(df['salary'].sort_values())\n",
    "print(df.salary.sort_values())\n",
    "\n",
    "print(\"Median: \", df['salary'].median())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Central Tendency (Mode)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Let's build a data set\n",
    "reviews_df = pd.DataFrame([\n",
    " { 'restaurant': 'Student Bryani', 'reviewer': 'Ahmed', 'stars': 2 } ,\n",
    " { 'restaurant': 'Student Bryani', 'reviewer': 'Saleh', 'stars': 2 },\n",
    " { 'restaurant': 'Student Bryani', 'reviewer': 'Hashir', 'stars': 4 },\n",
    " { 'restaurant': 'Student Bryani', 'reviewer': 'Qasim', 'stars': 4 },\n",
    " { 'restaurant': 'Student Bryani', 'reviewer': 'Junaid', 'stars': 4 },\n",
    " { 'restaurant': 'Student Bryani', 'reviewer': 'Ammaar', 'stars': 3 },\n",
    " { 'restaurant': 'Student Bryani', 'reviewer': 'Bill Gates', 'stars': 2 }    \n",
    "])\n",
    "reviews_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "reviews_df.columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "reviews_df.dtypes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#calling the mode function\n",
    "reviews_df['stars'].mode()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import HTML display from Amazon\n",
    "from IPython.display import display, HTML\n",
    "display(HTML('''<img src=\"data:image/png;base64,[encoded data, 135512 characters]\">'''))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#checking the mode\n",
    "display(HTML('''<img src=\"data:image/png;base64,[encoded data, 48812 characters]\">'''))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.tail()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Make it so we only have three decimal places\n",
    "pd.set_option('display.float_format', lambda x: '%.3f' % x)\n",
    "df.describe()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Box and Whiskers Plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#df.boxplot(column='salary', sym='o', return_type='axes')\n",
    "df.boxplot(column='salary', sym='o', whis = 'range', return_type='axes')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Let's try to remove Bill gates and see the plot\n",
    "mask = df['name'].isin(['Bill Gates'])\n",
    "df1 = df[~mask]\n",
    "\n",
    "df1.boxplot(column='salary', sym='o', whis = 'range', return_type='axes')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#if we want to set whiskers at max and min (to identify the outlier range)\n",
    "df1.boxplot(column='salary', sym='o', whis = 'range', return_type='axes')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Scatter Plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Create data\n",
    "N = 20\n",
    "x = np.random.rand(N)\n",
    "y = np.random.rand(N)\n",
    "print(x)\n",
    "print(y)\n",
    "colors = [[1,0,0]] #this means black. try out different values of RGB between 0 and 1\n",
    "area = np.pi*10 #this is the size of each dot. try out different values\n",
    "\n",
    "# Plot\n",
    "plt.scatter(x, y, s = area, c=colors, alpha=0.9) #alpha is the color intensity. try out different values\n",
    "plt.title('Scatter Plot Example')\n",
    "plt.xlabel('x label')\n",
    "plt.ylabel('y label')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.scatter(df1.age, df1.salary, s = area, c=colors, alpha=0.9) #alpha is the color intensity. try out different values"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Variance and Standard Deviation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "list_one = pd.Series([-10, 0, 10, 20, 30])\n",
    "list_two = pd.Series([8, 9, 10, 11, 12])\n",
    "print(\"list one is\")\n",
    "print(list_one)\n",
    "print(\"list two is\")\n",
    "print(list_two)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"The mean of list_one is\", list_one.mean())\n",
    "print(\"The mean of list_two is\", list_two.mean())\n",
    "print(\"The median of list_one is\", list_one.median())\n",
    "print(\"The median of list_two is\", list_two.median())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ddof=0 has to do with sample variance\n",
    "print(list_one.var(ddof=1))\n",
    "print(list_two.var(ddof=1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "# Data points [-10, 0, 10, 20, 30]\n",
    "# Variance: 200\n",
    "sd1 = math.sqrt(250)\n",
    "sd1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Data points [8, 9, 10, 11, 12]\n",
    "# Variance: 2.0\n",
    "sd2 = math.sqrt(2)\n",
    "sd2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#define a class to plot the normal distribution curve\n",
    "class norm1:\n",
    "    #constructor (mean, sd, density measure)\n",
    "    def __init__(self, a1, b1, c1):\n",
    "        self.a1 = a1\n",
    "        self.b1 = b1\n",
    "        self.c1 = c1\n",
    "    \n",
    "    #normal curve plotting function\n",
    "    def dist_curve(self):\n",
    "        plt.plot(self.c1, 1/(self.b1 * np.sqrt(2 * np.pi)) *\n",
    "            np.exp( - (self.c1 - self.a1)**2 / (2 * self.b1**2) ), linewidth=2, color='y')\n",
    "        plt.show()\n",
    "\n",
    "c = np.random.normal(list_one.mean(), sd1, 1000)\n",
    "w1, x1, z1 = plt.hist(c, 60, density=True) #hist\n",
    "print(x1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hist1 = norm1(list_one.mean(), sd1, x1)\n",
    "plot1 = hist1.dist_curve()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#let's see the salary plots\n",
    "varsal = df['salary'].var(ddof=1)\n",
    "sdsal = math.sqrt(varsal)\n",
    "\n",
    "print(sdsal)\n",
    "print(df['salary'].mean())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.salary.hist()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df1.salary.hist()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Outlier Analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rich_df = pd.read_excel(\"richpeople.xlsx\")\n",
    "rich_df.columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rich_df.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rich_df = rich_df[rich_df['year'] == 2014]\n",
    "rich_df.head(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rich_df.isna().sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#remove missing values and then describe\n",
    "rich_df = rich_df.dropna(subset=['age', 'networthusbillion', 'foundingdate'])\n",
    "rich_df.describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rich_df.boxplot(column='networthusbillion', whis='range', return_type='axes')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rich_df.boxplot(column='age', whis='range', return_type='axes')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rich_df['citizenship'].value_counts()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rich_df.sort_values(by='age', ascending=False).head(3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#some of the old people are rich.. let's check their standard deviation of ages\n",
    "rich_df['agestd'] = ((rich_df['age'] - rich_df['age'].mean()).apply(abs) / rich_df['age'].std())\n",
    "rich_df.sort_values(by='age', ascending=False).head(15)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#wealth standard dev\n",
    "rich_df['wealth_stdev'] = ((rich_df['networthusbillion'] - rich_df['networthusbillion'].mean()).apply(abs) / rich_df['networthusbillion'].std())\n",
    "rich_df.sort_values(by='networthusbillion', ascending=False).head(10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rich_df2 = rich_df['networthusbillion'].apply(lambda x: x>50)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#plot the histogram to find out\n",
    "rich_df['networthusbillion'].hist()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#some more bins\n",
    "rich_df['networthusbillion'].hist(bins=25)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rich_df['networthusbillion'].describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rich_df['age'].hist(bins=30)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Inferential Statistics"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## T-Test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy import stats\n",
    "\n",
    "## Define 2 random distributions\n",
    "#Sample Size\n",
    "N = 10\n",
    "\n",
    "#Gaussian distributed data with mean = 2 and var = 1\n",
    "a = np.random.randn(N)+2\n",
    "print(a)\n",
    "\n",
    "#Gaussian distributed data with with mean = 0 and var = 1\n",
    "b = np.random.randn(N)\n",
    "print(b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Calculate the Standard Deviation\n",
    "#Calculate the variance to get the standard deviation\n",
    "\n",
    "#For unbiased max likelihood estimate we have to divide the var by N-1, and therefore the parameter ddof = 1\n",
    "var_a = a.var(ddof=1)\n",
    "var_b = b.var(ddof=1)\n",
    "\n",
    "#std deviation\n",
    "s = np.sqrt((var_a + var_b)/2)\n",
    "\n",
    "## Calculate the t-statistics\n",
    "t = (a.mean() - b.mean())/(s*np.sqrt(2/N))\n",
    "\n",
    "## Compare with the critical t-value\n",
    "#Degrees of freedom\n",
    "df = 2*N - 2\n",
    "\n",
    "#p-value after comparison with the t \n",
    "p = 1 - stats.t.cdf(t,df=df)\n",
    "\n",
    "print(\"t = \" + str(t))\n",
    "print(\"p = \" + str(2*p))\n",
    "### You can see that after comparing the t statistic with the critical t value (computed internally) we get a good p value of 0.0005 and thus we reject the null hypothesis and thus it proves that the mean of the two distributions are different and statistically significant.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Cross Checking with the internal scipy function\n",
    "t2, p2 = stats.ttest_ind(a,b)\n",
    "print(\"t = \" + str(t2))\n",
    "print(\"p = \" + str(p2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1-Way ANOVA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd \n",
    "\n",
    "# load data file\n",
    "d = pd.read_csv(\"https://reneshbedre.github.io/myfiles/anova/onewayanova.txt\", sep=\"\\t\")\n",
    "print(d)\n",
    "# generate a boxplot to see the data distribution by treatments. Using boxplot, we can easily detect the differences \n",
    "# between different treatments\n",
    "d.boxplot(column=['A', 'B', 'C', 'D'], grid=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load packages\n",
    "import scipy.stats as stats\n",
    "# stats f_oneway functions takes the groups as input and returns F and P-value\n",
    "fvalue, pvalue = stats.f_oneway(d['A'], d['B'], d['C'], d['D'])\n",
    "print(fvalue, pvalue)\n",
    "# 17.492810457516338 2.639241146210922e-05\n",
    "\n",
    "# get ANOVA table as R like output\n",
    "import statsmodels.api as sm\n",
    "from statsmodels.formula.api import ols\n",
    "# reshape the d dataframe suitable for statsmodels package \n",
    "d_melt = pd.melt(d.reset_index(), id_vars=['index'], value_vars=['A', 'B', 'C', 'D'])\n",
    "# replace column names\n",
    "d_melt.columns = ['index', 'treatments', 'value']\n",
    "# Ordinary Least Squares (OLS) model\n",
    "model = ols('value ~ C(treatments)', data=d_melt).fit()\n",
    "anova_table = sm.stats.anova_lm(model, typ=2)\n",
    "anova_table"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load packages\n",
    "from statsmodels.stats.multicomp import pairwise_tukeyhsd\n",
    "# perform multiple pairwise comparison (Tukey HSD)\n",
    "m_comp = pairwise_tukeyhsd(endog=d_melt['value'], groups=d_melt['treatments'], alpha=0.05)\n",
    "print(m_comp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Assumption 1: Homogeneity of variances (variances are equal between treatment groups) (Levene or Bartlett Test)\n",
    "import scipy.stats as stats\n",
    "w, pvalue = stats.levene(d['A'], d['B'], d['C'], d['D'])\n",
    "print(w, pvalue)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Assumption 2: Residuals (experimental error) are normally distributed (Shapiro Wilks Test)\n",
    "import scipy.stats as stats\n",
    "w, pvalue = stats.shapiro(model.resid)\n",
    "print(w, pvalue)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2-Way ANOVA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load packages\n",
    "import pandas as pd\n",
    "import seaborn as sns\n",
    "# load data file\n",
    "d = pd.read_csv(\"https://reneshbedre.github.io/myfiles/anova/twowayanova.txt\", sep=\"\\t\")\n",
    "print(d)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# reshape the d dataframe suitable for statsmodels package \n",
    "# you do not need to reshape if your data is already in stacked format. Compare d and d_melt tables for detail \n",
    "# understanding \n",
    "d_melt = pd.melt(d, id_vars=['Genotype'], value_vars=['1_year', '2_year', '3_year'])\n",
    "# replace column names\n",
    "d_melt.columns = ['Genotype', 'years', 'value']\n",
    "# generate a boxplot to see the data distribution by genotypes and years. Using boxplot, we can easily detect the \n",
    "# differences between different groups\n",
    "sns.boxplot(x=\"Genotype\", y=\"value\", hue=\"years\", data=d_melt, palette=\"Set3\") "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load packages\n",
    "import statsmodels.api as sm\n",
    "from statsmodels.formula.api import ols\n",
    "# Ordinary Least Squares (OLS) model\n",
    "# C(Genotype):C(years) represent interaction term\n",
    "model = ols('value ~ C(Genotype) + C(years) + C(Genotype):C(years)', data=d_melt).fit()\n",
    "anova_table = sm.stats.anova_lm(model, typ=2)\n",
    "anova_table"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
