import pandas as pd
import numpy as np
from matplotlib.pyplot import subplots

data = pd.read_csv('family_heights.csv')
data.head(3)

	family	father	mother	midparentHeight	children	childNum	sex	childHeight
0	1	78.5	67.0	75.43	4	1	male	73.2
1	1	78.5	67.0	75.43	4	2	female	69.2
2	1	78.5	67.0	75.43	4	3	female	69.0

coefficients = np.polyfit(data['midparentHeight'], data['childHeight'], 1)
coefficients

array([ 0.6373609 , 22.63624055])

data1 = data.assign(regressionEstimate = coefficients[0]*data['midparentHeight'] + coefficients[1])
data1.head(10)

	family	father	mother	midparentHeight	children	childNum	sex	childHeight	regressionEstimate
0	1	78.5	67.0	75.43	4	1	male	73.2	70.712373
1	1	78.5	67.0	75.43	4	2	female	69.2	70.712373
2	1	78.5	67.0	75.43	4	3	female	69.0	70.712373
3	1	78.5	67.0	75.43	4	4	female	69.0	70.712373
4	2	75.5	66.5	73.66	4	1	male	73.5	69.584244
5	2	75.5	66.5	73.66	4	2	male	72.5	69.584244
6	2	75.5	66.5	73.66	4	3	female	65.5	69.584244
7	2	75.5	66.5	73.66	4	4	female	65.5	69.584244
8	3	75.0	64.0	72.06	2	1	male	71.0	68.564467
9	3	75.0	64.0	72.06	2	2	female	68.0	68.564467

data2 = data1.assign(residuals = data1['childHeight'] - data1['regressionEstimate'])
data2.head(5)

	family	father	mother	midparentHeight	children	childNum	sex	childHeight	regressionEstimate	residuals
0	1	78.5	67.0	75.43	4	1	male	73.2	70.712373	2.487627
1	1	78.5	67.0	75.43	4	2	female	69.2	70.712373	-1.512373
2	1	78.5	67.0	75.43	4	3	female	69.0	70.712373	-1.712373
3	1	78.5	67.0	75.43	4	4	female	69.0	70.712373	-1.712373
4	2	75.5	66.5	73.66	4	1	male	73.5	69.584244	3.915756

fig, ax = subplots(figsize = (5, 5))
ax.scatter(data2['midparentHeight'], data2['residuals'])
ax.axhline(y=0, color='red', linewidth = 3);

# We showed in class that the mean of the residuals equals 0. 
# Let us check this for the heights data
print('Mean of the residuals = ', data2['residuals'].mean())

Mean of the residuals =  -7.323002542212767e-14

Due to precision errors, the value does not come out to be exactly 0, but it is really close!

# We showed in class that the correlation between the predictor and the residuals is 0
# Let us check for the heights data
print('Correlation of midparentHeights and residuals = ', data2['midparentHeight'].corr(data2['residuals']))

Correlation of midparentHeights and residuals =  -1.3549018592547941e-15

Again, due to precision errors, the value does not come out to be exactly 0, but it is really close!

data2.shape

(934, 10)

# Let's study Birth Weight versus Gestational Days
baby = pd.read_csv('baby.csv')
fig, ax = subplots(figsize = (5, 5))
ax.scatter(baby['Gestational Days'], baby['Birth Weight'])
ax.set_xlabel('Gestational Days')
ax.set_ylabel('Birth Weight')

Text(0, 0.5, 'Birth Weight')

# Let's calculate the correlation
r = baby['Gestational Days'].corr(baby['Birth Weight'])
r

np.float64(0.40754279338885113)

# In this section, we use bootstrap to compute the 95% confidence interval on the true slope
# Note that the bootstrap works only under the regression model given in the class slides
num_repetitions = 5000

slopes = []
for i in range(num_repetitions):
    sample = baby.sample(frac=1, replace = True)
    coefficients = np.polyfit(sample['Gestational Days'], sample['Birth Weight'], 1)
    slopes.append(coefficients[0])

#plot the histogram of slopes and mark the 95% confidence interval
slopes = np.array(slopes)
left = np.percentile(slopes, 2.5)
right = np.percentile(slopes, 97.5)
fig, ax = subplots(figsize = (5, 5))
ax.hist(slopes, density = True, bins=50)
ax.hlines(0, left, right, color = 'yellow', linewidth = 7)
print('The 95% confidence interval on the true slope is: (', left, ',', right, ')')

The 95% confidence interval on the true slope is: ( 0.3845205020430253 , 0.5556891240274472 )

# We can use 95% confidence interval idea to test the hypothesis that there is no linear association between two variables
# Let's first check Birth weight versus maternal age which does not seem to be correlated
fig, ax = subplots(figsize = (5, 5))
ax.scatter(baby['Maternal Age'], baby['Birth Weight'])

print('Correlation = ', baby['Maternal Age'].corr(baby['Birth Weight']))

Correlation =  0.0269829110029295

We perform the following hypothesis test using 95% confidence intervals on slope:

• Null hypothesis: The slope of the true line is 0
• Alternate hypothesis: The slope of the true line is not 0

# Estimating 95% confidence interval on the true slope using bootstrapping
num_repetitions = 5000

slopes = []
for i in range(num_repetitions):
    sample = baby.sample(frac = 1, replace=True)
    coefficients = np.polyfit(sample['Maternal Age'], sample['Birth Weight'], 1)
    slopes.append(coefficients[0])

# plot a histogram of the slopes
#calculate the percentiles
slopes = np.array(slopes)
left = np.percentile(slopes, 2.5)
right = np.percentile(slopes, 97.5)
fig, ax = subplots(figsize = (5, 5))
ax.hist(slopes, density=True, bins=50);
ax.hlines(y=0, xmin=left, xmax=right, color = 'yellow', linewidth = 7)
print('The 95% confidence interval on slope of true line is: (', left, ', ', right, ')')

The 95% confidence interval on slope of true line is: ( -0.09600622904085086 ,  0.27336580374007435 )

Since the 95% confidence interval includes 0, we cannot reject the null hypothesis. So, conclusion is that guessing birth weight based on the maternal age is not a good idea.

# Suppose we want a prediction for a particular value of the predictor variable
# We can bootstrap the prediction to find the 95% confidence interval on the prediction
# Here we use the prediction of Birth Weight using the number of Gestational days

num_repetitions = 5000
x = 300 # The value of the predictor variable (gestational days)

predictions = []
for i in range(num_repetitions):
    sample = baby.sample(frac = 1, replace=True)
    coefficients = np.polyfit(sample['Gestational Days'], sample['Birth Weight'], 1)
    y = coefficients[0]*x + coefficients[1]
    # instead of collecting the slopes as we did earlier, we collect the predictions
    predictions.append(y)

# plot a histogram of the slopes
#calculate the percentiles
predictions = np.array(predictions)
left = np.percentile(predictions, 2.5)
right = np.percentile(predictions, 97.5)
fig, ax = subplots(figsize = (5, 5))
ax.hist(predictions, density=True, bins=50);
ax.hlines(0, left, right, color = 'yellow', linewidth = 7)
print('The 95% confidence interval on the prediction for Birth Weight at 300 gestational days:', '(', left, ',', right, ')')

The 95% confidence interval on the prediction for Birth Weight at 300 gestational days: ( 127.34584318751573 , 131.36402218780412 )

Classification

# We build the nearest neighbor classifier, which is the simplest classification alorithm
# We will work with the data on Chronic Kidney Disease (CKD)
ckd = pd.read_csv('ckd.csv')
ckd.head(3)

	Age	Blood Pressure	Specific Gravity	Albumin	Sugar	Red Blood Cells	Pus Cell	Pus Cell clumps	Bacteria	Blood Glucose Random	...	Packed Cell Volume	White Blood Cell Count	Red Blood Cell Count	Hypertension	Diabetes Mellitus	Coronary Artery Disease	Appetite	Pedal Edema	Anemia	Class
0	48	70	1.005	4	0	normal	abnormal	present	notpresent	117	...	32	6700	3.9	yes	no	no	poor	yes	yes	1
1	53	90	1.020	2	0	abnormal	abnormal	present	notpresent	70	...	29	12100	3.7	yes	yes	no	poor	no	yes	1
2	63	70	1.010	3	0	abnormal	abnormal	present	notpresent	380	...	32	4500	3.8	yes	yes	no	poor	yes	no	1

3 rows × 25 columns

# Rename the columns so that it becomes easier to handle later
ckd = ckd.rename(columns = {'Blood Glucose Random': 'Glucose'})
ckd = ckd.rename(columns = {'White Blood Cell Count': 'WBCC'})
ckd_test = ckd[['Glucose', 'Hemoglobin', 'WBCC', 'Class']]

# Add columns for readings in standard units
ckd_test = ckd_test.assign(Glucose_in_std = (ckd_test['Glucose'] - ckd_test['Glucose'].mean())/ckd_test['Glucose'].std())
ckd_test = ckd_test.assign(WBCC_in_std = (ckd_test['WBCC'] - ckd_test['WBCC'].mean())/ckd_test['WBCC'].std())
ckd_test = ckd_test.assign(Hemoglobin_in_std = (ckd_test['Hemoglobin'] - ckd_test['Hemoglobin'].mean())/ckd_test['Hemoglobin'].std())

ckd_test.head(3)

	Glucose	Hemoglobin	WBCC	Class	Glucose_in_std	WBCC_in_std	Hemoglobin_in_std
0	117	11.2	6700	1	-0.220847	-0.567962	-0.863000
1	70	9.5	12100	1	-0.944594	1.158999	-1.452826
2	380	10.8	4500	1	3.829056	-1.271539	-1.001783

# Let us see if we can make a prediction using the Glucose and Hemoglobin readings only
# First, we do a scatter plot where we color every data point with the class (ckd/no ckd)
# to see if there is a pattern that points to the success of the nearest neighbor classifier
fig, ax = subplots(figsize = (5, 5))
scatter = ax.scatter(ckd_test['Hemoglobin_in_std'], ckd_test['Glucose_in_std'], c = ckd_test['Class'], alpha = 0.5)
ax.set_xlabel('Hemoglobin')
ax.set_ylabel('Glucose')
handles, labels = scatter.legend_elements()
# You can define custom labels here if "1", "2", etc. are not descriptive enough
custom_labels = ["No CKD", "CKD"] 
# Add the legend using the generated handles and custom labels
ax.legend(handles = handles, labels=custom_labels, loc="upper right", title="Classes");

For a new observation, the nearest neighbor classifier based on Glucose and Hemoglobin, assigns the same label as that of its nearest neighbor (where the similarity is based on distance in the above 2-D graph) in the training set. Given this, a decision boundary gets defined in the graph, which is a boundary across which the class/label changes when classifying based on nearest neighbor. Try drawing the decision boundary for the above graph.

Let's try to see what decision boundary we get when classifying based on Glucose and WBCC.

fig, ax = subplots(figsize = (5, 5))
scatter = ax.scatter(ckd_test['WBCC_in_std'], ckd_test['Glucose_in_std'], c = ckd_test['Class'], alpha = 0.5)
ax.set_xlabel('WBCC')
ax.set_ylabel('Glucose')
handles, labels = scatter.legend_elements()
# You can define custom labels here if "1", "2", etc. are not descriptive enough
custom_labels = ["No CKD", "CKD"] 
# Add the legend using the generated handles and custom labels
ax.legend(handles = handles, labels=custom_labels, loc="upper right", title="Classes");

We see that there is a region of overlaps, which causes the decision boundary to be not as neat as when classifying based on Glucose and Hemoglobin. In such cases, the KNN classifier, where we take the majority vote from the k nearest neighbors, helps. More discussion on this to follow in the next class. Stay tuned...