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>
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.
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.
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
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.
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.
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