{
"cells": [
{
"cell_type": "markdown",
"id": "a45be4b8",
"metadata": {},
"source": [
"# Chapter 12: Modeling categorical relationships"
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "2c107588",
"metadata": {},
"outputs": [],
"source": [
"import pandas as pd\n",
"import sidetable\n",
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"import seaborn as sns\n",
"from scipy.stats import norm, t, binom, chi2\n",
"import pingouin as pg\n",
"import matplotlib\n",
"\n",
"import rpy2.robjects as ro\n",
"from rpy2.robjects.packages import importr\n",
"from rpy2.robjects import pandas2ri\n",
"pandas2ri.activate()\n",
"from rpy2.robjects.conversion import localconverter\n",
"%load_ext rpy2.ipython\n",
"\n",
"# import NHANES package\n",
"base = importr('NHANES')\n",
"\n",
"with localconverter(ro.default_converter + pandas2ri.converter):\n",
" NHANES = ro.conversion.rpy2py(ro.r['NHANES'])\n",
"\n",
" \n",
"NHANES = NHANES.drop_duplicates(subset='ID')\n",
"NHANES_adult = NHANES.dropna(subset=['Weight']).query('Age > 17 and BPSysAve > 0')\n",
"\n",
"rng = np.random.default_rng(123456)\n"
]
},
{
"cell_type": "markdown",
"id": "418f71a7",
"metadata": {},
"source": [
"## Table 12.1\n"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "0aee7fd3",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" | \n",
" Candy Type | \n",
" Count | \n",
" nullExpectation | \n",
" squaredDifference | \n",
"
\n",
" \n",
" \n",
" \n",
" | 0 | \n",
" chocolate | \n",
" 30 | \n",
" 33.333333 | \n",
" 11.111111 | \n",
"
\n",
" \n",
" | 1 | \n",
" licorice | \n",
" 33 | \n",
" 33.333333 | \n",
" 0.111111 | \n",
"
\n",
" \n",
" | 2 | \n",
" gumball | \n",
" 37 | \n",
" 33.333333 | \n",
" 13.444444 | \n",
"
\n",
" \n",
"
\n",
"
"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"x = np.arange(0.01, 20, .01)\n",
"dfvals = [1, 2, 4, 8]\n",
"\n",
"x_df = [(i, j) for i in x for j in dfvals]\n",
"\n",
"chisqDf = pd.DataFrame(x_df, columns = ['x', 'df'])\n",
"chisqDf['chisq'] = chi2.pdf(chisqDf.x, chisqDf.df)\n",
"for df in dfvals:\n",
" chisqDf.loc[chisqDf.df == df, 'chisq'] = chisqDf[chisqDf.df == df]['chisq'] / chisqDf[chisqDf.df == df]['chisq'].max()\n",
"fig, ax = plt.subplots(1, 2, figsize=(12,6))\n",
"\n",
"sns.lineplot(data=chisqDf, x='x', y='chisq', style='df', ax=ax[0])\n",
"ax[0].set_xlabel(\"Chi-squared values\")\n",
"ax[0].set_ylabel('Density')\n",
"\n",
"\n",
"# simulate 50,000 sums of 8 standard normal random variables and compare\n",
"# to theoretical chi-squared distribution\n",
"\n",
"# create a matrix with 50k columns of 8 rows of squared normal random variables\n",
"dSum = (rng.normal(size=(50000, 8)) ** 2).sum(axis=1)\n",
"\n",
"sns.histplot(dSum, ax=ax[1], stat='density')\n",
"ax[1].set_ylabel(\"Density\")\n",
"ax[1].set_xlabel(\"Sum of squared random normal variables\")\n",
"\n",
"csDf = pd.DataFrame({'x': np.arange(0.01, 30, 0.01)})\n",
"csDf['chisq'] = chi2.pdf(csDf.x, 8)\n",
"_ =sns.lineplot(data=csDf, x='x', y='chisq', ax=ax[1], color='black')\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "7df156ae",
"metadata": {},
"source": [
"## Table 12.2"
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "271a20c6",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"search_conducted False True All\n",
"driver_race \n",
"Black 36244 1219 37463\n",
"White 239241 3108 242349\n",
"All 275485 4327 279812\n",
"\n",
"Normalized:\n",
"search_conducted False True All\n",
"driver_race \n",
"Black 0.129530 0.004356 0.133886\n",
"White 0.855006 0.011107 0.866114\n",
"All 0.984536 0.015464 1.000000\n"
]
}
],
"source": [
"stopData = pd.read_csv('https://raw.githubusercontent.com/statsthinking21/statsthinking21-figures-data/main/CT_data_cleaned.csv')\n",
"table = pd.crosstab(stopData['driver_race'], stopData['search_conducted'], margins=True)\n",
"print(table)\n",
"n = stopData.shape[0]\n",
"print('')\n",
"print('Normalized:')\n",
"table_normalized = pd.crosstab(stopData['driver_race'], stopData['search_conducted'], margins=True, normalize='all')\n",
"print(table_normalized)"
]
},
{
"cell_type": "markdown",
"id": "6baec10f",
"metadata": {},
"source": [
"## Chi-squared test result"
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "4f372548",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"1.9001204194035138e-182"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"expected = np.outer(table_normalized.values[:2, 2], table_normalized.values[2, :2]) * n\n",
"actual = pd.crosstab(stopData['driver_race'], stopData['search_conducted'], margins=False)\n",
"diff = expected - actual\n",
"stdSqDiff = diff **2 / expected\n",
"chisq = stdSqDiff.sum().sum()\n",
"pval = chi2.pdf(chisq, 1)\n",
"pval"
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "4be02b50",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" | \n",
" test | \n",
" lambda | \n",
" chi2 | \n",
" dof | \n",
" pval | \n",
" cramer | \n",
" power | \n",
"
\n",
" \n",
" \n",
" \n",
" | 0 | \n",
" pearson | \n",
" 1.0 | \n",
" 827.005515 | \n",
" 1.0 | \n",
" 7.255988e-182 | \n",
" 0.054365 | \n",
" 1.0 | \n",
"
\n",
" \n",
"
\n",
"
\n",
"\n",
"
\n",
" \n",
" \n",
" | search_conducted | \n",
" False | \n",
" True | \n",
"
\n",
" \n",
" | driver_race | \n",
" | \n",
" | \n",
"
\n",
" \n",
" \n",
" \n",
" | Black | \n",
" 3.330746 | \n",
" -26.576456 | \n",
"
\n",
" \n",
" | White | \n",
" -1.309550 | \n",
" 10.449072 | \n",
"
\n",
" \n",
"
\n",
"