Descriptive Statistics

Bivariate Analysis

Introduction

This is our first lab when we are considering 2 dimensions and instead of calculating univariate statistics by groups (or factors) of other variable - we will measure their common relationships based on co-variance and correlation coefficients.

*Please be very careful when choosing the measure of correlation! In case of different measurument scales we have to recode one of the variables into weaker scale.

It would be nice to add some additional plots in the background. Feel free to add your own sections and use external packages.

Data

This time we are going to use a typical credit scoring data with predefined “default” variables and personal demografic and income data. Please take a look closer at headers and descriptions of each variable.

Scatterplots

First let’s visualize our quantitative relationships using scatterplots.

You can also normalize the skewed distribution of incomes using log:

We can add an estimated linear regression line:

## `geom_smooth()` using formula = 'y ~ x'

Scatterplots by groups

We can finally see if there any differences between risk status:

## `geom_smooth()` using formula = 'y ~ x'

We can also see more closely if there any differences between those two distributions adding their estimated density plots:

We can also put those plots together:

## Warning: The dot-dot notation (`..count..`) was deprecated in ggplot2 3.4.0.
## ℹ Please use `after_stat(count)` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

Scatterplots with density curves

We can also see more closely if there any differences between those two distributions adding their estimated density plots:

Correlation coefficients - Pearson’s linear correlation

Ok, let’s move to some calculations. In R, we can use the cor() function. It takes three arguments and the method: cor(x, y, method) For 2 quantitative data, with all assumptions met, we can calculate simple Pearson’s coefficient of linear correlation:

## [1] 0.574346

Ok, what about the percentage of the explained variability?

##    estimate      p.value statistic   n gp  Method
## 1 0.3194263 4.805085e-18  8.899323 700  1 pearson

So as we can see almost ??? of total log of incomes’ variability is explained by differences in age. The rest (???) is probably explained by other factors.

Partial and semipartial correlation

The partial and semi-partial (also known as part) correlations are used to express the specific portion of variance explained by eliminating the effect of other variables when assessing the correlation between two variables.

Partial correlation holds constant one variable when computing the relations to others. Suppose we want to know the correlation between X and Y holding Z constant for both X and Y. That would be the partial correlation between X and Y controlling for Z.

Semipartial correlation holds Z constant for either X or Y, but not both, so if we wanted to control X for Z, we could compute the semipartial correlation between X and Y holding Z constant for X.

Suppose we want to know the correlation between the log of income and age controlling for years of employment. How highly correlated are these after controlling for tenure?

**There can be more than one control variable.

##    estimate      p.value statistic   n gp  Method
## 1 0.3194263 4.805085e-18  8.899323 700  1 pearson

## [1] 0.1859717

How can we interpret the obtained partial correlation coefficient? What is the difference between that one and the semi-partial coefficient:

Rank correlation

For 2 different scales - like for example this pair of variables: income vs. education levels - we cannot use Pearson’s coefficient. The only possibility is to rank also incomes… and lose some more detailed information about them.

First, let’s see boxplots of income by education levels.

Now, let’s see Kendal’s coefficient of rank correlation (robust for ties).

## [1] 0.1577567

Point-biserial correlation

Let’s try to verify if there is a significant relationship between incomes and risk status. First, let’s take a look at the boxplot:

If you would like to compare 1 quantitative variable (income) and 1 dychotomous variable (default status - binary), then you can use point-biserial coefficient:

## [1] -0.07096966

Nonlinear correlation - eta coefficient

If you would like to check if there are any nonlinearities between 2 variables, the only possibility (beside transformations and linear analysis) is to calculate “eta” coefficient and compare it with the Pearson’s linear coefficient.

Correlation matrix

We can also prepare the correlation matrix for all quantitative variables stored in our data frame.

We can use ggcorr() function:

## corrplot 0.92 loaded

As you can see - the default correlation matrix is not the best idea for all measurement scales (including binary variable “default”).

That’s why now we can perform our bivariate analysis with ggpair with grouping.

Correlation matrix with scatterplots

Here is what we are about to calculate: - The correlation matrix between age, log_income, employ, address, debtinc, creddebt, and othdebt variable grouped by whether the person has a default status or not. - Plot the distribution of each variable by group - Display the scatter plot with the trend by group

## `geom_smooth()` using formula = 'y ~ x'

Qualitative data

In case of two variables measured on nominal or ordinal&nominal scale - we are forced to organize so called “contingency” table with frequencies and calculate some kind of the correlation coefficient based on them. This is so called “contingency analysis”.

Let’s consider one example based on our data: verify, if there is any significant correlation between education level and credit risk.

default <- c(bank$def)
educaion <- c(bank$ed)

contingency_table <- table(educaion, default)
print(contingency_table)

##         default
## educaion   0   1
##        1 293  79
##        2 139  59
##        3  57  30
##        4  24  14
##        5   4   1

chi_square <- chisq.test(contingency_table)

## Warning in chisq.test(contingency_table): Chi-squared approximation may be
## incorrect

print(chi_square)

## 
##  Pearson's Chi-squared test
## 
## data:  contingency_table
## X-squared = 11.492, df = 4, p-value = 0.02155

contingency_coefficient <- sqrt(chi_square$statistic / (sum(contingency_table) * min(nrow(contingency_table), ncol(contingency_table)) ))
print(contingency_coefficient)

##  X-squared 
## 0.09060251

Exercise 1. Contingency analysis.

Do you believe in the Afterlife? https://nationalpost.com/news/canada/millennials-do-you-believe-in-life-after-life A survey was conducted and a random sample of 1091 questionnaires is given in the form of the following contingency table:

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

Our task is to check if there is a significant relationship between the belief in the afterlife and gender. We can perform this procedure with the simple chi-square statistics and chosen qualitative correlation coefficient (two-way 2x2 table).

## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  dane
## X-squared = 0.11103, df = 1, p-value = 0.739

##         Believe
## Gender         Yes        No
##   Female 0.3987168 0.3437214
##   Male   0.1347388 0.1228231

As you can see we can calculate our chi-square statistic really quickly for two-way tables or larger. Now we can standardize this contingency measure to see if the relationship is significant.

## [1] 0.01218871

Exercise 2. Contingency analysis for the ‘Titanic’ data.

Let’s consider the titanic dataset which contains a complete list of passengers and crew members on the RMS Titanic. It includes a variable indicating whether a person did survive the sinking of the RMS Titanic on April 15, 1912. A data frame contains 2456 observations on 14 variables.

The website http://www.encyclopedia-titanica.org/ offers detailed information about passengers and crew members on the RMS Titanic. According to the website 1317 passengers and 890 crew member were aboard.

8 musicians and 9 employees of the shipyard company are listed as passengers, but travelled with a free ticket, which is why they have NA values in fare. In addition to that, fare is truely missing for a few regular passengers.

library(naniar)

titanic$Age[is.na(titanic$Age)] <- mean(titanic$Age, na.rm = TRUE)
titanic$Year.of.Birth[is.na(titanic$Year.of.Birth)] <- mean(titanic$Year.of.Birth, na.rm = TRUE)
titanic$Fare_GBP[is.na(titanic$Fare_GBP)] <- mean(titanic$Fare_GBP, na.rm = TRUE)
titanic$Fare_today[is.na(titanic$Fare_today)] <- mean(titanic$Fare_today, na.rm = TRUE)

titanic$Age <- floor(titanic$Age)
titanic$Year.of.Birth <- floor(titanic$Year.of.Birth)

vis_miss(titanic)

## Titanic Heatmap

Below is the heatmap for the correlations between numerical variables in the titanic dataset.

correlation_matrix_titanic <- cor(titanic[,sapply(titanic, is.numeric)])

corrplot(correlation_matrix_titanic, method = 'number')

Status and Age

titanic$Status <- replace(titanic$Status, titanic$Status == "", NA)
mode_status <- names(sort(table(titanic$Status), decreasing = TRUE))[1]
titanic$Status[is.na(titanic$Status)] <- mode_status

The most appropriate plot to visualize Status and Age variables is a box plot, because Status has only two unique observations. The median age of Survivors is slightly lower than the median age of Victims. It might indicate that younger individuals are more likely to survive.

ggplot(titanic, aes(x = Status, y = Age, group = Status)) +
  geom_boxplot()

To calculate the correlation between those two variables, the contingency table is created and then the contingency coefficient is calculated. A value of 0.21 indicates a weak corelation between Status and Age. The distribution of Status varaible is not strongly dependent on the Age variable.

contingency_table_titanic <- table(titanic$Age, titanic$Status)

chi_square_titanic <- chisq.test(contingency_table_titanic)
print(chi_square_titanic)

## 
##  Pearson's Chi-squared test
## 
## data:  contingency_table_titanic
## X-squared = 116.97, df = 71, p-value = 0.0004876

total_obs <- sum(contingency_table_titanic)

contingency_coefficient_titanic <- sqrt(chi_square_titanic$statistic / (chi_square_titanic$statistic + total_obs))
print(contingency_coefficient_titanic)

## X-squared 
##  0.213217

Status and Gender

The most appropriate plot to visualize Status and Gender variables is a grouped bar plot. The distribution of the Survivor status within Male gender is slightly smaller than within Female gender. It might indicate that Females were more likely to survive. There are significantly less Male survivors than Male victims, however the number of Survivors for Female is significantly greater than the number of Victims for Female. It suggests that gender could influence the likelihood of survival.

ggplot(titanic, aes(x = Gender, fill = Status)) +
  geom_bar()

For the Status and Gender the contingency table is also created and the contingency coefficient is calculated. A contingency coefficient of 0.42 suggests a moderate or strong correlation between Status and Gender. The variables are associated with each other. Knowing the gender of a person provides can be useful in predicting their Status.

contingency_table_titanic_2 <- table(titanic$Status, titanic$Gender)
print(contingency_table_titanic_2)

##           
##            Female Male
##   Survivor    359  352
##   Victim      145 1600

total_obs_2 <- sum(contingency_table_titanic_2)

chi_square_titanic_2 <- chisq.test(contingency_table_titanic_2)
print(chi_square_titanic_2)

## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  contingency_table_titanic_2
## X-squared = 548.55, df = 1, p-value < 2.2e-16

contingency_coefficient_titanic_2 <- sqrt(chi_square_titanic_2$statistic / (chi_square_titanic_2$statistic  + total_obs_2))
print(contingency_coefficient_titanic_2)

## X-squared 
## 0.4272847

Status and Crew.or.Passenger

ggplot(titanic, aes(x = titanic$Crew.or.Passenger, fill = Status)) +
  geom_bar() +
  labs(x = "Crew or Passenger")

A value 0.19 of a contingency coefficient between Status and Crew.or.Passenger variables indicates low correlation between these two variables.

contingency_table_titanic_3 <- table(titanic$Status, titanic$Crew.or.Passenger)
print(contingency_table_titanic_3)

##           
##            Crew Passenger
##   Survivor  211       500
##   Victim    893       852

total_obs_3 <- sum(contingency_table_titanic_3)

chi_square_titanic_3 <- chisq.test(contingency_table_titanic_3)
print(chi_square_titanic_3)

## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  contingency_table_titanic_3
## X-squared = 93.486, df = 1, p-value < 2.2e-16

contingency_coefficient_titanic_3 <- sqrt(chi_square_titanic_3$statistic / (chi_square_titanic_3$statistic  + total_obs_3))
print(contingency_coefficient_titanic_3)

## X-squared 
## 0.1914904

Status and Class…Department

ggplot(titanic, aes(x = titanic$Class...Department, fill = Status)) +
  geom_bar() +
  labs(x = "Class/Department")

A contingency coefficient of 0.30 indicates that there is a moderate/meaningful level of correlation between Status and Class…Department variables. The distribution of survival status is dependent on the class or department.

contingency_table_titanic_4 <- table(titanic$Status, titanic$Class...Department)
print(contingency_table_titanic_4)

##           
##            1st Class 2nd Class 3rd Class Deck Crew Engineering Crew
##   Survivor       201       119       180        43               71
##   Victim         149       175       528        56              432
##           
##            Restaurant Staff Victualling Crew
##   Survivor                3               94
##   Victim                 66              339

total_obs_4 <- sum(contingency_table_titanic_4)

chi_square_titanic_4 <- chisq.test(contingency_table_titanic_4)
print(chi_square_titanic_4)

## 
##  Pearson's Chi-squared test
## 
## data:  contingency_table_titanic_4
## X-squared = 256.53, df = 6, p-value < 2.2e-16

contingency_coefficient_titanic_4 <- sqrt(chi_square_titanic_4$statistic / (chi_square_titanic_4$statistic  + total_obs_4))
print(contingency_coefficient_titanic_4)

## X-squared 
## 0.3075266