library(tidyverse)
## Warning: package 'tidyverse' was built under R version 4.2.3
## Warning: package 'ggplot2' was built under R version 4.2.3
## Warning: package 'tibble' was built under R version 4.2.3
## Warning: package 'tidyr' was built under R version 4.2.3
## Warning: package 'readr' was built under R version 4.2.3
## Warning: package 'purrr' was built under R version 4.2.3
## Warning: package 'dplyr' was built under R version 4.2.3
## Warning: package 'stringr' was built under R version 4.2.3
## Warning: package 'forcats' was built under R version 4.2.3
## Warning: package 'lubridate' was built under R version 4.2.3
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.2     ✔ readr     2.1.4
## ✔ forcats   1.0.0     ✔ stringr   1.5.0
## ✔ ggplot2   3.4.2     ✔ tibble    3.2.0
## ✔ lubridate 1.9.2     ✔ tidyr     1.3.0
## ✔ purrr     1.0.1     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the ]8;;http://conflicted.r-lib.org/conflicted package]8;; to force all conflicts to become errors
library(ggplot2)
library(AER)
## Warning: package 'AER' was built under R version 4.2.3
## Loading required package: car
## Warning: package 'car' was built under R version 4.2.3
## Loading required package: carData
## Warning: package 'carData' was built under R version 4.2.3
## 
## Attaching package: 'car'
## 
## The following object is masked from 'package:dplyr':
## 
##     recode
## 
## The following object is masked from 'package:purrr':
## 
##     some
## 
## Loading required package: lmtest
## Warning: package 'lmtest' was built under R version 4.2.3
## Loading required package: zoo
## Warning: package 'zoo' was built under R version 4.2.3
## 
## Attaching package: 'zoo'
## 
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
## 
## Loading required package: sandwich
## Warning: package 'sandwich' was built under R version 4.2.3
## Loading required package: survival
## Warning: package 'survival' was built under R version 4.2.3
library(DT)
## Warning: package 'DT' was built under R version 4.2.3
data("NMES1988")
t <- NMES1988
datatable(t)
# 1. Phân tích mô tả

# So sánh số lần nhập viện trung bình giữa các nhóm bảo hiểm
t %>%
  group_by(insurance) %>%
  summarise(mean_hospital = mean(hospital, na.rm = TRUE))
## # A tibble: 2 × 2
##   insurance mean_hospital
##   <fct>             <dbl>
## 1 no                0.311
## 2 yes               0.292
t %>%
  group_by(medicaid) %>%
  summarise(mean_hospital = mean(hospital, na.rm = TRUE))
## # A tibble: 2 × 2
##   medicaid mean_hospital
##   <fct>            <dbl>
## 1 no               0.284
## 2 yes              0.415
# Biểu đồ so sánh
ggplot(t, aes(x = insurance, y = hospital)) +
  geom_boxplot() +
  labs(title = "Số lần nhập viện theo bảo hiểm cá nhân")

ggplot(t, aes(x = medicaid, y = hospital)) +
  geom_boxplot() +
  labs(title = "Số lần nhập viện theo Medicaid")

# 2. Phân tích thống kê

# T-test cho bảo hiểm tư nhân
t.test(hospital ~ insurance, data = t)
## 
##  Welch Two Sample t-test
## 
## data:  hospital by insurance
## t = 0.72203, df = 1664.9, p-value = 0.4704
## alternative hypothesis: true difference in means between group no and group yes is not equal to 0
## 95 percent confidence interval:
##  -0.03249703  0.07036169
## sample estimates:
##  mean in group no mean in group yes 
##         0.3106599         0.2917276
# T-test cho Medicaid
t.test(hospital ~ medicaid, data = t)
## 
##  Welch Two Sample t-test
## 
## data:  hospital by medicaid
## t = -3.0202, df = 464.63, p-value = 0.002665
## alternative hypothesis: true difference in means between group no and group yes is not equal to 0
## 95 percent confidence interval:
##  -0.21698777 -0.04592593
## sample estimates:
##  mean in group no mean in group yes 
##         0.2839660         0.4154229
# 3. Mô hình hồi quy

# Mô hình hồi quy đơn giản
model1 <- lm(hospital ~ insurance + medicaid, data = t)
summary(model1)
## 
## Call:
## lm(formula = hospital ~ insurance + medicaid, data = t)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.4420 -0.2890 -0.2890 -0.2577  7.7110 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   0.25770    0.02828   9.113  < 2e-16 ***
## insuranceyes  0.03130    0.03064   1.022 0.307033    
## medicaidyes   0.15297    0.04433   3.451 0.000564 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.7455 on 4403 degrees of freedom
## Multiple R-squared:  0.002809,   Adjusted R-squared:  0.002356 
## F-statistic: 6.201 on 2 and 4403 DF,  p-value: 0.002045
# Mô hình hồi quy đa biến
model2 <- lm(hospital ~ insurance + medicaid + age + chronic + health + income, data = t)
summary(model2)
## 
## Call:
## lm(formula = hospital ~ insurance + medicaid + age + chronic + 
##     health + income, data = t)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.0595 -0.3076 -0.1894 -0.0732  7.5696 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     -0.310744   0.131994  -2.354  0.01861 *  
## insuranceyes     0.052342   0.029768   1.758  0.07876 .  
## medicaidyes      0.060436   0.043182   1.400  0.16172    
## age              0.050090   0.017310   2.894  0.00383 ** 
## chronic          0.098685   0.008634  11.430  < 2e-16 ***
## healthpoor       0.302892   0.035157   8.616  < 2e-16 ***
## healthexcellent -0.083900   0.041162  -2.038  0.04158 *  
## income           0.002389   0.003780   0.632  0.52739    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.7185 on 4398 degrees of freedom
## Multiple R-squared:  0.07491,    Adjusted R-squared:  0.07343 
## F-statistic: 50.87 on 7 and 4398 DF,  p-value: < 2.2e-16
# Kiểm tra đa cộng tuyến
vif(model2)
##               GVIF Df GVIF^(1/(2*Df))
## insurance 1.312911  1        1.145823
## medicaid  1.319688  1        1.148777
## age       1.025862  1        1.012848
## chronic   1.158688  1        1.076424
## health    1.186497  2        1.043678
## income    1.043121  1        1.021333
# 4. Phân tích tương tác

# Tương tác giữa bảo hiểm và thu nhập
model3 <- lm(hospital ~ insurance * income + medicaid * income, data = t)
summary(model3)
## 
## Call:
## lm(formula = hospital ~ insurance * income + medicaid * income, 
##     data = t)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.5796 -0.2936 -0.2881 -0.2578  7.7136 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          2.579e-01  3.880e-02   6.647 3.35e-11 ***
## insuranceyes         3.969e-02  4.175e-02   0.951   0.3418    
## income              -8.776e-05  1.427e-02  -0.006   0.9951    
## medicaidyes          1.192e-01  5.792e-02   2.058   0.0396 *  
## insuranceyes:income -2.976e-03  1.479e-02  -0.201   0.8405    
## income:medicaidyes   2.889e-02  2.890e-02   1.000   0.3176    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.7456 on 4400 degrees of freedom
## Multiple R-squared:  0.003211,   Adjusted R-squared:  0.002078 
## F-statistic: 2.835 on 5 and 4400 DF,  p-value: 0.01465
# Tương tác giữa bảo hiểm và số bệnh mãn tính
model4 <- lm(hospital ~ insurance * chronic + medicaid * chronic, data = t)
summary(model4)
## 
## Call:
## lm(formula = hospital ~ insurance * chronic + medicaid * chronic, 
##     data = t)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.1358 -0.3536 -0.2233 -0.0929  7.6464 
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)           0.12735    0.03957   3.218  0.00130 ** 
## insuranceyes         -0.03444    0.04337  -0.794  0.42722    
## chronic               0.08841    0.01860   4.754 2.06e-06 ***
## medicaidyes          -0.04275    0.06914  -0.618  0.53637    
## insuranceyes:chronic  0.04195    0.02048   2.048  0.04059 *  
## chronic:medicaidyes   0.07701    0.02927   2.631  0.00854 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.7251 on 4400 degrees of freedom
## Multiple R-squared:  0.05725,    Adjusted R-squared:  0.05618 
## F-statistic: 53.44 on 5 and 4400 DF,  p-value: < 2.2e-16
# 5. Biểu đồ tương tác

# Biểu đồ tương tác giữa bảo hiểm và thu nhập
ggplot(t, aes(x = income, y = hospital, color = insurance)) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE) +
  labs(title = "Tương tác giữa bảo hiểm cá nhân và thu nhập")
## `geom_smooth()` using formula = 'y ~ x'

# Biểu đồ tương tác giữa bảo hiểm và số bệnh mãn tính
ggplot(t, aes(x = chronic, y = hospital, color = insurance)) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE) +
  labs(title = "Tương tác giữa bảo hiểm cá nhân và số bệnh mãn tính")
## `geom_smooth()` using formula = 'y ~ x'

# 1. So sánh số lần nhập viện trung bình giữa nhóm có và không có bảo hiểm cá nhân
insurance_comparison <- t %>%
  group_by(insurance) %>%
  summarise(
    mean_hospital = mean(hospital, na.rm = TRUE),
    sd_hospital = sd(hospital, na.rm = TRUE),
    n = n()
  )

print("So sánh số lần nhập viện theo bảo hiểm cá nhân:")
## [1] "So sánh số lần nhập viện theo bảo hiểm cá nhân:"
print(insurance_comparison)
## # A tibble: 2 × 4
##   insurance mean_hospital sd_hospital     n
##   <fct>             <dbl>       <dbl> <int>
## 1 no                0.311       0.716   985
## 2 yes               0.292       0.755  3421
# T-test cho bảo hiểm tư nhân
t_test_insurance <- t.test(hospital ~ insurance, data = t)
print(t_test_insurance)
## 
##  Welch Two Sample t-test
## 
## data:  hospital by insurance
## t = 0.72203, df = 1664.9, p-value = 0.4704
## alternative hypothesis: true difference in means between group no and group yes is not equal to 0
## 95 percent confidence interval:
##  -0.03249703  0.07036169
## sample estimates:
##  mean in group no mean in group yes 
##         0.3106599         0.2917276
# 2. So sánh số lần nhập viện trung bình giữa nhóm có và không có Medicaid
medicaid_comparison <- t %>%
  group_by(medicaid) %>%
  summarise(
    mean_hospital = mean(hospital, na.rm = TRUE),
    sd_hospital = sd(hospital, na.rm = TRUE),
    n = n()
  )

print("So sánh số lần nhập viện theo Medicaid:")
## [1] "So sánh số lần nhập viện theo Medicaid:"
print(medicaid_comparison)
## # A tibble: 2 × 4
##   medicaid mean_hospital sd_hospital     n
##   <fct>            <dbl>       <dbl> <int>
## 1 no               0.284       0.735  4004
## 2 yes              0.415       0.841   402
# T-test cho Medicaid
t_test_medicaid <- t.test(hospital ~ medicaid, data = t)
print(t_test_medicaid)
## 
##  Welch Two Sample t-test
## 
## data:  hospital by medicaid
## t = -3.0202, df = 464.63, p-value = 0.002665
## alternative hypothesis: true difference in means between group no and group yes is not equal to 0
## 95 percent confidence interval:
##  -0.21698777 -0.04592593
## sample estimates:
##  mean in group no mean in group yes 
##         0.2839660         0.4154229
# 3. Xem xét nhóm không có bất kỳ bảo hiểm nào
t <- t %>%
  mutate(no_insurance = ifelse(insurance == "no" & medicaid == "no", "khongBHCN", "coBHCN"))

no_insurance_comparison <- t %>%
  group_by(no_insurance) %>%
  summarise(
    mean_hospital = mean(hospital, na.rm = TRUE),
    sd_hospital = sd(hospital, na.rm = TRUE),
    n = n()
  )

print("So sánh số lần nhập viện cho nhóm không có bảo hiểm:")
## [1] "So sánh số lần nhập viện cho nhóm không có bảo hiểm:"
print(no_insurance_comparison)
## # A tibble: 2 × 4
##   no_insurance mean_hospital sd_hospital     n
##   <chr>                <dbl>       <dbl> <int>
## 1 coBHCN               0.303       0.763  3762
## 2 khongBHCN            0.255       0.641   644
# T-test cho nhóm không có bảo hiểm
t_test_no_insurance <- t.test(hospital ~ no_insurance, data = t)
print(t_test_no_insurance)
## 
##  Welch Two Sample t-test
## 
## data:  hospital by no_insurance
## t = 1.718, df = 982.64, p-value = 0.08612
## alternative hypothesis: true difference in means between group coBHCN and group khongBHCN is not equal to 0
## 95 percent confidence interval:
##  -0.00688189  0.10362573
## sample estimates:
##    mean in group coBHCN mean in group khongBHCN 
##               0.3030303               0.2546584
# Biểu đồ so sánh
ggplot(t, aes(x = insurance, y = hospital)) +
  geom_boxplot() +
  labs(title = "Số lần nhập viện theo bảo hiểm cá nhân")

ggplot(t, aes(x = medicaid, y = hospital)) +
  geom_boxplot() +
  labs(title = "Số lần nhập viện theo Medicaid")

ggplot(t, aes(x = no_insurance, y = hospital)) +
  geom_boxplot() +
  labs(title = "Số lần nhập viện theo tình trạng bảo hiểm")

# 1. Phân tích tương quan
correlation <- cor.test(t$chronic, t$hospital)
print("Tương quan giữa số bệnh mãn tính và số lần nhập viện:")
## [1] "Tương quan giữa số bệnh mãn tính và số lần nhập viện:"
print(correlation)
## 
##  Pearson's product-moment correlation
## 
## data:  t$chronic and t$hospital
## t = 15.938, df = 4404, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.2054118 0.2612516
## sample estimates:
##       cor 
## 0.2335242
# 2. Mô hình hồi quy tuyến tính đơn giản
model <- lm(hospital ~ chronic, data = t)
summary(model)
## 
## Call:
## lm(formula = hospital ~ chronic, data = t)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.1300 -0.3551 -0.2260 -0.0968  7.6449 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 0.096816   0.016604   5.831 5.91e-09 ***
## chronic     0.129148   0.008103  15.938  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.7258 on 4404 degrees of freedom
## Multiple R-squared:  0.05453,    Adjusted R-squared:  0.05432 
## F-statistic:   254 on 1 and 4404 DF,  p-value: < 2.2e-16
# 3. Biểu đồ phân tán
ggplot(t, aes(x = chronic, y = hospital)) +
  geom_point(alpha = 0.5) +
  geom_smooth(method = "lm", se = FALSE, color = "red") +
  labs(title = "Mối quan hệ giữa số bệnh mãn tính và số lần nhập viện",
       x = "Số bệnh mãn tính",
       y = "Số lần nhập viện")
## `geom_smooth()` using formula = 'y ~ x'

# 4. Kiểm tra giả định tuyến tính
residualPlot(model)

# 5. Phân tích theo nhóm
t %>%
  group_by(chronic) %>%
  summarise(
    mean_hospital = mean(hospital, na.rm = TRUE),
    sd_hospital = sd(hospital, na.rm = TRUE),
    n = n()
  ) %>%
  print(n = Inf)  # In tất cả các hàng
## # A tibble: 9 × 4
##   chronic mean_hospital sd_hospital     n
##     <int>         <dbl>       <dbl> <int>
## 1       0         0.127       0.432  1025
## 2       1         0.194       0.542  1498
## 3       2         0.379       0.851   968
## 4       3         0.467       0.950   525
## 5       4         0.609       0.985   220
## 6       5         0.764       1.23    127
## 7       6         0.706       1.55     34
## 8       7         2.17        1.47      6
## 9       8         1.33        2.31      3
# 6. Biểu đồ boxplot
ggplot(t, aes(x = factor(chronic), y = hospital)) +
  geom_boxplot() +
  labs(title = "Phân bố số lần nhập viện theo số bệnh mãn tính",
       x = "Số bệnh mãn tính",
       y = "Số lần nhập viện")

# 1. Mô hình hồi quy đa biến
model <- lm(hospital ~ chronic + insurance + medicaid, data = t)
summary(model)
## 
## Call:
## lm(formula = hospital ~ chronic + insurance + medicaid, data = t)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.1182 -0.3526 -0.2250 -0.0973  7.6474 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.067140   0.030089   2.231   0.0257 *  
## chronic      0.127609   0.008139  15.679   <2e-16 ***
## insuranceyes 0.030209   0.029819   1.013   0.3111    
## medicaidyes  0.094196   0.043307   2.175   0.0297 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.7256 on 4402 degrees of freedom
## Multiple R-squared:  0.05555,    Adjusted R-squared:  0.05491 
## F-statistic:  86.3 on 3 and 4402 DF,  p-value: < 2.2e-16
# 2. Kiểm tra đa cộng tuyến
vif(model)
##   chronic insurance  medicaid 
##  1.009509  1.291533  1.301277
# 3. Biểu đồ tương tác giữa số bệnh mãn tính và bảo hiểm
ggplot(t, aes(x = chronic, y = hospital, color = insurance)) +
  geom_point(alpha = 0.5) +
  geom_smooth(method = "lm", se = FALSE) +
  labs(title = "Mối quan hệ giữa số bệnh mãn tính và số lần nhập viện theo bảo hiểm tư nhân",
       x = "Số bệnh mãn tính", y = "Số lần nhập viện")
## `geom_smooth()` using formula = 'y ~ x'

# 4. Biểu đồ tương tác giữa số bệnh mãn tính và Medicaid
ggplot(t, aes(x = chronic, y = hospital, color = medicaid)) +
  geom_point(alpha = 0.5) +
  geom_smooth(method = "lm", se = FALSE) +
  labs(title = "Mối quan hệ giữa số bệnh mãn tính và số lần nhập viện theo Medicaid",
       x = "Số bệnh mãn tính", y = "Số lần nhập viện")
## `geom_smooth()` using formula = 'y ~ x'

# 5. Phân tích theo nhóm
t %>%
  group_by(insurance, medicaid) %>%
  summarise(
    mean_hospital = mean(hospital, na.rm = TRUE),
    mean_chronic = mean(chronic, na.rm = TRUE),
    n = n()
  ) %>%
  print(n = Inf)
## `summarise()` has grouped output by 'insurance'. You can override using the
## `.groups` argument.
## # A tibble: 4 × 5
## # Groups:   insurance [2]
##   insurance medicaid mean_hospital mean_chronic     n
##   <fct>     <fct>            <dbl>        <dbl> <int>
## 1 no        no               0.255         1.50   644
## 2 no        yes              0.416         1.95   341
## 3 yes       no               0.290         1.50  3360
## 4 yes       yes              0.410         1.98    61
# 6. Mô hình với tương tác
model_interaction <- lm(hospital ~ chronic * insurance * medicaid, data = t)
summary(model_interaction)
## 
## Call:
## lm(formula = hospital ~ chronic * insurance * medicaid, data = t)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.1488 -0.3555 -0.2233 -0.0911  7.6445 
## 
## Coefficients:
##                                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                       0.13399    0.04081   3.283  0.00103 ** 
## chronic                           0.08070    0.01948   4.141 3.52e-05 ***
## insuranceyes                     -0.04294    0.04505  -0.953  0.34056    
## medicaidyes                      -0.06636    0.07819  -0.849  0.39609    
## chronic:insuranceyes              0.05153    0.02172   2.373  0.01770 *  
## chronic:medicaidyes               0.09816    0.03382   2.902  0.00372 ** 
## insuranceyes:medicaidyes          0.08934    0.16816   0.531  0.59524    
## chronic:insuranceyes:medicaidyes -0.08126    0.06772  -1.200  0.23020    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.7251 on 4398 degrees of freedom
## Multiple R-squared:  0.05765,    Adjusted R-squared:  0.05615 
## F-statistic: 38.43 on 7 and 4398 DF,  p-value: < 2.2e-16
library(epitools)
## 
## Attaching package: 'epitools'
## The following object is masked from 'package:survival':
## 
##     ratetable
# Tạo biến nhị phân cho việc nhập viện
t <- t  %>%
  mutate(hospitalized = ifelse(hospital > 0, "Yes", "No"))

# 1. Tính RR và OR cho bảo hiểm tư nhân
table_insurance <- table(t$insurance, t$hospitalized)
rr_insurance <- riskratio(table_insurance)
or_insurance <- oddsratio(table_insurance)

print("Relative Risk và Odds Ratio cho bảo hiểm tư nhân:")
## [1] "Relative Risk và Odds Ratio cho bảo hiểm tư nhân:"
print(rr_insurance)
## $data
##        
##           No Yes Total
##   no     775 210   985
##   yes   2766 655  3421
##   Total 3541 865  4406
## 
## $measure
##      risk ratio with 95% C.I.
##        estimate    lower    upper
##   no  1.0000000       NA       NA
##   yes 0.8980596 0.782042 1.031289
## 
## $p.value
##      two-sided
##       midp.exact fisher.exact chi.square
##   no          NA           NA         NA
##   yes  0.1322232    0.1331772  0.1302475
## 
## $correction
## [1] FALSE
## 
## attr(,"method")
## [1] "Unconditional MLE & normal approximation (Wald) CI"
print(or_insurance)
## $data
##        
##           No Yes Total
##   no     775 210   985
##   yes   2766 655  3421
##   Total 3541 865  4406
## 
## $measure
##      odds ratio with 95% C.I.
##        estimate     lower    upper
##   no  1.0000000        NA       NA
##   yes 0.8736031 0.7345964 1.042055
## 
## $p.value
##      two-sided
##       midp.exact fisher.exact chi.square
##   no          NA           NA         NA
##   yes  0.1322232    0.1331772  0.1302475
## 
## $correction
## [1] FALSE
## 
## attr(,"method")
## [1] "median-unbiased estimate & mid-p exact CI"
# 2. Tính RR và OR cho Medicaid
table_medicaid <- table(t$medicaid, t$hospitalized)
rr_medicaid <- riskratio(table_medicaid)
or_medicaid <- oddsratio(table_medicaid)

print("Relative Risk và Odds Ratio cho Medicaid:")
## [1] "Relative Risk và Odds Ratio cho Medicaid:"
print(rr_medicaid)
## $data
##        
##           No Yes Total
##   no    3244 760  4004
##   yes    297 105   402
##   Total 3541 865  4406
## 
## $measure
##      risk ratio with 95% C.I.
##       estimate    lower    upper
##   no   1.00000       NA       NA
##   yes  1.37608 1.153519 1.641582
## 
## $p.value
##      two-sided
##         midp.exact fisher.exact   chi.square
##   no            NA           NA           NA
##   yes 0.0008763161 0.0009559402 0.0005928227
## 
## $correction
## [1] FALSE
## 
## attr(,"method")
## [1] "Unconditional MLE & normal approximation (Wald) CI"
print(or_medicaid)
## $data
##        
##           No Yes Total
##   no    3244 760  4004
##   yes    297 105   402
##   Total 3541 865  4406
## 
## $measure
##      odds ratio with 95% C.I.
##       estimate    lower    upper
##   no  1.000000       NA       NA
##   yes 1.510224 1.188013 1.906643
## 
## $p.value
##      two-sided
##         midp.exact fisher.exact   chi.square
##   no            NA           NA           NA
##   yes 0.0008763161 0.0009559402 0.0005928227
## 
## $correction
## [1] FALSE
## 
## attr(,"method")
## [1] "median-unbiased estimate & mid-p exact CI"
# 3. Tính RR và OR cho số bệnh mãn tính
# Chia số bệnh mãn tính thành hai nhóm: ít (<=2) và nhiều (>2)
t <- t %>%
  mutate(chronic_group = ifelse(chronic <= 2, "Low", "High"))

table_chronic <- table(t$chronic_group, t$hospitalized)
rr_chronic <- riskratio(table_chronic)
or_chronic <- oddsratio(table_chronic)

print("Relative Risk và Odds Ratio cho số bệnh mãn tính:")
## [1] "Relative Risk và Odds Ratio cho số bệnh mãn tính:"
print(rr_chronic)
## $data
##        
##           No Yes Total
##   High   608 307   915
##   Low   2933 558  3491
##   Total 3541 865  4406
## 
## $measure
##       risk ratio with 95% C.I.
##         estimate     lower     upper
##   High 1.0000000        NA        NA
##   Low  0.4763949 0.4230582 0.5364558
## 
## $p.value
##       two-sided
##        midp.exact fisher.exact   chi.square
##   High         NA           NA           NA
##   Low           0 7.777527e-30 1.068538e-32
## 
## $correction
## [1] FALSE
## 
## attr(,"method")
## [1] "Unconditional MLE & normal approximation (Wald) CI"
print(or_chronic)
## $data
##        
##           No Yes Total
##   High   608 307   915
##   Low   2933 558  3491
##   Total 3541 865  4406
## 
## $measure
##       odds ratio with 95% C.I.
##         estimate     lower     upper
##   High 1.0000000        NA        NA
##   Low  0.3768339 0.3198127 0.4444088
## 
## $p.value
##       two-sided
##        midp.exact fisher.exact   chi.square
##   High         NA           NA           NA
##   Low           0 7.777527e-30 1.068538e-32
## 
## $correction
## [1] FALSE
## 
## attr(,"method")
## [1] "median-unbiased estimate & mid-p exact CI"

Bảo hiểm tư nhân:

Relative Risk (RR): 0.8980596 (95% CI: 0.782042 - 1.031289) Odds Ratio (OR): 0.8736031 (95% CI: 0.7345964 - 1.042055)

RR và OR đều nhỏ hơn 1, cho thấy có xu hướng giảm nhẹ khả năng nhập viện ở những người có bảo hiểm tư nhân. Tuy nhiên, khoảng tin cậy 95% của cả RR và OR đều chứa 1, và p-value > 0.05, nên sự khác biệt này không có ý nghĩa thống kê.

Medicaid:

Relative Risk (RR): 1.37608 (95% CI: 1.153519 - 1.641582) Odds Ratio (OR): 1.510224 (95% CI: 1.188013 - 1.906643)

RR và OR đều lớn hơn 1 và khoảng tin cậy 95% không chứa 1. p-value < 0.001, cho thấy sự khác biệt có ý nghĩa thống kê. Người có Medicaid có nguy cơ nhập viện cao hơn 37.6% (RR) và có khả năng nhập viện cao hơn 51% (OR) so với người không có Medicaid.

Số bệnh mãn tính:

Relative Risk (RR): 0.4763949 (95% CI: 0.4230582 - 0.5364558) Odds Ratio (OR): 0.3768339 (95% CI: 0.3198127 - 0.4444088)

RR và OR đều nhỏ hơn 1 và khoảng tin cậy 95% không chứa 1. p-value rất nhỏ (< 0.001), cho thấy sự khác biệt có ý nghĩa thống kê cao. Người có ít bệnh mãn tính (≤2) có nguy cơ nhập viện thấp hơn 52.4% (1 - RR) và có khả năng nhập viện thấp hơn 62.3% (1 - OR) so với người có nhiều bệnh mãn tính (>2).

Kết luận:

Bảo hiểm tư nhân không có tác động đáng kể đến khả năng nhập viện. Có Medicaid làm tăng đáng kể khả năng nhập viện. Điều này có thể do những người có Medicaid thường có thu nhập thấp hơn và có thể có nhiều vấn đề sức khỏe hơn. Số bệnh mãn tính có tác động mạnh nhất: những người có nhiều bệnh mãn tính có khả năng nhập viện cao hơn đáng kể.