library(tidyverse)
## Warning: package 'tidyverse' was built under R version 4.1.1
## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --
## v ggplot2 3.4.2     v purrr   1.0.1
## v tibble  3.2.1     v dplyr   1.1.2
## v tidyr   1.3.0     v stringr 1.5.0
## v readr   2.1.2     v forcats 0.5.1
## Warning: package 'ggplot2' was built under R version 4.1.3
## Warning: package 'tibble' was built under R version 4.1.3
## Warning: package 'tidyr' was built under R version 4.1.3
## Warning: package 'readr' was built under R version 4.1.3
## Warning: package 'purrr' was built under R version 4.1.3
## Warning: package 'dplyr' was built under R version 4.1.3
## Warning: package 'stringr' was built under R version 4.1.3
## Warning: package 'forcats' was built under R version 4.1.1
## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
library(ggthemes)
## Warning: package 'ggthemes' was built under R version 4.1.1
data <- read.csv("https://vincentarelbundock.github.io/Rdatasets/csv/AER/DoctorVisits.csv")

# Q1
# The data described as: Cross-section data originating 
# from the 1977-1978 Australian Health Survey.

# There's a structure of the data, with all the variables and 
# the relevant information regarding the variables:

# data-description: https://vincentarelbundock.github.io/Rdatasets/doc/AER/DoctorVisits.html
# data is stored here: https://vincentarelbundock.github.io/Rdatasets/csv/AER/DoctorVisits.csv

str(data)
## 'data.frame':    5190 obs. of  13 variables:
##  $ rownames : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ visits   : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ gender   : chr  "female" "female" "male" "male" ...
##  $ age      : num  0.19 0.19 0.19 0.19 0.19 0.19 0.19 0.19 0.19 0.19 ...
##  $ income   : num  0.55 0.45 0.9 0.15 0.45 0.35 0.55 0.15 0.65 0.15 ...
##  $ illness  : int  1 1 3 1 2 5 4 3 2 1 ...
##  $ reduced  : int  4 2 0 0 5 1 0 0 0 0 ...
##  $ health   : int  1 1 0 0 1 9 2 6 5 0 ...
##  $ private  : chr  "yes" "yes" "no" "no" ...
##  $ freepoor : chr  "no" "no" "no" "no" ...
##  $ freerepat: chr  "no" "no" "no" "no" ...
##  $ nchronic : chr  "no" "no" "no" "no" ...
##  $ lchronic : chr  "no" "no" "no" "no" ...
# Q2
# With reproducible seed number:

vec_q2 <- c() #Empty Vector Generation
set.seed(2024)
sample_size_n <- floor((28/100)*nrow(data)) #Generation of sample size
for(i in 1:14000) {
  rows_for_sample <- sample(x = seq(1,nrow(data)), size = sample_size_n, replace = F)
  data_for_q2 <- data[rows_for_sample,c("income","age")]
  cor_value_i <- cor(data_for_q2$income, data_for_q2$age)
  vec_q2 <- c(vec_q2, cor_value_i)
  i = i+1
}

mean(vec_q2)
## [1] -0.2712032
paste("mean value is of the vector is: ", round(mean(vec_q2), 4)) %>% print(quote = F)
## [1] mean value is of the vector is:  -0.2712
sd(vec_q2)
## [1] 0.01887985
paste("Standard deviation of the vector is: ", 
      round(sd(vec_q2), 4)) %>% print(quote = F)
## [1] Standard deviation of the vector is:  0.0189
# Plotting the vector - with a density plot:
vec_q2 <- data.frame(vec_q2)
ggplot(vec_q2, aes(x = vec_q2)) +
  geom_density(fill = "purple1", size = 1.1) +
  theme(plot.title = element_text(hjust=0.5), 
        plot.background = element_rect(fill = "gray")) + 
  labs(title = "Density Plot of Vec_q2 values", 
     x = "Correlation", y = "Density")
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## i Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

# Q3

# I'll split the variable age based upon it's median value, so I'll get 
# about 50% for each part; In this case, the age is 32.

median(data$age)
## [1] 0.32
median_age <- median(data$age)
data$age_level <-
  ifelse(data$age < median_age, "low", "high")
table("Above/Below Median age(32)"=data$age_level)
## Above/Below Median age(32)
## high  low 
## 2702 2488
table("Above/Below Median age(32), in %"=data$age_level) %>% prop.table() %>% round(4)
## Above/Below Median age(32), in %
##   high    low 
## 0.5206 0.4794
# fills is a Color vector for the following histogram:
fills <- c(rep("purple1",10), rep("orange1",10))

ggplot(data, aes(x = income)) +
  geom_histogram(fill = fills, col = "black", bins = 10) +
  facet_grid(~age_level) +
  labs(title = "Income by Level of Low/High Age", x = "Age") + 
  theme(plot.title = element_text(hjust=0.5), 
        plot.background = element_rect(fill = "gray"))

# Q4

error_rate_q4 <- (min(table(data$age_level)) / nrow(data)) * 100

paste("percentage of error replies from the mode:", 
      round(error_rate_q4, 2), "%") %>% print(quote = F)
## [1] percentage of error replies from the mode: 47.94 %
# Q5
# I'll generate a variable based upon the median of the income variable.

data$income_above_median <- 
  ifelse(data$income <= median(data$income), 1, 0)

table("income's median"=data$income_above_median, 
      "age's median"=data$age_level) %>% prop.table() %>% round(4)
##                age's median
## income's median   high    low
##               0 0.1865 0.2418
##               1 0.3341 0.2376
# now testing the chi-square test:

table("income's median"=data$income_above_median, 
      "age's median"=data$age_level) %>% prop.table() %>% chisq.test()
## Warning in chisq.test(.): Chi-squared approximation may be incorrect
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  .
## X-squared = 2.3296e-32, df = 1, p-value = 1
print("according to the results, there's no significant dependency bewtween the variables.", quote = F)
## [1] according to the results, there's no significant dependency bewtween the variables.
# Q6
# I've decided to tag as NA the lowest 5% of income.

# This is the table before the removal 
table(data$income)
## 
##    0 0.01 0.06 0.15 0.25 0.35 0.45 0.55 0.65 0.75  0.9  1.1  1.3  1.5 
##   79   35   80  249 1195  462  400  467  455  441  589  361  162  215
# I'll copy the data to a new dataframe called data1 and remove the NA's

data1 <- data

data1$income <- 
  ifelse(data1$income <= quantile(data1$income, 0.05), NA, data1$income)

is.na(data1$income) %>% table()
## .
## FALSE  TRUE 
##  4747   443
# Number of new NA's: 443, which is 8.5% of the orignal data, but that's the data.
# According to the original data - the numbers match:

table(data$income <= 0.15)
## 
## FALSE  TRUE 
##  4747   443
# Removing data1's NA's:
data1 <- filter(data1, is.na(income) == F)
nrow(data1)
## [1] 4747
# As for the variances:
print(paste("Old variance is ",var(data$income)), quote = F)
## [1] Old variance is  0.13609215210329
print(paste("New variance is ",var(data1$income)), quote = F)
## [1] New variance is  0.124186551250547
print("The new value is lower than the old value. it makes sense because the new data has LESS extreme data so the variance is low", quote = F)
## [1] The new value is lower than the old value. it makes sense because the new data has LESS extreme data so the variance is low
# Q7; I'll continue my work on data1 - the trimmed data:

# First transformation: age_transform1: multiplication in 100 - according to the variable's description:
data1$age_transform1 <- data1$age * 100
# Second Transformation:
# Preparing normal data as requested
vec_norm <- rnorm(n = nrow(data1), mean = 0, sd = sd(data1$age))
vec_norm <- data.frame(vec_norm)

data1$age_transform2 <- data1$age_transform1 + vec_norm

data1$age_transform2 <- data1$age_transform2$vec_norm

# Original data with transform1:
cor(data1$age, data1$age_transform1)
## [1] 1
# Original data with transform2:
cor(data1$age, data1$age_transform2)
## [1] 0.9999494
# Conclusion:
#############

# The first correlation is exactly 1 because the transformation was a perfect linear transformation.
# The second transformation is a bit less than 1 because it wasn't a perfect linear transformation
# due to the normal value, thus the linear correlation isn't perfection which means it's lower than 1.


# Q8

y <- rexp(14000)
y <- data.frame(y)

ggplot(y, aes(x = y)) +
  geom_histogram(fill = "purple1", color = "black", size = 0.8, bins = 40) +
  ggtitle("Histogram of the variable y") + 
  theme(plot.title = element_text(hjust=0.5), 
        plot.background = element_rect(fill = "gray"))

# Based upon the shape of the distribution, it seems rightly/positively skewed.

means <- c()
for(i in 1:14000) {
  relevant_rows <- ceiling(runif(n = 5, min = 1, max = 14000))
  results_q8 <- y[relevant_rows,"y"]
  mean_value <- mean(results_q8)
  means <- c(means, mean_value)
  i = i+1
}
means <- data.frame(means)
ggplot(means, aes(x = means)) +
  geom_histogram(fill = "purple", color = "black", size = 1.05) +
  ggtitle("Histogram of the means (for sample size = 5)") + 
  theme(plot.title = element_text(hjust=0.5), 
        plot.background = element_rect(fill = "gray"))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

# In this case, for sample-size = 5, the distribution is positively/rightly skewed.


means <- c()
for(i in 1:14000) {
  relevant_rows <- ceiling(runif(n = 80, min = 1, max = 14000))
  results_q8 <- y[relevant_rows,"y"]
  mean_value <- mean(results_q8)
  means <- c(means, mean_value)
  i = i+1
}
means <- data.frame(means)
ggplot(means, aes(x = means)) +
  geom_histogram(fill = "purple", color = "black", size = 1.05) +
  ggtitle("Histogram of the means (for sample size = 80)") + 
  theme(plot.title = element_text(hjust=0.5), 
        plot.background = element_rect(fill = "gray"))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

# In this case, for sample-size = 80, the shape of the distribution, tends
# to approximate to a normal (symmetric) distribution.


# Q9:
###################################
# Note: the date I'm using for is #
# the original mean of the vector #
# y - that means: before sampling #
###################################


# For sample-size = 5

means <- c()
checkups <- c()
mean_value_of_y <- mean(y$y)
for(i in 1:14000) {
  relevant_rows <- ceiling(runif(n = 5, min = 1, max = 14000))
  data_q8 <- y[relevant_rows,"y"]
  mean_value <- mean(data_q8)
  diag_value <- ifelse(mean_value > mean_value_of_y, "bigger", "smaller")
  means <- c(means, mean_value)
  checkups <- c(checkups, diag_value)
  i = i+1
}
q9_df <- data.frame(means, checkups)
table(q9_df$checkups)
## 
##  bigger smaller 
##    6157    7843
# For sample-size = 5

means <- c()
checkups <- c()
mean_value_of_y <- mean(y$y)
for(i in 1:14000) {
  relevant_rows <- ceiling(runif(n = 80, min = 1, max = 14000))
  data_q8 <- y[relevant_rows,"y"]
  mean_value <- mean(data_q8)
  diag_value <- ifelse(mean_value > mean_value_of_y, "bigger", "smaller")
  means <- c(means, mean_value)
  checkups <- c(checkups, diag_value)
  i = i+1
}
q9_df <- data.frame(means, checkups)
table(q9_df$checkups)
## 
##  bigger smaller 
##    6746    7254
# We have more "Smaller" than "Bigger", both for sample-size = 80 and for sample size = 5.
# That makes sense because of the rightly-skewed shape of the initial exponential
# distribution - the mean is greater than the median, so there are less than
# 50% cases greater than the mean, and more than 50% cases lower than the mean,
# and that's exactly what can be seen.





# Q10
# Analysis 1:
############# 
# Difference in number of visits between patients that hold private insurance
# and patients that does not hold such an insurance (dependent variable: visits, 
# independent variable: private)

# Difference of variances:
var.test(visits ~ private, data)
## 
##  F test to compare two variances
## 
## data:  visits by private
## F = 1.123, num df = 2891, denom df = 2297, p-value = 0.00344
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  1.039075 1.213265
## sample estimates:
## ratio of variances 
##           1.122993
# Conclusion:
# Variance test indicated there's a significant difference between the two groups.

t.test(visits ~ private, data, var.equal = F)
## 
##  Welch Two Sample t-test
## 
## data:  visits by private
## t = 0.57747, df = 5038, p-value = 0.5636
## alternative hypothesis: true difference in means between group no and group yes is not equal to 0
## 95 percent confidence interval:
##  -0.03064385  0.05623529
## sample estimates:
##  mean in group no mean in group yes 
##         0.3073997         0.2946040
# Conclusion:
# According to the t-test, there's no significant difference between means from
# private = yes and private = no.




col_vector <- c(rep(times = 10, "orange1"), rep(times = 9, "purple1"))

ggplot(data, aes(x = visits)) + 
  geom_bar(fill = col_vector, color = "black") + 
  facet_wrap(~private) + 
  ggtitle("Number of visits per patients 
          with/without private insurance") + 
  theme(plot.title = element_text(hjust=0.5), 
        plot.background = element_rect(fill = "gray"))

# Analysis 2: Is there a correlation between age and income?
# There's the analysis:

cor.test(data$age, data$income)
## 
##  Pearson's product-moment correlation
## 
## data:  data$age and data$income
## t = -20.284, df = 5188, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.2960968 -0.2456782
## sample estimates:
##        cor 
## -0.2710734
# According to the results, there's a significant negative correlation between
# the variables, which means that an increase in one variable means a decrease
# in the other variable and vice verse, nonetheless, the association between
# the variables is rather week, yet it is significant.

# There's the relevant plot:

ggplot(data, aes(x = income, y = age)) + geom_point(color = "red", size = 2.1) + 
  ggtitle("Correlation between income and age") + 
  theme(plot.title = element_text(hjust=0.5), 
        plot.background = element_rect(fill = "gray"))

# Q11:
# Correlation plot between age and income

data_q11 <- data[,c("visits", "illness", "age", "income")]

library(corrplot)
## Warning: package 'corrplot' was built under R version 4.1.2
## corrplot 0.92 loaded
cor(data_q11) %>% round(4)
##          visits illness     age  income
## visits   1.0000  0.2236  0.1245 -0.0768
## illness  0.2236  1.0000  0.2050 -0.1488
## age      0.1245  0.2050  1.0000 -0.2711
## income  -0.0768 -0.1488 -0.2711  1.0000
corrplot(cor(data_q11))

## WRITE INTERPRETATION about this:


table(data$gender) %>% knitr::kable()
Var1 Freq
female 2702
male 2488 
# Boxplot of age by gender, using plotly package:
library(plotly)
## Warning: package 'plotly' was built under R version 4.1.3
## 
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## The following object is masked from 'package:stats':
## 
##     filter
## The following object is masked from 'package:graphics':
## 
##     layout
plot_ly(data, x = ~gender, y = ~age, type = "box", color = ~gender) %>%
  layout(title = "Age Distribution by Gender", xaxis = list(title = "Gender"), yaxis = list(title = "Age"))
## Warning in RColorBrewer::brewer.pal(N, "Set2"): minimal value for n is 3, returning requested palette with 3 different levels
## Warning in RColorBrewer::brewer.pal(N, "Set2"): minimal value for n is 3, returning requested palette with 3 different levels