Descriptive Statistics

Bivariate Analysis - Report_4: Python

Libraries

import requests
import pyreadstat
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
import seaborn as sns
from scipy.stats import pearsonr
import pingouin as pg
from scipy.stats import kendalltau
from scipy.stats import pointbiserialr
from scipy.stats import chi2_contingency
import statsmodels.api as sm
from statsmodels.graphics.mosaicplot import mosaic

Downloading the file

#import requests

# URL and destination file name for downloading the data
url = "https://github.com/kflisikowski/ds/blob/master/bank_defaults.sav?raw=true"
destfile = "bank_defaults.sav"

# Download the file
response = requests.get(url)
if response.status_code == 200:
    with open(destfile, 'wb') as file:
        file.write(response.content)
    print("File downloaded successfully.")
else:
    print("Failed to download the file. Status code:", response.status_code)

## 53033
## File downloaded successfully.

Reading the .sav file on Python

# import pyreadstat

bank_defaults, meta = pyreadstat.read_sav(destfile)

Dropping missing values

#import pandas as pd

bank = bank_defaults.dropna()

Converting default column to a categorical type.

#Also we’re checking the data with .head() & .dtypes functions / instances from Python.

# Convert 'default' column to a categorical type
bank['def'] = bank['default'].astype('category')

## <string>:3: SettingWithCopyWarning: 
## A value is trying to be set on a copy of a slice from a DataFrame.
## Try using .loc[row_indexer,col_indexer] = value instead
## 
## See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy

# Print the first few rows and data types of the DataFrame
print(bank.head())

##     age   ed  employ  address  ...  preddef1  preddef2  preddef3  def
## 0  41.0  3.0    17.0     12.0  ...  0.808394  0.788640  0.213043  1.0
## 1  27.0  1.0    10.0      6.0  ...  0.198297  0.128445  0.436903  0.0
## 2  40.0  1.0    15.0     14.0  ...  0.010036  0.002987  0.141023  0.0
## 3  41.0  1.0    15.0     14.0  ...  0.022138  0.010273  0.104422  0.0
## 4  24.0  2.0     2.0      0.0  ...  0.781588  0.737885  0.436903  1.0
## 
## [5 rows x 13 columns]

print(bank.dtypes)

## age          float64
## ed           float64
## employ       float64
## address      float64
## income       float64
## debtinc      float64
## creddebt     float64
## othdebt      float64
## default      float64
## preddef1     float64
## preddef2     float64
## preddef3     float64
## def         category
## dtype: object

#We can see the plot:

# Scatter plot of original data
plt.figure(figsize=(7, 6))
plt.scatter(bank['income'], bank['debtinc'], label='Data Points')


# Add title and axis labels
plt.title("Income vs. Debt Income")
plt.xlabel("Income")
plt.ylabel("Debt Income")

# Show the plot
plt.show()

Applying logarithmic transformation:

#import numpy as np

#Apply logarithmic transformation
bank.loc[:, 'log_income'] = np.log(bank.loc[:, 'income'])

## <string>:4: SettingWithCopyWarning: 
## A value is trying to be set on a copy of a slice from a DataFrame.
## Try using .loc[row_indexer,col_indexer] = value instead
## 
## See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy

bank.loc[:, 'log_debtinc'] = np.log(bank.loc[:, 'debtinc'])

## <string>:1: SettingWithCopyWarning: 
## A value is trying to be set on a copy of a slice from a DataFrame.
## Try using .loc[row_indexer,col_indexer] = value instead
## 
## See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy

Normalizing the skewed distribution of incomes using log:

# Scatter plot of original data
plt.figure(figsize=(7, 6))
plt.scatter(bank['log_income'], bank['log_debtinc'], label='Data Points')


# Add title and axis labels
plt.title("Log Income vs. Log Debt Income")
plt.xlabel("Log Income")
plt.ylabel("Log Debt Income")

# Show the plot
plt.show()

Creating a scatter plot with the logarithmically transformed data

# Scatter plot of logarithmically transformed data
plt.figure(figsize=(7, 6))
plt.scatter(bank['log_income'], bank['log_debtinc'], label='Data Points')

# Perform linear regression on logarithmically transformed data
slope_log, intercept_log, r_value_log, p_value_log, std_err_log = stats.linregress(bank['log_income'], bank['log_debtinc'])

# Generate values for the regression line
x_log = np.linspace(bank['log_income'].min(), bank['log_income'].max(), 100)
y_log = slope_log * x_log + intercept_log

# Plot the regression line
plt.plot(x_log, y_log, color='red', label='Regression Line')

# Add title and axis labels
plt.title("Logarithmic Income vs. Logarithmic Debt Income with Regression Line")
plt.xlabel("Logarithmic Income")
plt.ylabel("Logarithmic Debt Income")

# Add a legend
plt.legend()

# Show the plot
plt.show()

Scatterplots by groups

# Grouping the data by education level
grouped_data = bank.groupby('ed')

# Creating separate scatterplots for each education level
for name, group in grouped_data:
    plt.figure(figsize=(7, 6))
    plt.scatter(group['log_income'], group['log_debtinc'], label='Data Points')
    
    # Perform linear regression
    slope, intercept, r_value, p_value, std_err = stats.linregress(group['log_income'], group['log_debtinc'])
    
    # Generate values for the regression line
    x = np.linspace(group['log_income'].min(), group['log_income'].max(), 100)
    y = slope * x + intercept
    
    # Plot the regression line
    plt.plot(x, y, color='red', label='Regression Line')
    
    # Add title and axis labels
    plt.title(f"Logarithmic Income vs. Logarithmic Debt Income (Education Level: {name})")
    plt.xlabel("Logarithmic Income")
    plt.ylabel("Logarithmic Debt Income")
    
    # Add a legend
    plt.legend()
    
    # Show the plot
    plt.show()

Estimated Density Plotes

#To see more closely if there any differences between those two distributions adding their estimated density plots

# Scatter plot with color by education level
plt.figure(figsize=(10, 6))
sns.scatterplot(x=bank['income'], y=bank['debtinc'], hue=bank['ed'], palette='viridis')

# Add title and axis labels
plt.title("Income vs. Debt Income")
plt.xlabel("Income")
plt.ylabel("Debt Income")

# Show the plot
plt.show()

# Marginal density plot of age
plt.figure(figsize=(8, 6))
sns.histplot(bank['age'], kde=True)

# Add title and axis labels
plt.title("Density Plot of Age")
plt.xlabel("Age")
plt.ylabel("Density")

# Show the plot
plt.show()

# Marginal density plot of y
plt.figure(figsize=(8, 6))
sns.histplot(bank_defaults['age'], kde=True)

# Add title and axis labels
plt.title("Density Plot of Y")
plt.xlabel("Y")
plt.ylabel("Density")

# Show the plot
plt.show()

Plots Together:

import matplotlib.pyplot as plt
import seaborn as sns

# Set up the layout
plt.figure(figsize=(12, 10))  # Adjust the figure size as needed

# Scatter plot
plt.subplot(2, 2, 1)
sns.scatterplot(x=bank['income'], y=bank['debtinc'], hue=bank['ed'], palette='viridis')
plt.title("Income vs. Debt Income")
plt.xlabel("Income")
plt.ylabel("Debt Income")



# Marginal density plot of age (top panel)
plt.subplot(2, 2, 2)
sns.histplot(bank['age'], kde=True)
plt.title("Density Plot of Age")
plt.xlabel("Age")
plt.ylabel("Density")

# Create space for the third plot
plt.subplot(2, 2, 3)
plt.axis('off')  # Turn off the axis for this subplot

## (0.0, 1.0, 0.0, 1.0)

# Marginal density plot of y (right panel)
plt.subplot(2, 2, 4)
sns.histplot(bank_defaults['age'], kde=True)
plt.title("Density Plot of Y")
plt.xlabel("Y")
plt.ylabel("Density")

# Adjust layout
plt.tight_layout()

# Show the plot
plt.show()

##Giving Density Curves to Scatterplots

# Scatter plot
plt.figure(figsize=(10, 6))
plt.scatter(bank['income'], bank['debtinc'], c=bank['ed'], cmap='viridis')
plt.title("Income vs. Debt Income")
plt.xlabel("Income")
plt.ylabel("Debt Income")

# Density plot for income
sns.kdeplot(bank['income'], color='blue', linestyle='-')

# Density plot for debtinc
sns.kdeplot(bank['debtinc'], color='red', linestyle='-')

# Show the plot
plt.show()

Correlation coefficients - P.L correlation:

#from scipy.stats import pearsonr

# Calculate Pearson's correlation coefficient
correlation, _ = pearsonr(bank['income'], bank['debtinc'])

# Print the correlation coefficient
print("Pearson's correlation coefficient:", correlation)

## Pearson's correlation coefficient: -0.026777293424863746

Percentage of the explained variability:

#from scipy.stats import pearsonr

# Calculate Pearson's correlation coefficient
correlation, _ = pearsonr(bank['income'], bank['debtinc'])

# Calculate the percentage of explained variability
explained_variability = correlation**2 * 100

# Print the percentage of explained variability
print("Percentage of explained variability:", explained_variability)

## Percentage of explained variability: 0.07170234431612511

Partial and semipartial correlation:

#import pingouin as pg

# Ensure there are no missing values in the variables
complete_cases = bank.iloc[:, [bank.columns.get_loc('income'), bank.columns.get_loc('age'), bank.columns.get_loc('employ')]].dropna()

# Extract variables
income = complete_cases.iloc[:, complete_cases.columns.get_loc('income')]
age = complete_cases.iloc[:, complete_cases.columns.get_loc('age')]
tenure = complete_cases.iloc[:, complete_cases.columns.get_loc('employ')]

# Compute partial correlation
partial_corr = pg.partial_corr(data=complete_cases, x='income', y='age', covar='employ', method='pearson')

# Print the partial correlation coefficient
print("Partial correlation coefficient:", partial_corr['r'][0])

## <string>:3: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`
## Partial correlation coefficient: 0.220812408872544

Difference between that one and the S-P coefficient:

# Ensure there are no missing values in the variables
complete_cases = bank[['income', 'age', 'employ']].dropna()

# Compute semi-partial (part) correlation
semi_partial_corr = pg.partial_corr(data=complete_cases, x='income', y='age', covar='employ', method='pearson')

# Print the semi-partial correlation coefficient
print("Semi-partial (part) correlation coefficient:", semi_partial_corr['r'].iloc[0])

## Semi-partial (part) correlation coefficient: 0.220812408872544

Rank correlation

# Create the boxplot
plt.figure(figsize=(10, 6))

## <Figure size 1000x600 with 0 Axes>

sns.boxplot(x='ed', y='income', data=bank, hue='ed', palette='light:#5A9', dodge=False)

## <Axes: xlabel='ed', ylabel='income'>

plt.legend([],[], frameon=False)

## <matplotlib.legend.Legend object at 0x0000020654B26D20>

# Customize the plot
plt.title('Income by Education Levels')

## Text(0.5, 1.0, 'Income by Education Levels')

plt.xlabel('Education Level')

## Text(0.5, 0, 'Education Level')

plt.ylabel('Income')

## Text(0, 0.5, 'Income')

# Show the plot
plt.show()

Kendal’s coefficient of rank correlation:

#(robust for ties)

#from scipy.stats import kendalltau

kendall_corr, _ = kendalltau(bank['income'], bank['debtinc'])

# Print the Kendall's coefficient of rank correlation
print("Kendall's coefficient of rank correlation:", kendall_corr)

## Kendall's coefficient of rank correlation: -0.012242089172559398

Point-biserial correlation

#I'm not sure in these codes. I hope it's true.

# Define risk_status based on debtinc thresholds
bank.loc[bank['debtinc'] <= 15, 'risk_status'] = 'Low'

## <string>:5: SettingWithCopyWarning: 
## A value is trying to be set on a copy of a slice from a DataFrame.
## Try using .loc[row_indexer,col_indexer] = value instead
## 
## See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy

bank.loc[bank['debtinc'] > 15, 'risk_status'] = 'High'

# Create the boxplot
plt.figure(figsize=(10, 6))

## <Figure size 1000x600 with 0 Axes>

sns.boxplot(x='risk_status', y='income', data=bank, hue='risk_status', palette='light:#5A9', dodge=False)

## <Axes: xlabel='risk_status', ylabel='income'>

plt.legend([],[], frameon=False)

## <matplotlib.legend.Legend object at 0x00000206544BA0C0>

# Customize the plot
plt.title('Income by Risk Status')

## Text(0.5, 1.0, 'Income by Risk Status')

plt.xlabel('Risk Status')

## Text(0.5, 0, 'Risk Status')

plt.ylabel('Income')

## Text(0, 0.5, 'Income')

# Show the plot
plt.show()

Comparing QVar and 1 DVar

#from scipy.stats import pointbiserialr

point_biserial_corr, _ = pointbiserialr(bank['income'], bank['default'])

# Print the point-biserial correlation coefficient
print("Point-biserial correlation coefficient:", point_biserial_corr)

## Point-biserial correlation coefficient: -0.07096965661922319

Eta Coefficient:

# Calculating Pearson's linear correlation coefficient
pearson_corr, _ = pearsonr(bank['creddebt'], bank['othdebt'])

# Calculate Kendall's rank correlation coefficient
kendall_corr, _ = kendalltau(bank['creddebt'], bank['othdebt'])

# Calculate the eta coefficient
eta_coefficient = np.sqrt(1 - (kendall_corr**2 / pearson_corr**2))

# Print the eta coefficient
print("Eta coefficient:", eta_coefficient)

## Eta coefficient: 0.7083779812539175

Correlation matrix

# Ensuring all variables are numeric
numeric_columns = bank.select_dtypes(include=['number'])

# Calculate the correlation matrix for quantitative variables
correlation_matrix = numeric_columns.corr()

# Print the correlation matrix
print("Correlation Matrix:")

## Correlation Matrix:

print(correlation_matrix)

##                   age        ed    employ  ...  preddef3  log_income  log_debtinc
## age          1.000000  0.022325  0.536497  ...  0.013103    0.574346     0.020615
## ed           0.022325  1.000000 -0.153621  ...  0.007231    0.223493     0.009844
## employ       0.536497 -0.153621  1.000000  ... -0.025159    0.723908    -0.044799
## address      0.597591  0.056919  0.322334  ...  0.007748    0.368778     0.017506
## income       0.478710  0.235190  0.619681  ... -0.026165    0.893931    -0.023819
## debtinc      0.016398  0.008838 -0.031182  ...  0.989753   -0.014848     0.903224
## creddebt     0.295207  0.088274  0.403694  ...  0.496677    0.529251     0.442019
## othdebt      0.340217  0.165459  0.406091  ...  0.575755    0.601675     0.535911
## default     -0.137657  0.114676 -0.282978  ...  0.390154   -0.135226     0.340456
## preddef1    -0.314406  0.165519 -0.489568  ...  0.676668   -0.240918     0.587082
## preddef2    -0.228460  0.190319 -0.469639  ...  0.650348   -0.220286     0.559461
## preddef3     0.013103  0.007231 -0.025159  ...  1.000000   -0.012468     0.843109
## log_income   0.574346  0.223493  0.723908  ... -0.012468    1.000000    -0.020960
## log_debtinc  0.020615  0.009844 -0.044799  ...  0.843109   -0.020960     1.000000
## 
## [14 rows x 14 columns]

# Plotting the correlation matrix
plt.figure(figsize=(10, 8))

## <Figure size 1000x800 with 0 Axes>

sns.heatmap(correlation_matrix, annot=True, cmap='coolwarm', fmt=".2f")

## <Axes: >

plt.title('Correlation Matrix')

## Text(0.5, 1.0, 'Correlation Matrix')

plt.show()

Correlation matrix with scatterplots

# Step 1: Calculating the correlation matrix by group
correlation_matrix = bank_defaults.groupby('default').corr()

# Step 2: Plot the distribution of each variable by group
fig, axes = plt.subplots(nrows=len(bank_defaults.columns) - 1, ncols=1, figsize=(8, 6))
for i, col in enumerate(bank_defaults.columns[:-1]):
    sns.histplot(data=bank_defaults, x=col, hue='default', ax=axes[i], kde=True)
    axes[i].set_title(f'Distribution of {col} by default status')

## <Axes: xlabel='age', ylabel='Count'>
## Text(0.5, 1.0, 'Distribution of age by default status')
## <Axes: xlabel='ed', ylabel='Count'>
## Text(0.5, 1.0, 'Distribution of ed by default status')
## <Axes: xlabel='employ', ylabel='Count'>
## Text(0.5, 1.0, 'Distribution of employ by default status')
## <Axes: xlabel='address', ylabel='Count'>
## Text(0.5, 1.0, 'Distribution of address by default status')
## <Axes: xlabel='income', ylabel='Count'>
## Text(0.5, 1.0, 'Distribution of income by default status')
## <Axes: xlabel='debtinc', ylabel='Count'>
## Text(0.5, 1.0, 'Distribution of debtinc by default status')
## <Axes: xlabel='creddebt', ylabel='Count'>
## Text(0.5, 1.0, 'Distribution of creddebt by default status')
## <Axes: xlabel='othdebt', ylabel='Count'>
## Text(0.5, 1.0, 'Distribution of othdebt by default status')
## <Axes: xlabel='default', ylabel='Count'>
## Text(0.5, 1.0, 'Distribution of default by default status')
## <Axes: xlabel='preddef1', ylabel='Count'>
## Text(0.5, 1.0, 'Distribution of preddef1 by default status')
## <Axes: xlabel='preddef2', ylabel='Count'>
## Text(0.5, 1.0, 'Distribution of preddef2 by default status')

plt.tight_layout()

## <string>:2: UserWarning: Tight layout not applied. tight_layout cannot make Axes height small enough to accommodate all Axes decorations.

plt.show()

# Step 3: Displaying scatter plots with trend by group
sns.set(style="ticks")
sns.pairplot(data=bank_defaults, hue='default', diag_kind='kde')

## <seaborn.axisgrid.PairGrid object at 0x0000020651BE7E90>

plt.show()

Qualitative data

#from scipy.stats import chi2_contingency


# Creating a contingency table
contingency_table = pd.crosstab(bank_defaults['ed'], bank_defaults['default'])

# Perform chi-square test of independence
chi2, p, dof, expected = chi2_contingency(contingency_table)

# Print the results
print("Chi-square statistic:", chi2)

## Chi-square statistic: 11.492341703403657

print("P-value:", p)

## P-value: 0.02155395852574417

print("Degrees of freedom:", dof)

## Degrees of freedom: 4

print("Expected frequencies table:")

## Expected frequencies table:

print(expected)

## [[274.74857143  97.25142857]
##  [146.23714286  51.76285714]
##  [ 64.25571429  22.74428571]
##  [ 28.06571429   9.93428571]
##  [  3.69285714   1.30714286]]

Exercise 1

# Creating the contingency table
data = np.array([[435, 375], [147, 134]])
dane = pd.DataFrame(data, index=['Female', 'Male'], columns=['Yes', 'No'])

# Displaying the contingency table
print("Contingency Table:")

## Contingency Table:

print(dane)

##         Yes   No
## Female  435  375
## Male    147  134

# Chi-square test of independence
chi2_stat, p_val, dof, ex = chi2_contingency(dane)

print("\nChi-Square Test:")

## 
## Chi-Square Test:

print(f"Chi2 Stat: {chi2_stat}")

## Chi2 Stat: 0.1110272160868229

print(f"P-value: {p_val}")

## P-value: 0.7389776820172238

print(f"Degrees of Freedom: {dof}")

## Degrees of Freedom: 1

print("Expected Frequencies:")

## Expected Frequencies:

print(ex)

## [[432.09899175 377.90100825]
##  [149.90100825 131.09899175]]

# Phi Coefficient
n = np.sum(data)
phi = np.sqrt(chi2_stat / n)
print(f"\nPhi Coefficient: {phi}")

## 
## Phi Coefficient: 0.010087936733575699

fig, ax = plt.subplots(figsize=(8,6))
mosaic(dane.stack(), title='Mosaic Plot of Belief in Afterlife by Gender', ax=ax)

## (<Figure size 800x600 with 3 Axes>, {('Female', 'Yes'): (0.0, 0.0, 0.7387444081152442, 0.5352528608342562), ('Female', 'No'): (0.0, 0.5385751199704688, 0.7387444081152442, 0.4614248800295311), ('Male', 'Yes'): (0.7437195324933537, 0.0, 0.2562804675066464, 0.5213936936191342), ('Male', 'No'): (0.7437195324933537, 0.5247159527553469, 0.2562804675066464, 0.47528404724465306)})

plt.show()

# Bar plot
dane.plot(kind='bar', stacked=True)

## <Axes: >

plt.title('Belief in Afterlife by Gender')

## Text(0.5, 1.0, 'Belief in Afterlife by Gender')

plt.xlabel('Gender')

## Text(0.5, 0, 'Gender')

plt.ylabel('Count')

## Text(0, 0.5, 'Count')

plt.show()

# Proportion table
prop_table = dane / dane.sum().sum()
print("\nProportion Table:")

## 
## Proportion Table:

print(prop_table)

##              Yes        No
## Female  0.398717  0.343721
## Male    0.134739  0.122823

##Exercise 2

# Downloading the CSV file
url = "https://github.com/kflisikowski/ds/blob/master/titanic.csv?raw=true"
titanic = pd.read_csv(url, sep=";", index_col=0)

#Dropping rows:

titanic.dropna(inplace=True)

# Convert columns to factors
titanic['Status'] = pd.Categorical(titanic['Status'])
titanic['Gender'] = pd.Categorical(titanic['Gender'])

#Creating Contingency Table

contingency_table = pd.crosstab(titanic['Gender'], titanic['Status'])

#Performing Chi-Square Test

chi2, p, dof, expected = chi2_contingency(contingency_table)

print("Chi-Square Test:")

## Chi-Square Test:

print("Chi2 Statistic:", chi2)

## Chi2 Statistic: 60.862575878954544

print("P-value:", p)

## P-value: 6.120129207752533e-15

print("Degrees of Freedom:", dof)

## Degrees of Freedom: 1

#Calculating Phi coefficient

n = contingency_table.sum().sum()
phi = np.sqrt(chi2 / n)

print("\nPhi Coefficient:", phi)

## 
## Phi Coefficient: 0.3264803973833301

#Calculating Cramer’s Value:

min_dim = min(contingency_table.shape) - 1
v = np.sqrt(chi2 / (n * min_dim))

print("Cramer's V:", v)

## Cramer's V: 0.3264803973833301

#Ploting mosaic plot:

plt.figure(figsize=(8, 6))

## <Figure size 800x600 with 0 Axes>

mosaic_plot = contingency_table.plot(kind='bar', stacked=True)
plt.title('Mosaic Plot')

## Text(0.5, 1.0, 'Mosaic Plot')

plt.xlabel('Gender')

## Text(0.5, 0, 'Gender')

plt.ylabel('Count')

## Text(0, 0.5, 'Count')

plt.show()

#Ploting bar plot:

plt.figure(figsize=(8, 6))

## <Figure size 800x600 with 0 Axes>

bar_plot = contingency_table.plot(kind='bar')
plt.title('Bar Plot')

## Text(0.5, 1.0, 'Bar Plot')

plt.xlabel('Gender')

## Text(0.5, 0, 'Gender')

plt.ylabel('Count')

## Text(0, 0.5, 'Count')

plt.show()