실습

1. 분석을 위해 필요한 라이브러리 및 데이터 불러오기

rm(list = ls())를 이용해 R의 environment를 초기화 할 수 있습니다.

rm(list = ls())

분석에 필요한 패키지들을 불러옵니다.

library(dplyr)        # 데이터 처리

## 
## 다음의 패키지를 부착합니다: 'dplyr'

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

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

library(ggcorrplot)   # 상관관계 시각화

## Warning: 패키지 'ggcorrplot'는 R 버전 4.3.3에서 작성되었습니다

## 필요한 패키지를 로딩중입니다: ggplot2

library(ggplot2)      # 시각화
library(rpart)        # 결정나무
library(rpart.plot)   # 결정나무 시각화
library(randomForest) # 배깅 및 랜덤포레스트

## Warning: 패키지 'randomForest'는 R 버전 4.3.3에서 작성되었습니다

## randomForest 4.7-1.2

## Type rfNews() to see new features/changes/bug fixes.

## 
## 다음의 패키지를 부착합니다: 'randomForest'

## The following object is masked from 'package:ggplot2':
## 
##     margin

## The following object is masked from 'package:dplyr':
## 
##     combine

library(neuralnet)    # 인공신경망

## Warning: 패키지 'neuralnet'는 R 버전 4.3.3에서 작성되었습니다

## 
## 다음의 패키지를 부착합니다: 'neuralnet'

## The following object is masked from 'package:dplyr':
## 
##     compute

library(caret)        # 데이터 분리 및 모형 평가

## Warning: 패키지 'caret'는 R 버전 4.3.3에서 작성되었습니다

## 필요한 패키지를 로딩중입니다: lattice

library(tidyr)        # 데이터 전처리
library(pROC)         # ROC 계산

## Warning: 패키지 'pROC'는 R 버전 4.3.3에서 작성되었습니다

## Type 'citation("pROC")' for a citation.

## 
## 다음의 패키지를 부착합니다: 'pROC'

## The following objects are masked from 'package:stats':
## 
##     cov, smooth, var

KNHANES(국민건강영양조사) 2019년도 데이터를 R로 불러옵니다.
실습 편의상 미리 준비된 HN19_lec3.csv 파일을 불러 옵니다.
head()함수를 통해 불러온 데이터의 앞부분을 확인 가능하며, n = 5와 같이 표시할 행의 수를 명시할 수 있습니다.

hn19 <- read.csv("HN19_lec3.csv")
head(hn19, n = 5) # 데이터의 첫 5개 행 확인

##   age   HE_BMI HE_HP HE_DM_HbA1c
## 1  61 25.98739     3           3
## 2  28 16.90094     1           1
## 3  53 19.78183     2           1
## 4  50 26.63165     1           2
## 5  16       NA    NA          NA

2. 변수 생성 및 결측 처리

주어진 데이터로부터 BMI(체질량지수), 고혈압 여부, 당뇨 여부를 나타내는 범주형 변수를 생성합니다.
본래 고혈압 여부는 HE_HP에 1 = 정상, 2 = 고혈압 전단계, 3 = 고혈압으로 코딩돼 있는데 이를 hp에 0 = 정상, 1 = 고혈압 전단계 또는 고혈압으로 변환하고, 당뇨 여부는 HE_DM_HbA1c에 1 = 정상, 2 = 당뇨병 전단계, 3 = 당뇨병으로 코딩돼 있는데 이를 diabetes에 0 = 정상, 1 = 당뇨병 전단계 또는 고혈압으로 코딩합니다.
결측치(NA)가 포함된 행을 제거합니다.
분석에 필요한 변수만 선택합니다.
head()와 동일한 방법으로, tail()을 이용해 데이터의 뒷부분을 확인할 수 있으며, n = 3과 같이 표시할 행의 수를 명시할 수 있습니다.

hn19 <- hn19 %>%
  mutate(
    bmi = HE_BMI,
    hp = as.numeric(HE_HP %in% c(2, 3)),
    diabetes = as.numeric(HE_DM_HbA1c %in% c(2, 3))
  ) %>%
  filter(!is.na(hp), !is.na(bmi), !is.na(diabetes)) %>%
  select(age, bmi, hp, diabetes)
tail(hn19, n = 3) # 데이터의 마지막 5개 행 확인

##      age      bmi hp diabetes
## 7668  43 15.37515  0        0
## 7669  16 14.84497  0        0
## 7670   8 14.52800  0        0

3. 상관관계 분석 및 시각화

변수 간 상관관계를 확인하면 데이터 구조와 변수 간 관계를 이해하는 데 도움이 됩니다.
상관계수(correlation coefficient)는 두 변수 간 선형 관계를 -1 과 1 사이의 값으로 나타냅니다.
R에서 상관관계는 수치형 설명변수끼리만 구할 수 있습니다.
우선 반응변수 y와 수치형 설명변수 age, bmi, hp, diabetes의 상관계수를 구하고 이를 시각화합니다. hp, diabetes는 실질적으로는 범주형이긴 하지만 as.factor() 변환을 하지 않은 상태이므로 수치형으로 인식돼 상관계수를 구할 수 있습니다.
상관계수를 구할 변수들인 y, age, bmi, hp, diabetes를 select()함수를 이용해서 뽑아내고, cor() 함수를 이용해 상관계수 행렬을 만듭니다.

cor_mat <- hn19 %>%
  select(age, bmi, hp, diabetes) %>%
  cor()

상관계수 행렬을 시각화하기 위해 ggcorrplot패키지의 ggcorrplot() 함수를 이용합니다.
lab = TRUE: 각 셀에 상관계수 숫자를 표시
lab_size = 5: 상관계수 숫자의 크기 지정
tl.cex = 10: 변수명 글자 크기 지정
colors: 상관계수 크기에 따른 색상 팔레트 설정

ggcorrplot(
  cor_mat,
  lab = TRUE,
  lab_size = 5,
  tl.cex = 10,
  colors = c("#6D9EC1", "white", "#E46726")
)

4. BMI와 당뇨병 관계 시각화

변수 간 관계를 시각화하면 데이터의 패턴을 쉽게 확인할 수 있습니다.
ggplot2를 이용해 반응변수와 설명변수에 대한 산점도를 그립니다.
설명변수 중 BMI를 선택해 x축으로 하고, 반응변수인 diabetes를 y축으로 하는 산점도를 그립니다.
반응변수인 diabetes가 0 또는 1이므로 그냥 geom_point를 사용하면 y=0 또는 y=1에 겹쳐지게 됩니다. 따라서 geom_jitter를 사용해 점들이 y=0 또는 y=1 주변으로 흩뿌려지게 합니다.

ggplot(hn19, aes(x = bmi, y = diabetes)) +
  geom_jitter(
    color = "#6D9EC1",
    height = 0.05,
    width = 0,
    alpha = 0.4
  ) +
  labs(
    title = "BMI vs 당뇨병",
    x = "BMI",
    y = "당뇨병"
  ) +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5))

5. 로지스틱 회귀 (2변수 모형) 및 예측 확률 시각화

로지스틱 회귀는 종속변수가 이항형(0/1)일 때 확률을 예측하는데 사용됩니다.
주어진 데이터에서 우선 설명변수를 2개 사용하여 age, bmi를 설명변수로, diabetes를 반응변수로 하는 로지스틱 회귀모형을 적합합니다.
R에서 로지스틱 회귀모형을 적합하는 경우 glm() 함수를 이용합니다. 문법은 lm과 대동소이합니다. 함수 호출 시 전달해야 하는 인자는 다음과 같습니다.
- formula: 적합할 회귀모형의 형태를 나타냅니다. 반응변수 ~ 설명변수1 + \(\cdots\) + 설명변수p의 형태와 같이 반응변수와 설명변수는 ~로 연결하고, 설명변수끼리는 +로 연결합니다. 이 예제의 경우 반응변수는 diabetes, 설명변수는 age, bmi이므로 다음과 같이 인자를 전달할 수 있습니다. formula = diabetes ~ age + bmi
- data: 사용할 데이터를 명시합니다. 앞의 formula에서 명시된 변수들이 인자로 전달되는 데이터에 모두 포함되어 있어야 합니다. 우리의 경우 데이터를 hn19라는 변수명으로 저장하였으므로 data = hn19와 같이 인자를 전달합니다.
- family: 로지스틱 회귀모형의 경우에는 "binomial"을 사용합니다.
glm() 함수를 이용해 선형회귀모형 적합을 완료한 후, summary()함수를 이용하여 회귀모형 적합 결과를 확인할 수 있습니다.
summary()의 결과를 통해 개별 회귀계수 추정치와 이에 대한 t-검정 경과는 물론 전체 로지스틱 회귀모형의 유의성에 대한 F-검정 결과, R-squared와 adjusted R-squared를 모두 확인할 수 있습니다.

model_logistic_2d <- glm(
  formula = diabetes ~ age + bmi,
  data = hn19,
  family = "binomial"
)
summary(model_logistic_2d)

## 
## Call:
## glm(formula = diabetes ~ age + bmi, family = "binomial", data = hn19)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -6.994652   0.217239  -32.20   <2e-16 ***
## age          0.062048   0.001672   37.11   <2e-16 ***
## bmi          0.160854   0.008059   19.96   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 10488.7  on 7669  degrees of freedom
## Residual deviance:  7472.3  on 7667  degrees of freedom
## AIC: 7478.3
## 
## Number of Fisher Scoring iterations: 5

시각화를 위해, 격자 데이터를 grid에 생성합니다.
격자의 각 지점별 당뇨 확률을 계산해 pred_logistic_2d에 저장합니다.
ggplot()을 이용해 색상으로 시각화합니다.

age_seq <- seq(min(hn19$age), max(hn19$age), length.out = 100)
bmi_seq <- seq(min(hn19$bmi), max(hn19$bmi), length.out = 100)
grid <- expand.grid(age = age_seq, bmi = bmi_seq)

pred_logistic_2d <- predict(
  model_logistic_2d,
  newdata = grid,
  type = "response"
)

ggplot() +
  geom_tile(
    aes(
      x = grid[["age"]],
      y = grid[["bmi"]],
      fill = pred_logistic_2d
    ),
    alpha = 0.4
  ) +
  scale_fill_gradient(low = "lightblue", high = "red") +
  labs(
    title = "로지스틱 회귀의 당뇨병 예측 확률",
    x = "나이 (Age)",
    y = "BMI",
    fill = "P(당뇨병=1)"
  ) +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5))

6. 로지스틱 회귀 (다변량 모형) 및 잔차 진단

이번에는 설명변수를 3개 사용하여 age, bmi, hp를 설명변수로, diabetes를 반응변수로 하는 로지스틱 회귀모형을 적합합니다.
summary() 함수를 통해 계수 및 유의성, 모델 적합도를 확인합니다.

model_logistic <- glm(
  formula = diabetes ~ age + bmi + hp,
  data = hn19,
  family = "binomial"
)
summary(model_logistic)

## 
## Call:
## glm(formula = diabetes ~ age + bmi + hp, family = "binomial", 
##     data = hn19)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -6.741318   0.220841 -30.526  < 2e-16 ***
## age          0.057853   0.001831  31.589  < 2e-16 ***
## bmi          0.150998   0.008249  18.304  < 2e-16 ***
## hp           0.341577   0.063894   5.346 8.99e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 10489  on 7669  degrees of freedom
## Residual deviance:  7444  on 7666  degrees of freedom
## AIC: 7452
## 
## Number of Fisher Scoring iterations: 5

적합한 회귀모형 model_logistic에 함수 plot()을 적용함으로써 잔차 및 진단 그래프를 확인할 수 있습니다.

plot(model_logistic)

7. 결정 나무 (2변수와 다변량 모델 + 가지치기)

결정나무(Decision Tree)는 변수 공간을 여러 영역으로 나누어 각 영역마다 같은 예측값을 갖도록 하는 모델입니다.
해석이 쉽고 비선형 관계를 잡아낼 수 있지만, 깊게 성장하면 과적합(overfitting)이 발생할 수 있습니다.
이를 방지하기 위해 가지치기(pruning) 를 수행하여 최적의 복잡도를 선택합니다
결정 나무의 경우에는 분류로 사용하려면 반응변수에 factor 형태의 변수를 넣어야 합니다. as.factor()를 사용하여 당뇨병 여부를 factor로 변환하고 이를 diabetes_factor에 저장합니다.

hn19 <- hn19 %>% mutate(diabetes_factor = as.factor(diabetes))

주어진 데이터에서 우선 설명변수를 2개 사용하여 age, bmi를 설명변수로, diabetes를 반응변수로 하는 분류 결정나무를 적합합니다.
R에서 결정나무를 적합하는 경우 rpart 패키지의 rpart()함수를 이용합니다. 문법은lm`과 대동소이합니다. 함수 호출 시 전달해야 하는 인자는 다음과 같습니다.
- formula: 적합할 회귀모형의 형태를 나타냅니다. 반응변수 ~ 설명변수1 + \(\cdots\) + 설명변수p의 형태와 같이 반응변수와 설명변수는 ~로 연결하고, 설명변수끼리는 +로 연결합니다. 이 예제의 경우 반응변수는 diabetes, 설명변수는 age, bmi이므로 다음과 같이 인자를 전달할 수 있습니다. formula = diabetes ~ age + bmi
- data: 사용할 데이터를 명시합니다. 앞의 formula에서 명시된 변수들이 인자로 전달되는 데이터에 모두 포함되어 있어야 합니다. 우리의 경우 데이터를 hn19라는 변수명으로 저장하였으므로 data = hn19와 같이 인자를 전달합니다.
- method: 분류 결정나무를 사용하는 경우 "class"을 사용합니다.
rpart.plot() 함수를 사용해 적합한 결정나무를 시각화합니다.

model_tree_2d <- rpart(
  diabetes_factor ~ age + bmi,
  data = hn19,
  method = "class"
)
rpart.plot(model_tree_2d, type = 2, extra = 106, cex = 1.2)

앞에서와 마찬가지로 시각화를 합니다. 격자 데이터는 앞에서 생성한 grid를 이용합니다.
격자의 각 지점별 당뇨 확률을 계산해 pred_tree_2d에 저장합니다.
ggplot()을 이용해 색상으로 시각화합니다.
결정 나무의 경우에는 경계선이 항상 축과 평행하게 나옴을 알 수 있습니다.

pred_tree_2d <-
  predict(model_tree_2d, newdata = grid, type = "prob")[, 2]

ggplot() +
  geom_tile(
    aes(
      x = grid[["age"]],
      y = grid[["bmi"]],
      fill = pred_tree_2d
    ),
    alpha = 0.4
  ) +
  scale_fill_gradient(low = "lightblue", high = "red") +
  labs(
    title = "결정 나무의 영역별 당뇨병 예측 확률",
    x = "나이(Age)",
    y = "BMI",
    fill = "P(당뇨병=1)"
  ) +
  theme_minimal() + 
  theme(plot.title = element_text(hjust = 0.5))

이번에는 설명변수를 3개 사용하여 age, bmi, hp를 설명변수로, diabetes를 반응변수로 하는 분류 결정나무를 적합하고 가지치기(pruning)을 합니다.
control은 결정나무를 적합하는 데 추가적으로 사용하는 옵션을 지정해 줍니다. 나중에 가지치기를 하기 위해, 우선 cp = 0.001로 설정하여 우선 결정나무를 충분히 복잡하게 만듭니다.

model_tree <- rpart(
  diabetes_factor ~ age + bmi + hp,
  data = hn19,
  method = "class",
  control = rpart.control(cp = 0.001)
)

printcp() 함수로 각 단계에서의 cp 값을 출력합니다.
best_cp에는 교차검증 오차가 최소가 되는 cp 값이 저장되게 됩니다.
prune() 함수로 가지치기(pruning)을 해서 최적의 결정나무를 model_ptree에 저장합니다.
rpart.plot() 함수를 사용해 가지치기한 결정나무를 시각화합니다.

printcp(model_tree)

## 
## Classification tree:
## rpart(formula = diabetes_factor ~ age + bmi + hp, data = hn19, 
##     method = "class", control = rpart.control(cp = 0.001))
## 
## Variables actually used in tree construction:
## [1] age bmi hp 
## 
## Root node error: 3310/7670 = 0.43155
## 
## n= 7670 
## 
##          CP nsplit rel error  xerror     xstd
## 1 0.3927492      0   1.00000 1.00000 0.013105
## 2 0.0253776      1   0.60725 0.60906 0.011647
## 3 0.0021148      3   0.55650 0.56979 0.011394
## 4 0.0015710      8   0.54411 0.56949 0.011392
## 5 0.0013595     15   0.53263 0.56284 0.011346
## 6 0.0012085     17   0.52991 0.55559 0.011296
## 7 0.0010070     32   0.50906 0.55529 0.011294
## 8 0.0010000     35   0.50604 0.55498 0.011292

best_cp <- model_tree$cptable[which.min(model_tree$cptable[, "xerror"]), "CP"]
model_ptree <- prune(model_tree, cp = best_cp)
rpart.plot(model_ptree, type = 2, extra = 106, cex = 1)

8. 앙상블 모형(배깅, 랜덤 포레스트)

배깅(Bagging) 은 여러 부트스트랩 샘플에 대해 모델(결정나무)을 학습하고 예측을 평균/다수결로 결합합니다.
랜덤 포레스트(Random Forest) 는 배깅에 변수 무작위 선택을 추가하여 상관성을 줄이고 성능을 향상시킵니다.
배깅과 랜덤 포레스트는 모두 randomForest 패키지의 randomForest() 함수로 적합할 수 있습니다. 이는 랜덤 포레스트에서 변수의 개수를 전부 다 쓰면 그게 배깅이 되기 때문입니다.
앞에서와 마찬가지로 설명변수를 3개 사용하여 age, bmi, hp를 설명변수로, diabetes를 반응변수로 하는 앙상블 모형을 생각합니다.
우선 배깅 모형을 적합합니다. 이는 randomForest() 함수에서 사용하는 설명변수의 개수를 전체 설명변수의 개수인 3으로 지정하면 됩니다 (mtry=3).
print() 함수로 적합한 배깅 모형의 요약 정보를 출력합니다.
plot() 함수로 적합한 트리 수에 따른 오차 변화를 확인합니다.

model_bag <- randomForest(
  diabetes_factor ~ age + bmi + hp,
  data = hn19,
  mtry = 3,
  ntree = 100,
  importance = TRUE
)
print(model_bag)

## 
## Call:
##  randomForest(formula = diabetes_factor ~ age + bmi + hp, data = hn19,      mtry = 3, ntree = 100, importance = TRUE) 
##                Type of random forest: classification
##                      Number of trees: 100
## No. of variables tried at each split: 3
## 
##         OOB estimate of  error rate: 28.45%
## Confusion matrix:
##      0    1 class.error
## 0 3220 1140   0.2614679
## 1 1042 2268   0.3148036

plot(model_bag)

다음으로 랜덤포레스트 모형을 적합합니다. 이는 randomForest() 함수에서 사용하는 설명변수의 개수를 전체 설명변수의 개수보다 적게 지정하면 됩니다. (mtry=2).
print() 함수로 적합한 배깅 모형의 요약 정보를 출력합니다.
plot() 함수로 적합한 트리 수에 따른 오차 변화를 확인합니다.

model_rf <- randomForest(
  diabetes_factor ~ age + bmi + hp,
  data = hn19,
  mtry = 2,
  ntree = 100,
  importance = TRUE
)
print(model_rf)

## 
## Call:
##  randomForest(formula = diabetes_factor ~ age + bmi + hp, data = hn19,      mtry = 2, ntree = 100, importance = TRUE) 
##                Type of random forest: classification
##                      Number of trees: 100
## No. of variables tried at each split: 2
## 
##         OOB estimate of  error rate: 27.46%
## Confusion matrix:
##      0    1 class.error
## 0 3222 1138   0.2610092
## 1  968 2342   0.2924471

plot(model_rf)

9. 신경망 모형 (Neural Network)

신경망(Neural Network)은 입력층, 은닉층, 출력층으로 구성되며, 비선형 관계 학습에 강점이 있습니다.
충분한 데이터와 계산 자원이 필요하지만, 복잡한 패턴을 학습할 수 있습니다.
앞에서와 마찬가지로 설명변수를 3개 사용하여 age, bmi, hp를 설명변수로, diabetes를 반응변수로 하는, 은닉층(hidden layer)이 하나인 신경망 모형을 생각합니다.
신경망 모형을 적합하는 패키지는 여러 가지 있지만, 여기서는 가장 간단하게 할 수 있는 neuralnet 패키지의 neuralnet() 함수를 소개합니다. 다만 이 패키지는 느린 편이어서, 더 복잡한 자료에서는 keras 패키지 내지는 torch 패키지를 사용하면 좀 더 빠르게 실행할 수 있습니다.
속도가 느린 관계로 hn19[1:100, ] 명령어로 처음 100개의 행만 골라냅니다.
hidden은 은닉층(hidden layer)의 뉴런 수입니다. 아래에서는 hidden = 2 를 사용했습니다.
출력층은 로지스틱 회귀처럼 0~1 범위의 확률을 출력하도록 linear.output = FALSE 사용합니다.
plot() 함수로 신경망 구조를 시각화합니다.

model_nn <- neuralnet(
  diabetes ~ age + bmi + hp,
  data = hn19[1:100, ],
  hidden = 2,
  linear.output = FALSE
)
plot(model_nn)

10. 모델 평가 (데이터 분리, Threshold 성능 분석, ROC Curve)

혼동행렬, 민감도, 특이도(Slide 60~64)로 모델 성능을 평가합니다.
threshold 변화에 따른 민감도·특이도 변화(ROC Curve)를 시각화하고, AUC를 계산해 분류기 성능을 종합적으로 판단합니다.
학습용/테스트용 데이터 분리를 통해 과대적합을 방지합니다.
우선 학습용/테스트용 데이터(training data/test data)를 10% 대 90%로 분리합니다. caret 패키지의 createDataPartition() 함수를 사용합니다.

set.seed(2025)
id_train <- createDataPartition(hn19[["diabetes"]], p = 0.1, list = FALSE)
hn19_train <- hn19[id_train, ]
hn19_test <- hn19[-id_train, ]

다음으로,threshold 별 평가 지표(accuracy, sensitivity, specificity)를 한 번에 계산하는 코드를 함수로 작성합니다.
thresholds 에는 threshold 의 값이 벡터로 들어옵니다.
actual은 실제 반응변수의 값, predicted는 예측된 반응변수의 값이 들어옵니다. 단, 여기서 predicted는 확률로 값이 들어가야 합니다.

calc_metrics <- function(thresholds, actual, predicted) {

  # threshold 하나에 대해 평가 지표를 계산
  metrics_one <- function(threshold, actual, predicted) {
    # 예측 확률이 threshold보다 크면 1(양성), 작으면 0(음성)
    pred_class <- ifelse(predicted > threshold, 1, 0)

    # 혼동 행렬 요소 계산
    TP <- sum(pred_class == 1 & actual == 1) # 참 양성
    TN <- sum(pred_class == 0 & actual == 0) # 참 음성
    FP <- sum(pred_class == 1 & actual == 0) # 거짓 양성
    FN <- sum(pred_class == 0 & actual == 1) # 거짓 음성

    # 평가 지표 계산
    accuracy <- (TP + TN) / (TP + TN + FP + FN)
    sensitivity <- TP / (TP + FN)   # 재현율
    specificity <- TN / (TN + FP)   # 특이도
    c(accuracy, sensitivity, specificity)
  }
  
  # 모든 threshold에 대해 metrics_one 계산
  results <- data.frame(
    threshold = thresholds,
    t(sapply(thresholds, function(t) metrics_one(t, actual, predicted)))
  )
  colnames(results) <- c("threshold", "Accuracy", "Sensitivity", "Specificity")
  results
}

thresholds 후보를 0부터 1까지 0.01 단위의 벡터로 만듭니다.

thresholds <- seq(0, 1, by = 0.01)
id_05 <- which(thresholds == 0.5)

앞과 마찬가지로 age, bmi, hp를 설명변수로, diabetes를 반응변수로 하는 모형을 생각합니다.
아래 코드에서는 로지스틱 회귀모형을 생각했습니다. 학습 데이터(training data)에서 모형을 학습하고, 테스트 데이터(test data)에서 예측 확률을 산출합니다.

model_logistic <- glm(
  diabetes ~ age + bmi + hp,
  data = hn19_test,
  family = "binomial"
)
pred_logistic <- predict(model_logistic, newdata = hn19_test, type = "response")

앞에서 작성한 calc_metrics() 함수로 threshold 변화 시 accuracy, sensitivity, specificity 를 계산합니다. 특히나, accuracy 에서 threshold = 0.5에 해당하는 값이 우리가 일반적으로 생각하는 정확도입니다.

metric_logistic <- calc_metrics(thresholds, hn19_test[["diabetes"]], pred_logistic)
cat("Logistic Accuracy:", round(metric_logistic[["Accuracy"]][id_05], 3), "\n")

## Logistic Accuracy: 0.756

위에서 계산한 값을 바탕으로 accuracy, sensitivity, specificity 를 시각화합니다.

metric_logistic_long <- metric_logistic %>%
  pivot_longer(
    cols = c("Accuracy", "Sensitivity", "Specificity"),
    names_to = "Metric",
    values_to = "Value"
  )
ggplot(metric_logistic_long, aes(x = threshold, y = Value, color = Metric)) +
  geom_line(size = 1.2) +
  labs(title = "Accuracy / Sensitivity / Specificity", x = "Threshold", y = "Metric Value") +
  theme_minimal() + 
  theme(plot.title = element_text(hjust = 0.5))

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

ROC curve 를 계산 및 시각화하고 AUC 를 계산합니다. pROC 패키지의 roc() 함수를 이용하면 됩니다.

roc_obj <- roc(hn19_test[["diabetes"]], pred_logistic)

## Setting levels: control = 0, case = 1

## Setting direction: controls < cases

plot(roc_obj, main = paste("ROC Curve (AUC =", round(roc_obj$auc, 3), ")"))

연습문제

똑같은 데이터를 사용합니다.

1. 로지스틱 회귀 모형

설명변수를 1개 사용하여 bmi를 설명변수로, diabetes를 반응변수로 하는 로지스틱 회귀모형을 적합합니다.
4단계에서 그린 BMI와 당뇨병 관계 시각화에 더해서 BMI 값 별로 로지스틱 회귀모형의 당뇨병 예측 확률을 곡선으로 그립니다. geom_line(aes(y = pred_logistic))과 같은 형태로 할 수 있습니다.

풀이

우선 아래와 같이 BMI와 당뇨병 관계를 시각화는 코드를 가져옵니다.

p <- ggplot(hn19, aes(x = bmi, y = diabetes)) +
  geom_jitter(
    color = "#6D9EC1",
    height = 0.05,
    width = 0,
    alpha = 0.4
  ) +
  labs(
    title = "BMI vs 당뇨병",
    x = "BMI",
    y = "당뇨병"
  ) +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5))

p

이제 bmi를 설명변수로, diabetes를 반응변수로 하는 로지스틱 회귀모형을 적합합니다.

model_logistic_bmi <- glm(
  formula = diabetes ~ bmi,
  data = hn19,
  family = "binomial"
)

bmi에 해당하는 격자를 grid_bmi로 생성합니다.
격자의 각 지점별 당뇨 확률을 계산해 grid_bmi[["pred_logistic"]]에 저장합니다.
ggplot()의 geom_line()을 이용해 회귀곡선을 시각화합니다. grid_bmi의 bmi에 격자가, pred_logistic에 당뇨 확률이 계산돼 있으므로 aes(x = bmi, y = pred_logistic) 같은 식으로 코드를 작성하면 됩니다.

grid_bmi <- data.frame(
  bmi = seq(min(hn19[["bmi"]]), max(hn19[["bmi"]]), length.out = 100))

grid_bmi[["pred_logistic"]] <- predict(
  model_logistic_bmi,
  newdata = grid_bmi,
  type = "response"
)

p <- p +
  geom_line(data = grid_bmi, aes(x = bmi, y = pred_logistic),
    color = "blue", size = 1.2)
p

2. 결정나무

설명변수를 1개 사용하여 bmi를 설명변수로, diabetes를 반응변수로 하는 분류 결정나무를 적합합니다. 가지치기는 따로 안해도 됩니다.
4단계에서 그린 BMI와 당뇨병 관계 시각화에 더해서 BMI 값 별로 분류 결정나무의 당뇨병 예측 확률을 곡선으로 그립니다. geom_line(aes(y = pred_tree))과 같은 형태로 할 수 있습니다.

풀이

우선 아래와 같이 BMI와 당뇨병 관계를 시각화는 코드를 가져옵니다.

p <- ggplot(hn19, aes(x = bmi, y = diabetes)) +
  geom_jitter(
    color = "#6D9EC1",
    height = 0.05,
    width = 0,
    alpha = 0.4
  ) +
  labs(
    title = "BMI vs 당뇨병",
    x = "BMI",
    y = "당뇨병"
  ) +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5))

p

이제 bmi를 설명변수로, diabetes를 반응변수로 하는 분류 결정나무를 적합합니다.

model_tree_bmi <- rpart(
  diabetes_factor ~ bmi,
  data = hn19,
  method = "class"
)

앞 문제의 grid_bmi를 bmi에 해당하는 격자로 활용합니다.
격자의 각 지점별 당뇨 확률을 계산해 grid_bmi[["pred_tree"]]에 저장합니다.
ggplot()의 geom_line()을 이용해 회귀곡선을 시각화합니다. grid_bmi의 bmi에 격자가, pred_tree에 당뇨 확률이 계산돼 있으므로 aes(x = bmi, y = pred_tree) 같은 식으로 코드를 작성하면 됩니다.

grid_bmi[["pred_tree"]] <- predict(
  model_tree_bmi,
  newdata = grid_bmi,
  type = "prob"
)[, 2]

p <- p +
  geom_line(data = grid_bmi, aes(x = bmi, y = pred_tree),
    color = "blue", size = 1.2)
p

3. 배깅

설명변수를 2개 사용하여 age, bmi를 설명변수로, diabetes를 반응변수로 하는 배깅 모형을 적합합니다.
5단계 마지막에서 했던 것처럼 age 및 BMI 값 별로 배깅 모형의 당뇨병 예측 확률을 색칠해 나타냅니다.

풀이

우선 설명변수를 2개 사용하여 age, bmi를 설명변수로, diabetes를 반응변수로 하는 배깅 모형을 적합합니다.

model_bag_2d <- randomForest(
  diabetes_factor ~ age + bmi,
  data = hn19,
  mtry = 2,
  ntree = 100,
  importance = TRUE
)

그러고서는 5단계 마지막에 했던 시각화를 반복합니다.
우선 격자 데이터를 grid에 생성합니다.
격자의 각 지점별 당뇨 확률을 계산해 pred_bag_2d에 저장합니다.
ggplot()을 이용해 색상으로 시각화합니다.

age_seq <- seq(min(hn19$age), max(hn19$age), length.out = 100)
bmi_seq <- seq(min(hn19$bmi), max(hn19$bmi), length.out = 100)
grid <- expand.grid(age = age_seq, bmi = bmi_seq)

pred_bag_2d <- predict(
  model_bag_2d,
  newdata = grid,
  type = "prob"
)[, 2]

ggplot() +
  geom_tile(
    aes(
      x = grid[["age"]],
      y = grid[["bmi"]],
      fill = pred_bag_2d
    ),
    alpha = 0.4
  ) +
  scale_fill_gradient(low = "lightblue", high = "red") +
  labs(
    title = "배깅의 당뇨병 예측 확률",
    x = "나이 (Age)",
    y = "BMI",
    fill = "P(당뇨병=1)"
  ) +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5))

4. 랜덤 포레스트

설명변수를 2개 사용하여 age, bmi를 설명변수로, diabetes를 반응변수로 하는 랜덤 포레스트 모형을 적합합니다.
5단계 마지막에서 했던 것처럼 age 및 BMI 값 별로 랜덤 포레스트 모형의 당뇨병 예측 확률을 색칠해 나타냅니다.

풀이

우선 설명변수를 2개 사용하여 age, bmi를 설명변수로, diabetes를 반응변수로 하는 랜덤 포레스트 모형을 적합합니다.

model_rf_2d <- randomForest(
  diabetes_factor ~ age + bmi,
  data = hn19,
  mtry = 1,
  ntree = 100,
  importance = TRUE
)

그러고서는 5단계 마지막에 했던 시각화를 반복합니다.
앞 문제에서 생성한 grid를 격자로 활용합니다.
격자의 각 지점별 당뇨 확률을 계산해 pred_bag_2d에 저장합니다.
ggplot()을 이용해 색상으로 시각화합니다.

age_seq <- seq(min(hn19$age), max(hn19$age), length.out = 100)
bmi_seq <- seq(min(hn19$bmi), max(hn19$bmi), length.out = 100)
grid <- expand.grid(age = age_seq, bmi = bmi_seq)

pred_rf_2d <- predict(
  model_rf_2d,
  newdata = grid,
  type = "prob"
)[, 2]

ggplot() +
  geom_tile(
    aes(
      x = grid[["age"]],
      y = grid[["bmi"]],
      fill = pred_rf_2d
    ),
    alpha = 0.4
  ) +
  scale_fill_gradient(low = "lightblue", high = "red") +
  labs(
    title = "랜덤 포레스트의 당뇨병 예측 확률",
    x = "나이 (Age)",
    y = "BMI",
    fill = "P(당뇨병=1)"
  ) +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5))

5. 인공 신경망

설명변수를 2개 사용하여 age, bmi를 설명변수로, diabetes를 반응변수로 하는 인공 신경망을 적합합니다.
5단계 마지막에서 했던 것처럼 age 및 BMI 값 별로 인공 신경망 모형의 당뇨병 예측 확률을 색칠해 나타냅니다.

풀이

우선 설명변수를 2개 사용하여 age, bmi를 설명변수로, diabetes를 반응변수로 하는 인공 신경망을 적합합니다.

model_nn_2d <- neuralnet(
  diabetes ~ age + bmi,
  data = hn19[1:100, ],
  hidden = 2,
  linear.output = FALSE
)
plot(model_nn_2d)

그러고서는 5단계 마지막에 했던 시각화를 반복합니다.
앞 문제에서 생성한 grid를 격자로 활용합니다.
격자의 각 지점별 당뇨 확률을 계산해 pred_bag_2d에 저장합니다.
ggplot()을 이용해 색상으로 시각화합니다.

pred_nn_2d <- predict(
  model_nn_2d,
  newdata = grid,
  type = "prob"
)

ggplot() +
  geom_tile(
    aes(
      x = grid[["age"]],
      y = grid[["bmi"]],
      fill = pred_nn_2d
    ),
    alpha = 0.4
  ) +
  scale_fill_gradient(low = "lightblue", high = "red") +
  labs(
    title = "인공 신경망의 당뇨병 예측 확률",
    x = "나이 (Age)",
    y = "BMI",
    fill = "P(당뇨병=1)"
  ) +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5))

6. 모델 평가

설명변수를 3개 사용하여 age, bmi, hp를 설명변수로, diabetes를 반응변수로 하는 모형을 생각합니다.
위 실습에서는 로지스틱 회귀모형에 대해서 accuracy, sensitivity, specificity 의 곡선을 그리고 roc curve 를 그렸습니다. 이를 결정나무, 배깅, 랜덤 포레스트, 인공신경망에 대해 반복합니다.

풀이

앞에서처럼 학습용/테스트용 데이터(training data/test data)를 10% 대 90%로 분리합니다. caret 패키지의 createDataPartition() 함수를 사용합니다.

set.seed(2025)
id_train <- createDataPartition(hn19[["diabetes"]], p = 0.1, list = FALSE)
hn19_train <- hn19[id_train, ]
hn19_test <- hn19[-id_train, ]

앞에서 작성했던 threshold 별 평가 지표(accuracy, sensitivity, specificity)를 한 번에 계산하는 calc_metrics 함수를 가져옵니다.

calc_metrics <- function(thresholds, actual, predicted) {

  # threshold 하나에 대해 평가 지표를 계산
  metrics_one <- function(threshold, actual, predicted) {
    # 예측 확률이 threshold보다 크면 1(양성), 작으면 0(음성)
    pred_class <- ifelse(predicted > threshold, 1, 0)

    # 혼동 행렬 요소 계산
    TP <- sum(pred_class == 1 & actual == 1) # 참 양성
    TN <- sum(pred_class == 0 & actual == 0) # 참 음성
    FP <- sum(pred_class == 1 & actual == 0) # 거짓 양성
    FN <- sum(pred_class == 0 & actual == 1) # 거짓 음성

    # 평가 지표 계산
    accuracy <- (TP + TN) / (TP + TN + FP + FN)
    sensitivity <- TP / (TP + FN)   # 재현율
    specificity <- TN / (TN + FP)   # 특이도
    c(accuracy, sensitivity, specificity)
  }
  
  # 모든 threshold에 대해 metrics_one 계산
  results <- data.frame(
    threshold = thresholds,
    t(sapply(thresholds, function(t) metrics_one(t, actual, predicted)))
  )
  colnames(results) <- c("threshold", "Accuracy", "Sensitivity", "Specificity")
  results
}

여기에 더해, accuracy, sensitivity, specificity 를 시각화하는 함수, roc curve 를 그리는 함수를 각각 만듭니다.

plot_metrics <- function(thresholds, actual, predicted) {

  metric <- calc_metrics(thresholds, actual, predicted)

  metric_long <- metric %>%
  pivot_longer(
    cols = c("Accuracy", "Sensitivity", "Specificity"),
    names_to = "Metric",
    values_to = "Value"
  )
  
  ggplot(metric_long, aes(x = threshold, y = Value, color = Metric)) +
  geom_line(size = 1.2) +
  labs(title = "Accuracy / Sensitivity / Specificity", x = "Threshold", y = "Metric Value") +
  theme_minimal() + 
  theme(plot.title = element_text(hjust = 0.5))
}

plot_roc <- function(actual, predicted) {

  roc_obj <- roc(actual, predicted)
  plot(roc_obj, main = paste("ROC Curve (AUC =", round(roc_obj[["auc"]], 3), ")"))
}

앞에서처럼 thresholds 후보를 0부터 1까지 0.01 단위의 벡터로 만듭니다.

thresholds <- seq(0, 1, by = 0.01)
id_05 <- which(thresholds == 0.5)

이제 결정나무에 대해 accuracy, sensitivity, specificity 의 곡선을 그리고 roc curve 를 그립니다.
학습 데이터(training data)에서 모형을 학습하고, 테스트 데이터(test data)에서 예측 확률을 산출합니다.

model_tree <- rpart(
  diabetes_factor ~ age + bmi + hp,
  data = hn19,
  method = "class"
)
pred_tree <- predict(model_tree, newdata = hn19_test, type = "prob")[, 2]

앞에서 작성한 plot_metrics() 함수로 threshold 변화 시 accuracy, sensitivity, specificity 를 시각화합니다.

plot_metrics(thresholds, hn19_test[["diabetes"]], pred_tree)

ROC curve 를 계산 및 시각화하고 AUC 를 계산합니다.

plot_roc(hn19_test[["diabetes"]], pred_tree)

## Setting levels: control = 0, case = 1

## Setting direction: controls < cases

이제 배깅에 대해 accuracy, sensitivity, specificity 의 곡선을 그리고 roc curve 를 그립니다.
학습 데이터(training data)에서 모형을 학습하고, 테스트 데이터(test data)에서 예측 확률을 산출합니다.

model_bag <-  randomForest(
  diabetes_factor ~ age + bmi + hp,
  data = hn19_train,
  mtry = 3,
  ntree = 100,
  importance = TRUE
)
pred_bag <- predict(model_bag, newdata = hn19_test, type = "prob")[, 2]

앞에서 작성한 plot_metrics() 함수로 threshold 변화 시 accuracy, sensitivity, specificity 를 시각화합니다.

plot_metrics(thresholds, hn19_test[["diabetes"]], pred_bag)

ROC curve 를 계산 및 시각화하고 AUC 를 계산합니다.

plot_roc(hn19_test[["diabetes"]], pred_bag)

## Setting levels: control = 0, case = 1

## Setting direction: controls < cases

이제 랜덤포레스트에 대해 accuracy, sensitivity, specificity 의 곡선을 그리고 roc curve 를 그립니다.
학습 데이터(training data)에서 모형을 학습하고, 테스트 데이터(test data)에서 예측 확률을 산출합니다.

model_rf <-  randomForest(
  diabetes_factor ~ age + bmi + hp,
  data = hn19_train,
  mtry = 2,
  ntree = 100,
  importance = TRUE
)
pred_rf <- predict(model_rf, newdata = hn19_test, type = "prob")[, 2]

앞에서 작성한 plot_metrics() 함수로 threshold 변화 시 accuracy, sensitivity, specificity 를 시각화합니다.

plot_metrics(thresholds, hn19_test[["diabetes"]], pred_rf)

ROC curve 를 계산 및 시각화하고 AUC 를 계산합니다.

plot_roc(hn19_test[["diabetes"]], pred_rf)

## Setting levels: control = 0, case = 1

## Setting direction: controls < cases

마지막으로 인공신경망에 대해 accuracy, sensitivity, specificity 의 곡선을 그리고 roc curve 를 그립니다.
학습 데이터(training data)에서 모형을 학습하고, 테스트 데이터(test data)에서 예측 확률을 산출합니다.

model_nn <-  neuralnet(
  diabetes ~ age + bmi + hp,
  data = hn19_train,
  hidden = 2,
  linear.output = FALSE
)
pred_nn <- predict(model_nn, newdata = hn19_test, type = "prob")

앞에서 작성한 plot_metrics() 함수로 threshold 변화 시 accuracy, sensitivity, specificity 를 시각화합니다.

plot_metrics(thresholds, hn19_test[["diabetes"]], pred_nn)

ROC curve 를 계산 및 시각화하고 AUC 를 계산합니다.

plot_roc(hn19_test[["diabetes"]], pred_nn)

## Setting levels: control = 0, case = 1

## Warning in roc.default(actual, predicted): Deprecated use a matrix as
## predictor. Unexpected results may be produced, please pass a numeric vector.

## Setting direction: controls < cases

모델링2(R)

서울대학교 과학데이터혁신연구소

2025년 8월

실습

1. 분석을 위해 필요한 라이브러리 및 데이터 불러오기

2. 변수 생성 및 결측 처리

3. 상관관계 분석 및 시각화

4. BMI와 당뇨병 관계 시각화

5. 로지스틱 회귀 (2변수 모형) 및 예측 확률 시각화

6. 로지스틱 회귀 (다변량 모형) 및 잔차 진단

7. 결정 나무 (2변수와 다변량 모델 + 가지치기)

8. 앙상블 모형(배깅, 랜덤 포레스트)

9. 신경망 모형 (Neural Network)

10. 모델 평가 (데이터 분리, Threshold 성능 분석, ROC Curve)

연습문제

1. 로지스틱 회귀 모형

풀이

2. 결정나무

풀이

3. 배깅

풀이

4. 랜덤 포레스트

풀이

5. 인공 신경망

풀이

6. 모델 평가

풀이