library(readxl)
addTaskCallback(function(...) {set.seed(123);TRUE})

##real rawdata
local<-read_xlsx("C:/Users/bkiook/Dropbox/locali/locali0607.xlsx")
local
local<-as.data.frame(local)
library(randomForest)
library(caret)
head(local)

local$emphysema<-as.factor(local$emphysema)
local$solid<-as.factor(local$solid)
local$hemorrhage<-as.factor(local$hemorrhage)
local$dislog<-as.factor(local$dislog)
local$pneumo<-as.factor(local$pneumo)
local$penetration<-as.factor(local$penetration)
local$sex<-as.factor(local$sex)
local$multiple<-as.factor(local$multiple)

head(local)
local<-na.omit(local)
data<-local



# 만휘트니 검정 수행 (wilcox.test)
wilcox_test_result <- wilcox.test(gap  ~ dislog, data = data)
print(wilcox_test_result)
require(moonBook)
require(ztable)
#install.packages("devtools")
require(devtools)
#devtools::install_github("cardiomoon/moonBook") 
#install.packages("devtools")

#devtools::install_github("cardiomoon/ztable")
#devtools::install_github("johanez/rfTools")


require(ztable)
require(rfTools)
table1<-mytable(dislog~., data=data, method=1, show.total=TRUE)
table1

table2<-mytable(dislog~., data=data, method=2, show.total=TRUE)
table2
ztable(table1)




install.packages("DMwR")
library(DMwR)
library(dplyr)
set.seed(123)

balanced_data <- SMOTE(dislog ~ ., data = data, perc.over = 600, perc.under = 150)

balanced_data
table(balanced_data$dislog)
balance<-mytable(dislog~., data=balanced_data, show.total=TRUE)
mycsv(balance, "c:/Users/bkiook/Dropbox/locali/balance.csv", row.names=FALSE)
balance
?mycsv


table(data$dislog)
table(balanced_data$dislog)

balanced_data<-na.omit(balanced_data)
length(balanced_data$sex)

balanced_data

# balanced_data 데이터프레임에서 dislog 변수를 제외한 나머지 변수 선택
x_data <- balanced_data[, !names(balanced_data) %in% c("dislog")]

# 결과 확인
head(x_data)
###랜덤포레스트
#rf_model <- randomForest(dislog ~ ., data=data_rose, importance=TRUE)

rf_model
summary(rf_model)


confusion_matrix <- table(rf_model$predicted, balanced_data$dislog)
print(confusion_matrix)


conf_matrix <- confusionMatrix(rf_model$predicted, balanced_data$dislog, positive="0")
print(conf_matrix)


accuracy <- conf_matrix$overall['Accuracy']
accuracy_CI <- conf_matrix$overall['AccuracyLower'], conf_matrix$overall['AccuracyUpper']
print(paste("Accuracy:", accuracy))
print(paste("95% Confidence Interval for Accuracy:", accuracy_CI[1], "-", accuracy_CI[2]))


# 정확도 계산
accuracy <- sum(diag(confusion_matrix)) / sum(confusion_matrix)
print(paste("Accuracy:", accuracy))
confusionMatrix(rf_model$predicted, balanced_data$dislog)


#지표들

# Confusion Matrix 직접 계산
TP <- 130  # True Positives
FN <- 5    # False Negatives
FP <- 9    # False Positives
TN <- 96   # True Negatives

# 주요 지표 계산
sensitivity <- TP / (TP + FN)
specificity <- TN / (TN + FP)
ppv <- TP / (TP + FP)  # Positive Predictive Value
npv <- TN / (TN + FN)  # Negative Predictive Value

print(paste("Sensitivity:", sensitivity))
print(paste("Specificity:", specificity))
print(paste("PPV:", ppv))
print(paste("NPV:", npv))


# 신뢰구간 계산
alpha <- 0.05  # 95% 신뢰수준

# 민감도 신뢰구간
sens_CI <- binom.test(TP, TP + FN, conf.level = 1 - alpha)$conf.int
print(paste("Sensitivity CI:", round(sens_CI[1], 3), "-", round(sens_CI[2], 3)))

# 특이도 신뢰구간
spec_CI <- binom.test(TN, TN + FP, conf.level = 1 - alpha)$conf.int
print(paste("Specificity CI:", round(spec_CI[1], 3), "-", round(spec_CI[2], 3)))

# 양성예측도 신뢰구간
ppv_CI <- binom.test(TP, TP + FP, conf.level = 1 - alpha)$conf.int
print(paste("PPV CI:", round(ppv_CI[1], 3), "-", round(ppv_CI[2], 3)))

# 음성예측도 신뢰구간
npv_CI <- binom.test(TN, TN + FN, conf.level = 1 - alpha)$conf.int
print(paste("NPV CI:", round(npv_CI[1], 3), "-", round(npv_CI[2], 3)))



# Confusion Matrix
library(caret)



# Confusion Matrix 생성
conf_matrix <- confusionMatrix(rf_model$predicted, balanced_data$dislog)
print(conf_matrix)
# 혼동 행렬 값 추출
TP <- conf_matrix$table[2, 2]  # True Positive
FP <- conf_matrix$table[1, 2]  # False Positive
FN <- conf_matrix$table[2, 1]  # False Negative

# Precision과 Recall 계산
precision <- TP / (TP + FP)  # 양성예측도
recall <- TP / (TP + FN)     # 민감도(재현율)

# F1 Score 계산
f1_score <- 2 * (precision * recall) / (precision + recall)
print(paste("F1 Score:", round(f1_score, 3)))







# 변수 중요도 추출
importance <- importance(rf_model)
varImportance <- data.frame(Variable = rownames(importance), Importance = importance[, 1])

# 변수 중요도 정렬
varImportance <- varImportance[order(varImportance$Importance, decreasing = TRUE), ]




# 변수 중요도 정렬
varImportance <- varImportance[order(varImportance$Importance, decreasing = TRUE), ]
varImportance
plot(varImportance)
barplot(varImportance$Importance, names.arg = varImportance$Variables, las = 2, col = 'lightblue', main = "Variable Importance")
varImpPlot(rf_model)

#save(rf_model, file="c:/Users/bkiook/Dropbox/locali/random_forest_model.RData")

a<-load("c:/Users/bkiook/Dropbox/locali/random_forest_model.RData")
a
rf_model
summary(rf_model)
##행렬
predict_balanced<-predict(rf_model, newdata=balanced_data)
predict_balanced
balanced_data$dislog

# 혼동 행렬 생성
confusion_matrix <- table(predict_balanced, balanced_data$dislog)
print(confusion_matrix)

# 정확도 계산
accuracy <- sum(diag(confusion_matrix)) / sum(confusion_matrix)
print(paste("Accuracy:", accuracy))
confusionMatrix(predict_balanced, balanced_data$dislog)










#복제하지 않은 데이터에 적용
predictions <- predict(rf_model, newdata=data)
print(predictions)

# 혼동 행렬 생성
confusion_matrix <- table(predictions, data$dislog)
print(confusion_matrix)

# 정확도 계산
accuracy <- sum(diag(confusion_matrix)) / sum(confusion_matrix)
print(paste("Accuracy:", accuracy))
conf<-confusionMatrix(predictions, data$dislog, positive="1")
conf$table
conf



# Confusion Matrix 직접 계산
TP <- 101  # True Positives
FN <- 7    # False Negatives
FP <- 0    # False Positives
TN <- 15   # True Negatives

# 주요 지표 계산
sensitivity <- TP / (TP + FN)
specificity <- TN / (TN + FP)
ppv <- TP / (TP + FP)  # Positive Predictive Value
npv <- TN / (TN + FN)  # Negative Predictive Value

print(paste("Sensitivity:", sensitivity))
print(paste("Specificity:", specificity))
print(paste("PPV:", ppv))
print(paste("NPV:", npv))


# 신뢰구간 계산
alpha <- 0.05  # 95% 신뢰수준

# 민감도 신뢰구간
sens_CI <- binom.test(TP, TP + FN, conf.level = 1 - alpha)$conf.int
print(paste("Sensitivity CI:", round(sens_CI[1], 3), "-", round(sens_CI[2], 3)))

# 특이도 신뢰구간
spec_CI <- binom.test(TN, TN + FP, conf.level = 1 - alpha)$conf.int
print(paste("Specificity CI:", round(spec_CI[1], 3), "-", round(spec_CI[2], 3)))

# 양성예측도 신뢰구간
ppv_CI <- binom.test(TP, TP + FP, conf.level = 1 - alpha)$conf.int
print(paste("PPV CI:", round(ppv_CI[1], 3), "-", round(ppv_CI[2], 3)))

# 음성예측도 신뢰구간
npv_CI <- binom.test(TN, TN + FN, conf.level = 1 - alpha)$conf.int
print(paste("NPV CI:", round(npv_CI[1], 3), "-", round(npv_CI[2], 3)))








###parameters


# 트리 깊이 계산 함수 tree depth
calculate_tree_depth <- function(tree) {
  # 각 노드의 깊이 추적
  node_depths <- numeric(nrow(tree))
  node_depths[1] <- 1  # 루트 노드는 깊이 1
  
  for (i in 1:nrow(tree)) {
    if (tree$`left daughter`[i] != 0) {  # 왼쪽 자식 노드
      node_depths[tree$`left daughter`[i]] <- node_depths[i] + 1
    }
    if (tree$`right daughter`[i] != 0) {  # 오른쪽 자식 노드
      node_depths[tree$`right daughter`[i]] <- node_depths[i] + 1
    }
  }
  return(max(node_depths))
}

# 모든 트리의 깊이 계산
tree_depths <- sapply(1:rf_model$ntree, function(k) {
  tree <- getTree(rf_model, k = k, labelVar = TRUE)
  calculate_tree_depth(tree)
})

# 결과 출력
cat("Average tree depth:", mean(tree_depths), "\n")
cat("Maximum tree depth:", max(tree_depths), "\n")
cat("Minimum tree depth:", min(tree_depths), "\n")




##k-fold

library(caret)

# 데이터 준비
data(iris)
set.seed(123)

# trainControl로 cross-validation 설정
control <- trainControl(method = "cv", number = 10)  # 5-fold CV

# 랜덤 포레스트 모델 학습
model <- train(dislog ~ ., 
               data = balanced_data, 
               method = "rf", 
               trControl = control, 
               tuneGrid = expand.grid(mtry = 4)) 
# Resampling 결과 확인

model
resample_results <- model$resample
resample_results
# Accuracy와 Kappa의 95% 신뢰구간 계산
accuracy_ci <- t.test(resample_results$Accuracy)$conf.int
kappa_ci <- t.test(resample_results$Kappa)$conf.int

# 결과 출력
cat("Accuracy 95% Confidence Interval:", accuracy_ci, "\n")
cat("Kappa 95% Confidence Interval:", kappa_ci, "\n")




##xgboost
# 필요한 라이브러리 설치 및 로드
if (!require("xgboost")) install.packages("xgboost")
if (!require("Matrix")) install.packages("Matrix")

library(xgboost)
library(Matrix)


# 범주형 변수 처리: factor로 변환 (One-Hot Encoding을 위해)
balanced_data$emphysema <- as.factor(balanced_data$emphysema)
balanced_data$solid <- as.factor(balanced_data$solid)  # 3개 범주
balanced_data$penetration <- as.factor(balanced_data$penetration)
balanced_data$pneumo <- as.factor(balanced_data$pneumo)
balanced_data$hemorrhage <- as.factor(balanced_data$hemorrhage)
balanced_data$multiple <- as.factor(balanced_data$multiple)
balanced_data <- balanced_data %>%
  mutate(across(c(emphysema, solid, penetration, pneumo, hemorrhage, multiple), 
                ~as.factor(.)))

dummy_data <- dummy.data.frame(balanced_data, names = c("emphysema", "solid", "penetration", "pneumo", "hemorrhage", "multiple"))


# 연속형 변수와 범주형 변수 분리
continuous_vars <- c("age", "size", "totaldepth", "walldepth", "muscledepth", 
                     "nodulepleura", "needle", "protime", "gap")

# One-Hot Encoding을 적용 (범주형 변수에 대해)
X <- model.matrix(~ emphysema + solid + penetration + pneumo + hemorrhage + multiple - 1, 
                  data = balanced_data)
head(X)

# 연속형 변수 처리 (연속형 변수는 그대로 사용)
continuous_vars <- c("age", "size", "totaldepth", "walldepth", "muscledepth", 
                     "nodulepleura", "needle", "protime", "gap", "sex")

# 연속형 변수 추가
X <- cbind(X, balanced_data[, continuous_vars])

# Step 2: 종속 변수 분리
y <- balanced_data$dislog  # 종속 변수

X<-as.matrix(X)
head(X)
X <- apply(X, 2, as.numeric)
head(X)

# Step 3: XGBoost용 데이터 생성
dtrain <- xgb.DMatrix(data = X, label=y)


# Step 4: XGBoost 모델 학습
xgb_model <- xgboost(
  data = dtrain,
  max_depth = 3,       # 트리의 최대 깊이
  eta = 0.01,          # 학습률
  nrounds = 1000,      # 부스팅 반복 횟수
  objective = "binary:logistic",  # 이진 분류
  verbose = 0,
  min_child_weight = 3,
  subsample = 0.8,     # 0.7 ~ 0.9 범위에서 실험 가능
  colsample_bytree = 0.8, # 0.7 ~ 0.9 범위에서 실험 가능
  early_stopping_rounds = 20
)




# Step 5: 변수 중요도 계산
importance_matrix <- xgb.importance(feature_names = colnames(X), model = xgb_model)



# Step 6: 결과 출력 및 시각화
print(importance_matrix)

# 변수 중요도 시각화
xgb.plot.importance(importance_matrix)


# 예측값 계산
pred <- predict(xgb_model, newdata = as.matrix(X))

# 이진 예측값으로 변환 (0.5를 임계값으로)
pred_binary <- ifelse(pred > 0.5, 1, 0)

table (pred_binary, y)





# 정확도, 정밀도, 재현율, F1-score 출력
cat("Accuracy:", conf_matrix$overall["Accuracy"], "\n")
cat("Precision:", conf_matrix$byClass["Pos Pred Value"], "\n")
cat("Recall:", conf_matrix$byClass["Sensitivity"], "\n")
cat("F1-score:", conf_matrix$byClass["F1"], "\n")





##저장
# 패키지 설치
install.packages("officer")
install.packages("rvg")

# 패키지 로드
library(officer)
library(rvg)


figuresupp<-xgb.plot.importance(importance_matrix)