mydata <- read.table("C:/Users/eneja/OneDrive/Namizje/IMB R/archive/insurance.csv",
                     fill = TRUE,
                     header = TRUE,
                     sep = ",",
                     dec = ".")

Because we have a lot of units in the sample, we randomly picked 100 units.

set.seed(11)
n<-100
mydata1 <-mydata%>%
  sample_n(n)

head(mydata1)
##   age    sex    bmi children smoker    region   charges
## 1  63   male 28.310        0     no northwest 13770.098
## 2  26 female 40.185        0     no northwest  3201.245
## 3  62 female 25.000        0     no southwest 13451.122
## 4  29   male 29.735        2     no northwest 18157.876
## 5  44 female 25.000        1     no southwest  7623.518
## 6  18 female 37.290        1     no southeast  2219.445
str(mydata1)
## 'data.frame':    100 obs. of  7 variables:
##  $ age     : int  63 26 62 29 44 18 26 63 34 25 ...
##  $ sex     : chr  "male" "female" "female" "male" ...
##  $ bmi     : num  28.3 40.2 25 29.7 25 ...
##  $ children: int  0 0 0 2 1 1 2 0 1 4 ...
##  $ smoker  : chr  "no" "no" "no" "no" ...
##  $ region  : chr  "northwest" "northwest" "southwest" "northwest" ...
##  $ charges : num  13770 3201 13451 18158 7624 ...

Description:

Unit of observation: Person

Sample size: 1338 people later we sampled 100 people

Source: https://www.kaggle.com/datasets/mirichoi0218/insurance?resource=download

Hypothesis about the population arithmetic mean

Average BMI in slovenia in 2016 was 26.75 based on: https://www.euronews.com/2019/05/09/which-country-has-the-highest-average-bmi-in-europe?fbclid=IwAR0YAPvbeLfxQ_cXNfiIhFwJwNVGSN3pjCo47YHdIv7zYzUHgpnhLIvn48g

How did average BMI changed in 2018 compared to 2016?

H0:µ= 26.75

H1:µ≠ 26.75

ggplot(mydata1, aes(x=bmi))+
  geom_histogram(binwidth = 1, colour = "white", fill = "skyblue")+
  ylab("Frequency")+
  xlab("BMI")

shapiro.test(mydata1$bmi)
## 
##  Shapiro-Wilk normality test
## 
## data:  mydata1$bmi
## W = 0.98547, p-value = 0.3434

H0: Variable is normally distributed.

H1: Variable is not normally distributed.

Based on p-value we can’t reject null hypothesis. So we can say that distribution is normal.

mean(mydata1$bmi)
## [1] 30.01695
sd(mydata1$bmi)
## [1] 5.96881

Confidence interval

ybar = mean(mydata1$bmi); sd= sd(mydata1$bmi); n = nrow(mydata1)

se = sd/sqrt(n)

ybar_lower_t = ybar + qt(0.025, df=n-1)*se
ybar_upper_t = ybar + qt(0.975, df=n-1)*se

This upper and lower boundary doesn’t include our arithmetic mean.

mydata1 %>%
  identify_outliers(bmi)
## [1] age        sex        bmi        children   smoker     region     charges    is.outlier is.extreme
## <0 rows> (or 0-length row.names)
t.test(mydata1$bmi)
## 
##  One Sample t-test
## 
## data:  mydata1$bmi
## t = 50.29, df = 99, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  28.83261 31.20129
## sample estimates:
## mean of x 
##  30.01695

We are prepared to reject null hypothesis at p<0.001. Therefore we can say that people in 2018 have higher BMI score compared to 2016.

effectsize::cohens_d(mydata1$bmi, mu = 26.75)
## Cohen's d |       95% CI
## ------------------------
## 0.55      | [0.34, 0.76]
## 
## - Deviation from a difference of 26.75.
interpret_cohens_d(0.55, rules = "sawilowsky2009")
## [1] "medium"
## (Rules: sawilowsky2009)

Based on sample data, we found that the average BMI score in 2018 was 30.02 and it has increased compared to average in 2016 (p < 0.001, d = 0.55- medium sized effect).

HYPOTHESIS ABOUT THE DIFFERENCE BETWEEN TWO POPULATION AEITHMETIC MEANS

Research questions: Was there a difference in BMI score between Female and Male in 2018?

Hypothesis: Males and Females have different BMI score on average.

H0: µ male = µ female or µ male- µ female=0

H1: µ male ≠ µ female or µ male- µ female≠0

mydata1$SexF <- factor(mydata1$sex,
                              levels = c( "female", "male"),
                              labels = c("F","M"))
describeBy(mydata1$bmi,mydata1$SexF)
## 
##  Descriptive statistics by group 
## group: F
##    vars  n  mean   sd median trimmed  mad   min   max range skew kurtosis   se
## X1    1 47 28.01 5.66  26.89   27.81 5.95 17.29 40.19  22.9 0.37    -0.65 0.83
## ------------------------------------------------------------------------------------------ 
## group: M
##    vars  n mean   sd median trimmed  mad   min   max range skew kurtosis   se
## X1    1 53 31.8 5.71  31.35    31.7 5.71 18.34 43.01 24.67 0.08    -0.46 0.78

Based on means for groups of females and males we can say, that women have better BMI score on average.

First we are going to use ggplot to check distribution.

Female <- ggplot(mydata1[mydata1$SexF=="F", ], aes(x = bmi)) +
  theme_linedraw() +
  geom_histogram() +
  ggtitle("Female")

Male <- ggplot(mydata1[mydata1$SexF=="M", ], aes(x = bmi)) +
  theme_linedraw() +
  geom_histogram() +
  ggtitle("Male")

library(ggpubr)
ggarrange(Female, Male, ncol = 2, nrow = 1)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Based on histogram we can say that distribution looks normal.

library(ggpubr)
ggqqplot(mydata1,
         "bmi",
         facet.by = "SexF")
## Warning: The following aesthetics were dropped during statistical transformation: sample
## ℹ This can happen when ggplot fails to infer the correct grouping structure in the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical variable into a factor?
## The following aesthetics were dropped during statistical transformation: sample
## ℹ This can happen when ggplot fails to infer the correct grouping structure in the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical variable into a factor?
## The following aesthetics were dropped during statistical transformation: sample
## ℹ This can happen when ggplot fails to infer the correct grouping structure in the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical variable into a factor?
## The following aesthetics were dropped during statistical transformation: sample
## ℹ This can happen when ggplot fails to infer the correct grouping structure in the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical variable into a factor?

We also did Quantile quantile plot to check if variable is normally distributed.

mydata1 %>%
  group_by(SexF) %>%
  shapiro_test(bmi)
## # A tibble: 2 × 4
##   SexF  variable statistic     p
##   <fct> <chr>        <dbl> <dbl>
## 1 F     bmi          0.971 0.299
## 2 M     bmi          0.983 0.635

H0:Distribution is normal

H1:Distribution is not normal

Based on shapiro test (p>10%) we can say that distribution is normal.

mydata1 %>%
  group_by(SexF) %>%
  identify_outliers(bmi)
##  [1] SexF       age        sex        bmi        children   smoker     region     charges    is.outlier is.extreme
## <0 rows> (or 0-length row.names)
t.test(mydata1$bmi ~ mydata1$SexF, 
       paired = FALSE,
       var.equal = FALSE,
       alternative = "two.sided")
## 
##  Welch Two Sample t-test
## 
## data:  mydata1$bmi by mydata1$SexF
## t = -3.3262, df = 96.749, p-value = 0.001245
## alternative hypothesis: true difference in means between group F and group M is not equal to 0
## 95 percent confidence interval:
##  -6.048826 -1.527782
## sample estimates:
## mean in group F mean in group M 
##        28.00915        31.79745

Based on p-value we can reject null hypothesis. We can say that there are differences in BMI score for Female and Male.

library(effectsize)
effectsize::cohens_d(mydata1$bmi ~ mydata1$SexF,
                     pooled_sd = FALSE)
## Cohen's d |         95% CI
## --------------------------
## -0.67     | [-1.07, -0.26]
## 
## - Estimated using un-pooled SD.
interpret_cohens_d(0.03, rules = "sawilowsky2009")
## [1] "tiny"
## (Rules: sawilowsky2009)

Effect size tells us how big those differences are. In this case we have tiny differences.

Using the sample data we can say that there were differences in BMI scores for Female and Male in 2018 (the effect size is tiny, r=0.03). Women had better BMI score on average.

One way analysis of variance - ANOVA

We want to find out if number of children and BMI are related. Is average BMI different for people with different number of children in 2018?

H0: µ (0 children) = µ (1 child) = µ (2 children)= µ (3 children) = µ (4 children)

H1: at least one is different

head(mydata1[order(decreasing= TRUE, mydata1$children),], 5)
##    age  sex    bmi children smoker    region   charges SexF
## 29  41 male 29.640        5     no northeast  9222.403    M
## 10  25 male 33.660        4     no southeast  4504.662    M
## 43  40 male 30.875        4     no northwest  8162.716    M
## 68  51 male 24.415        4     no northwest 11520.100    M
## 71  33 male 29.400        4     no southwest  6059.173    M
mydata1$childrenF <- factor(mydata1$children,
                            levels = c(0, 1, 2, 3, 4),
                            labels = c("No children", "1 child", "2 children", "3 children", "4 children"))

We have only one person with 5 children so we remove it:

mydata2 <- mydata1[-c(29), ]
mydata2 %>%
  group_by(childrenF) %>%
  identify_outliers(bmi)
## # A tibble: 1 × 11
##   childrenF    age sex      bmi children smoker region    charges SexF  is.outlier is.extreme
##   <fct>      <int> <chr>  <dbl>    <int> <chr>  <chr>       <dbl> <fct> <lgl>      <lgl>     
## 1 2 children    37 female  17.3        2 no     northeast   6878. F     TRUE       FALSE

The assumption of having no outliers has been violated. Therefore we are going to use nonparametric test

library(psych)
psych::describe(mydata2$bmi)
##    vars  n  mean sd median trimmed  mad   min   max range skew kurtosis  se
## X1    1 99 30.02  6  29.74   29.86 6.48 17.29 43.01 25.72 0.19    -0.61 0.6

General arithmetic mean is 30.02.

psych::describeBy(x = mydata2$bmi, group = mydata2$childrenF)
## 
##  Descriptive statistics by group 
## group: No children
##    vars  n  mean   sd median trimmed  mad   min   max range skew kurtosis   se
## X1    1 37 29.68 6.13  28.31   29.41 5.63 18.34 43.01 24.67  0.4    -0.77 1.01
## ------------------------------------------------------------------------------------------ 
## group: 1 child
##    vars  n  mean   sd median trimmed  mad  min  max range skew kurtosis   se
## X1    1 25 31.05 7.16  29.15   31.05 6.46 18.5 42.9  24.4 0.16     -1.2 1.43
## ------------------------------------------------------------------------------------------ 
## group: 2 children
##    vars  n  mean   sd median trimmed  mad   min   max range  skew kurtosis   se
## X1    1 19 30.36 5.38   32.3   30.67 3.85 17.29 38.09  20.8 -0.73    -0.28 1.23
## ------------------------------------------------------------------------------------------ 
## group: 3 children
##    vars  n  mean   sd median trimmed  mad   min   max range  skew kurtosis   se
## X1    1 13 28.35 4.91  29.83   28.38 6.62 20.05 36.29 16.24 -0.03    -1.15 1.36
## ------------------------------------------------------------------------------------------ 
## group: 4 children
##    vars n  mean   sd median trimmed  mad   min  max range  skew kurtosis   se
## X1    1 5 30.49 3.92  30.88   30.49 4.13 24.41 34.1  9.69 -0.49    -1.61 1.75

First assumption that we must test is Lavene test for homogeneity of variances.

library(car)
leveneTest(mydata2$bmi, group = mydata2$childrenF)
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  4  1.3951 0.2417
##       94

H0: (homogeneity) Variances of BMI are the same in all groups

H1: (heterogeneity) Variances of BMI are not the same in all groups

Based on Levene test we cannot reject H0(p>10%). Therefore, we can assume homogeneity.

mydata2 %>%
  group_by(childrenF) %>%
  shapiro_test(bmi)
## # A tibble: 5 × 4
##   childrenF   variable statistic     p
##   <fct>       <chr>        <dbl> <dbl>
## 1 No children bmi          0.962 0.228
## 2 1 child     bmi          0.949 0.243
## 3 2 children  bmi          0.931 0.184
## 4 3 children  bmi          0.956 0.691
## 5 4 children  bmi          0.909 0.460

H0: BMI is normally distributed

H1: BMI is not normally distributed

Based on Shapiro test of normality we cannot reject null hypothesis (p>10%). We can assume BMI is normally distributed.

Because assumption of no outliers was violated we are going to perform nonparametric test.

H0: Distribution locations of BMI are the same

H1: Distribution locations of BMI are not the same

kruskal.test(bmi ~ childrenF, 
             data = mydata2)
## 
##  Kruskal-Wallis rank sum test
## 
## data:  bmi by childrenF
## Kruskal-Wallis chi-squared = 1.9024, df = 4, p-value = 0.7537

Based on p-value we cannot reject null hypothesis. Therefore, we can say that number of children does not affect BMI score of a person.

kruskal_effsize(bmi ~ childrenF, 
                data = mydata2)
## # A tibble: 1 × 5
##   .y.       n effsize method  magnitude
## * <chr> <int>   <dbl> <chr>   <ord>    
## 1 bmi      99 -0.0223 eta2[H] small
groups_nonpar <- wilcox_test(bmi ~ childrenF,
                             paired = FALSE,
                             p.adjust.method = "bonferroni",
                             data = mydata2)

groups_nonpar
## # A tibble: 10 × 9
##    .y.   group1      group2        n1    n2 statistic     p p.adj p.adj.signif
##  * <chr> <chr>       <chr>      <int> <int>     <dbl> <dbl> <dbl> <chr>       
##  1 bmi   No children 1 child       37    25      402. 0.393     1 ns          
##  2 bmi   No children 2 children    37    19      312  0.5       1 ns          
##  3 bmi   No children 3 children    37    13      268  0.55      1 ns          
##  4 bmi   No children 4 children    37     5       84  0.756     1 ns          
##  5 bmi   1 child     2 children    25    19      237  1         1 ns          
##  6 bmi   1 child     3 children    25    13      200  0.259     1 ns          
##  7 bmi   1 child     4 children    25     5       64  0.957     1 ns          
##  8 bmi   2 children  3 children    19    13      149  0.337     1 ns          
##  9 bmi   2 children  4 children    19     5       47  1         1 ns          
## 10 bmi   3 children  4 children    13     5       25  0.503     1 ns

After correction we can see that p values for every possible pair are not statistically significant or they are statistically different.

pwc <- mydata2 %>%
  pairwise_t_test(bmi ~ childrenF, 
                  paired = FALSE,
                  p.adjust.method = "bonferroni")

Kruskal_results <- kruskal_test(bmi ~ childrenF,
                                data = mydata2)

library(rstatix)
pwc <- pwc %>% 
       add_y_position(fun = "median", step.increase = 0.35)

library(ggpubr)
ggboxplot(mydata2, x = "childrenF", y = "bmi", add = "point", ylim=c(15, 50)) +
  stat_pvalue_manual(pwc, hide.ns = FALSE) +
  stat_summary(fun = mean, geom = "point", shape = 16, size = 4,
               aes(group = childrenF), color = "darkred",
               position = position_dodge(width = 0.8)) +
  stat_summary(fun = mean, colour = "darkred", 
               position = position_dodge(width = 0.8),
               geom = "text", vjust = -0.5, hjust = -1,
               aes(label = round(after_stat(y), digits = 2), group = childrenF)) +
  labs(subtitle = get_test_label(Kruskal_results,  detailed = TRUE),
       caption = get_pwc_label(pwc))

Based on sample data we didn’t find statistically significant relationship between BMI score and number of children in 2018 (χ2=1.9, p = 0.75, n=99).