Untitled

library(readxl)
library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(glmnet)

## Loading required package: Matrix

## Loaded glmnet 4.1-8

# Đọc dữ liệu
TOCAhq <- read_excel("C:/dl/data/R-Rstudio/TOCAhq.xlsx", 
    col_types = c("text", "numeric", "numeric", 
        "text", "text", "numeric", "numeric", 
        "text", "text", "numeric", "numeric", 
        "numeric", "text", "numeric", "numeric", 
        "numeric", "numeric", "text", "numeric", 
        "text", "text", "numeric", "numeric", "text", "numeric", "numeric", "text",  
        "numeric", "numeric", "numeric", 
        "text", "text", "text","text", "text", "text", "text", 
        "text", "text", "numeric", "numeric"))

# Hàm để thay thế giá trị thiếu bằng số trung bình
replace_na_numeric <- function(df, vars) {
  for (var in vars) {
    if (var %in% colnames(df)) {
      mean_value <- mean(df[[var]], na.rm = TRUE)
      df[[var]][is.na(df[[var]])] <- mean_value
    }
  }
  return(df)
}

# Hàm để thay thế giá trị thiếu bằng mode
replace_na_factor <- function(df, vars) {
  for (var in vars) {
    if (var %in% colnames(df)) {
      mode_value <- names(sort(table(df[[var]]), decreasing = TRUE))[1]
      df[[var]][is.na(df[[var]])] <- mode_value
    }
  }
  return(df)
}

# Thay thế giá trị thiếu cho các biến định lượng
TOCAhq <- replace_na_numeric(TOCAhq, c("cTnI-hs2", "dapt", "HDL-C"))

# Thay thế giá trị thiếu cho các biến định tính
TOCAhq <- replace_na_factor(TOCAhq, c("THBH", "RENTROP"))


# Tạo ma trận thiết kế
xvars <- model.matrix(TOCA ~ LAD + LCx + THBH , data = TOCAhq)[,-1]  # Loại bỏ cột intercept
yvar <- TOCAhq$TOCA

# Hồi quy Ridge và Lasso
ridge <- cv.glmnet(xvars, yvar, alpha = 0)
plot(ridge)

coef(ridge)

## 4 x 1 sparse Matrix of class "dgCMatrix"
##                       s1
## (Intercept)  0.240930415
## LAD1        -0.007886670
## LCx1         0.005657099
## THBH1        0.012214880

lasso <- cv.glmnet(xvars, yvar, alpha = 1)
plot(lasso)

coef(lasso)

## 4 x 1 sparse Matrix of class "dgCMatrix"
##                      s1
## (Intercept)  0.24969249
## LAD1        -0.04466726
## LCx1         .         
## THBH1        0.09036207

# Kết quả LASSO cho thấy có các biến như Culprit là nhánh LAD và THBH tiên lượng tốt mà không có hiện tượng đa cộng tuyến ## Đưa các biến  LAD và THBH vào chạy đa biến
# Tạo biến giả từ biến phân loại Culprit
xvars_full <- model.matrix(~ Culprit - 1, data = TOCAhq)

# Chọn biến Culprit2 và 3 từ ma trận thiết kế
xvars_selected <- xvars_full[, "Culprit2", drop = FALSE]
xvars_selected1 <- xvars_full[, "Culprit3", drop = FALSE]
temp1 <- glm(TOCA ~ LAD + THBH, family="binomial", data=TOCAhq)
search1 <- step(temp1)

## Start:  AIC=152.13
## TOCA ~ LAD + THBH
## 
##        Df Deviance    AIC
## <none>      146.13 152.13
## - THBH  1   158.15 162.15
## - LAD   1   159.55 163.55

summary(search1)

## 
## Call:
## glm(formula = TOCA ~ LAD + THBH, family = "binomial", data = TOCAhq)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -0.8785     0.2621  -3.351 0.000805 ***
## LAD1         -1.5581     0.4566  -3.413 0.000643 ***
## THBH1         1.7673     0.5179   3.413 0.000643 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 174.33  on 157  degrees of freedom
## Residual deviance: 146.13  on 155  degrees of freedom
## AIC: 152.13
## 
## Number of Fisher Scoring iterations: 5

# Sử dụng hàm step để chọn biến
search <- step(temp1)

## Start:  AIC=152.13
## TOCA ~ LAD + THBH
## 
##        Df Deviance    AIC
## <none>      146.13 152.13
## - THBH  1   158.15 162.15
## - LAD   1   159.55 163.55

# Tính OR và khoảng tin cậy 95%
OR <- exp(coef(temp1))
conf_int <- exp(confint(temp1))

## Waiting for profiling to be done...

results <- cbind(OR, conf_int)
colnames(results) <- c("OR", "2.5%", "97.5%")
results

##                    OR       2.5%      97.5%
## (Intercept) 0.4154118 0.24368619  0.6846304
## LAD1        0.2105268 0.08109568  0.4952271
## THBH1       5.8547521 2.15648897 16.7536104