INFO-H510 Week 7 Data Dive

library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.4
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.4.4     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.0
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(ggplot2)
library(ggthemes)
library(ggrepel)
library(effsize)
library(pwrss)
## 
## Attaching package: 'pwrss'
## 
## The following object is masked from 'package:stats':
## 
##     power.t.test
#Loading data
data <- read_delim("../Churn_Modelling.csv", delim = ",")
## Rows: 10002 Columns: 14
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr  (3): Surname, Geography, Gender
## dbl (11): RowNumber, CustomerId, CreditScore, Age, Tenure, Balance, NumOfPro...
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
str(data)
## spc_tbl_ [10,002 × 14] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
##  $ RowNumber      : num [1:10002] 1 2 3 4 5 6 7 8 9 10 ...
##  $ CustomerId     : num [1:10002] 15634602 15647311 15619304 15701354 15737888 ...
##  $ Surname        : chr [1:10002] "Hargrave" "Hill" "Onio" "Boni" ...
##  $ CreditScore    : num [1:10002] 619 608 502 699 850 645 822 376 501 684 ...
##  $ Geography      : chr [1:10002] "France" "Spain" "France" "France" ...
##  $ Gender         : chr [1:10002] "Female" "Female" "Female" "Female" ...
##  $ Age            : num [1:10002] 42 41 42 39 43 44 50 29 44 NA ...
##  $ Tenure         : num [1:10002] 2 1 8 1 2 8 7 4 4 2 ...
##  $ Balance        : num [1:10002] 0 83808 159661 0 125511 ...
##  $ NumOfProducts  : num [1:10002] 1 1 3 2 1 2 2 4 2 1 ...
##  $ HasCrCard      : num [1:10002] 1 0 1 0 NA 1 1 1 0 1 ...
##  $ IsActiveMember : num [1:10002] 1 1 0 0 1 0 1 0 NA 1 ...
##  $ EstimatedSalary: num [1:10002] 101349 112543 113932 93827 79084 ...
##  $ Exited         : num [1:10002] 1 0 1 0 0 1 0 1 0 0 ...
##  - attr(*, "spec")=
##   .. cols(
##   ..   RowNumber = col_double(),
##   ..   CustomerId = col_double(),
##   ..   Surname = col_character(),
##   ..   CreditScore = col_double(),
##   ..   Geography = col_character(),
##   ..   Gender = col_character(),
##   ..   Age = col_double(),
##   ..   Tenure = col_double(),
##   ..   Balance = col_double(),
##   ..   NumOfProducts = col_double(),
##   ..   HasCrCard = col_double(),
##   ..   IsActiveMember = col_double(),
##   ..   EstimatedSalary = col_double(),
##   ..   Exited = col_double()
##   .. )
##  - attr(*, "problems")=<externalptr>

Null Hypothesis 1:

Hypothesis: There is no significant difference in the average credit score between customers who have churned and those who haven’t.

data |>
  ggplot() +
  geom_boxplot(mapping = 
                 aes(x = CreditScore, 
                     y = factor(Exited, levels = c(0, 1),
                                labels = c("Yes", "No")))) +
  labs(title = "Credit Score effect on Customer Churning",
       x = "Credit Score",
       y = "Churn Status")

From the visualization it is clear that the difference of mean for both samples is not large.

avg_credit_scores <- data |>
  group_by(Exited) |>
  summarize(avg_credit_score = mean(CreditScore)) |>
  arrange(Exited)

avg_credit_scores
## # A tibble: 2 × 2
##   Exited avg_credit_score
##    <dbl>            <dbl>
## 1      0             652.
## 2      1             645.

Neyman-Pearson Hypothesis Test

From this is clear that the difference in the means of both samples is very low, therefore it is vital to choose alpha and power level values. Since, mean difference is low, we need to choose values such that the Type I and Type II errors need to less. So, to limit the Type I error rate, a value smaller than the general value need to chosen, therefore alpha = 0.05. Similarly, as the power level indicates the probability at which we can correctly reject the null hypothesis, when it’s false, we need to choose a value large than the general value and close to 1, therefore power=0.9.

library(pwrss)

test <- pwrss.t.2means(mu1 = 100, 
                       sd1 = sd(pluck(data, "CreditScore")),
                       kappa = 1,
                       power = .9, alpha = 0.05, 
                       alternative = "not equal")
##  Difference between Two means 
##  (Independent Samples t Test) 
##  H0: mu1 = mu2 
##  HA: mu1 != mu2 
##  ------------------------------ 
##   Statistical power = 0.9 
##   n1 = 21 
##   n2 = 21 
##  ------------------------------ 
##  Alternative = "not equal" 
##  Degrees of freedom = 40 
##  Non-centrality parameter = 3.352 
##  Type I error rate = 0.05 
##  Type II error rate = 0.1
str(test)
## List of 6
##  $ parms:List of 12
##   ..$ mu1        : num 100
##   ..$ mu2        : num 0
##   ..$ sd1        : num 96.7
##   ..$ sd2        : num 96.7
##   ..$ kappa      : num 1
##   ..$ welch.df   : logi FALSE
##   ..$ paired     : logi FALSE
##   ..$ paired.r   : num 0.5
##   ..$ alpha      : num 0.05
##   ..$ margin     : num 0
##   ..$ alternative: chr "not equal"
##   ..$ verbose    : logi TRUE
##  $ test : chr "t"
##  $ df   : num 40
##  $ ncp  : num 3.35
##  $ power: num 0.9
##  $ n    : Named num [1:2] 21 21
##   ..- attr(*, "names")= chr [1:2] "n1" "n2"
##  - attr(*, "class")= chr [1:3] "pwrss" "t" "2means"
plot(test)

cohen.d(d = filter(data, Exited == 0) |> pluck("CreditScore"),
        f = filter(data, Exited == 1) |> pluck("CreditScore"))
## 
## Cohen's d
## 
## d estimate: 0.06682248 (negligible)
## 95 percent confidence interval:
##      lower      upper 
## 0.01815319 0.11549177

Here the d value is minimum, which indicates that the spread among the values difference is not significant, which means that it results in confident hypothesis.

Fisher’s Style Test

bootstrap <- function(x, func = mean, n_iter = 10000) {
  func_values <- numeric(n_iter)
  
  for (i in 1:n_iter) {
    x_sample <- sample(x, size = length(x), replace = TRUE)
    func_values[i] <- func(x_sample)
  }
  return(func_values)
}
diff_in_avg <- function(x_data) {
  avg_values <- mean(x_data, na.rm=TRUE)
  return(avg_values)
}

bootstrap_credit_score <- bootstrap(data$CreditScore, diff_in_avg, n_iter = 10000)
hist(bootstrap_credit_score, breaks = 30, col = "skyblue", 
     main = "Bootstrap Distribution of Credit Score",
     xlab = "CreditScore")

observed_diff_credit_score <- mean(data$CreditScore[data$Exited == 0]) - 
                         mean(data$CreditScore[data$Exited == 1])

p_value_credit_score <- mean(abs(bootstrap_credit_score) >= abs(observed_diff_credit_score))
p_value_credit_score
## [1] 1

Null Hypothesis 2:

Hypothesis: There is no significant difference in the distribution of ages between male and female customers.

data |>
  ggplot() +
  geom_boxplot(mapping = 
                 aes(x = Age, 
                     y = Gender)) +
  labs(x = "Age",
       y = "Gender")
## Warning: Removed 1 rows containing non-finite values (`stat_boxplot()`).

avg_ages <- data |>
  group_by(Gender) |>
  summarize(avg_age = mean(Age, na.rm=TRUE)) |>
  arrange(Gender)

avg_ages
## # A tibble: 2 × 2
##   Gender avg_age
##   <chr>    <dbl>
## 1 Female    39.2
## 2 Male      38.7

Similar to earlier hypothesis the difference is means is very low in this case as well, in fact much lower than the earlier value, therefore we need to choose much critical values for alpha and power to mitigate the Type I and Type II errors. So, alpha=0.01 and power=0.95, which such as small difference, and low Type I error rate (alpha) and high probability in reject a false null hypothesis (power) it could result in a hypothesis high good confidence. However, it can be only be verified upon calculating the d value.

Neyman-Pearson Hypothesis Test

test <- pwrss.t.2means(mu1 = 100, 
                       sd1 = sd(pluck(data, "Age"), na.rm=TRUE),
                       kappa = 1,
                       power = .9, alpha = 0.05, 
                       alternative = "not equal")
##  Difference between Two means 
##  (Independent Samples t Test) 
##  H0: mu1 = mu2 
##  HA: mu1 != mu2 
##  ------------------------------ 
##   Statistical power = 0.9 
##   n1 = 4 
##   n2 = 4 
##  ------------------------------ 
##  Alternative = "not equal" 
##  Degrees of freedom = 6 
##  Non-centrality parameter = 13.485 
##  Type I error rate = 0.05 
##  Type II error rate = 0.1
str(test)
## List of 6
##  $ parms:List of 12
##   ..$ mu1        : num 100
##   ..$ mu2        : num 0
##   ..$ sd1        : num 10.5
##   ..$ sd2        : num 10.5
##   ..$ kappa      : num 1
##   ..$ welch.df   : logi FALSE
##   ..$ paired     : logi FALSE
##   ..$ paired.r   : num 0.5
##   ..$ alpha      : num 0.05
##   ..$ margin     : num 0
##   ..$ alternative: chr "not equal"
##   ..$ verbose    : logi TRUE
##  $ test : chr "t"
##  $ df   : num 6
##  $ ncp  : num 13.5
##  $ power: num 0.9
##  $ n    : Named num [1:2] 4 4
##   ..- attr(*, "names")= chr [1:2] "n1" "n2"
##  - attr(*, "class")= chr [1:3] "pwrss" "t" "2means"
plot(test)

An interesting observation from the visualization is that the beta value is zero. Which essentially means that the rate at is which a Type II error is observed is almost 0, which is a result of choosing strict power value.

cohen.d(d = filter(data, Gender == "Female") |> pluck("Age"),
        f = filter(data, Gender == "Male") |> pluck("Age"), na.rm=TRUE)
## 
## Cohen's d
## 
## d estimate: 0.0548348 (negligible)
## 95 percent confidence interval:
##      lower      upper 
## 0.01546101 0.09420858

As anticipated the d value is also very low, which indicates that this a confident hypothesis, and therefore it not needs to be warned.

Fisher’s Style Test

bootstrap_age <- bootstrap(data$Age, diff_in_avg, n_iter = 10000)

hist(bootstrap_age, breaks = 30, col = "skyblue", 
     main = "Bootstrap Distribution of Age",
     xlab = "Age")

observed_diff_age <- mean(data$Age[data$Gender == "Female"], na.rm=TRUE) - 
                         mean(data$Age[data$Gender == "Male"], na.rm=TRUE)

p_value_age <- mean(abs(bootstrap_age) >= abs(observed_diff_age))
p_value_age
## [1] 1