생존분석 패키지 survival, survminer, muhaz 설치 및 예제

   필자가 근무하고 있는 부서는  보전(maintenance Dept') 부서이다.  생존분석이  매우 중요하다.   고장확률은 아래와 같다. 

    고장확률 =  1 - 생존확를  

  센서에서 고장에 대한 Event 알람 신호를 받아,  알람 신호기간 동안 고장 누적 분포 함수를 구할 수 있다.   

     아래 데이터는 Worcester Heart Attack Study  건으로  심장마비 후 생존 시간에 영향을 미칠 수 있는 연령, 성별 BMI 등 여러가지 요인 조사하였다.  모든 참가자는 심장마비 후 병원 입원시 시작되어 후속 조치에 따른 손실 또는 사망으로 끝난다.    

 아래 사용된 변수는 아래와 같다. 

• lenfol :  추적기간 (사망 또는  후속조치에 따른 중단)
• fstat :  후속조치에 다른 중단 변수 ( 손실 = 0 ,   사망 = 1)
• age : 병원에 있을때 나이
• hr : 초기 심박수
• gender : 남성 = 0, 여성 =1

1. survaival packages   

• 설치 : R base 4.05 버전에 내장 되어 있다.
• 사용 : 생존 분석 모델링 ( Kaplan-Meier 분포) 

2. survminer packages

• 설치 : install.packages("survminer")   
• 사용 : ggplot2와 연계되어 생존 분석 그래프 표현

3. muhaz packages

• 설치 : install.packages("muhaz")
• 사용 : 위험 함수 모델링

4. 생존분석 package 예제 

• 생존 데이터 로딩 후에  survaival 패키지를 이용한  Kaplan-Meier 분포를 모델링 하여, survminer 패키지로 그래프 그린다. 

# 라이브러리 불러오기 
  library(rio)  
  library(survival)
  library(survminer)
  library(muhaz)
  library(dplyr)

library(readxl)
library(httr)

# 파일을 불러오는 path를 url 변수로 정의      

url = "https://docs.google.com/uc?export=download&id=1O8on1RzEE_-CSi-9zYHR7Hk_lMOCGxo4"
GET(url, write_disk(tf <- tempfile(fileext = ".xlsx")))

whass <- read_excel(tf, 1L)      

# dplyr의 select 함수를 이용하여 필요한 변수만 선택
  
  whass <- whass %>% 
            select(lenfol, fstat, age, hr, gender)
  head(whass)      

##   lenfol fstat age hr gender
## 1   2178     0  83 89      0
## 2   2172     0  49 84      0
## 3   2190     0  70 83      1
## 4    297     1  70 65      0
## 5   2131     0  70 63      0
## 6      1     1  70 76      0

# kaplan meier test
  s_fit <- survfit(Surv(lenfol, fstat)~ gender, data = whass)
  ggsurvplot(s_fit, surv.median.line = "hv", risk.table = T, 
             main = "Survival Curve with Censored Data")

  생존 함수의 그래프를 분석해 보면,  처음에는 남자 300명, 여자200명이 들어 온다.  500시간 되는 시점에 약 60%의 인원이 생존 되다가,   1000시간 까지 어느 정도 유지 되다가  급격히 사망인원이 증가 한다.   이것은 생존에 대한 것만  보일 뿐 어느 시간대가 위험 한지 알 수 없다. 

# hazard ratio fstat =1  
  whass_1 <- whass %>%  filter(fstat == 1)
  fit_haz_nc <- muhaz(whass_1$lenfol, whass_1$fstat, bw.method ="g", bw.grid=200)
  plot(fit_haz_nc, main="Hazard Rate with No Censored Data")

위에 있는 위험 함수를 보면 초기에 사망률이 높다가, 1400시간 이후 급격하게 사망률이 높아지는 것을 볼 수 있다.  

   설비관리에서는 욕조 커브라고 부른다.  설비 처음에 들어 올 때는 설비 유동고장이 많고,  교체 직전에는 마모 고장이 많다는 것이다.  생존분포의 위험 함수는 즉 고장위험 함수와 맥락이 같다고 보면 된다.