로지스틱 회귀분석-환자들 생존예측 (예제 따라하기)

해당 코드와 각주의 설명은 아래 유튜브 영상을 토대로 작성했습니다.

https://www.youtube.com/watch?v=Fgl2cOXpuyM&t=608s

 

혼자서 독학을 하려고 하다보면, 특히 해당 전공분야가 아닐경우 더더욱 어려움을 겪게 되는데

이런 강의들 너무 소중합니다.

 

예제에 사용된 데이터는 Survival이라는 패키지에 담겨있는 colon데이터 입니다.

 

각주는 기본적으로 유튜브 강의 내 설명을 기초로 하고, 제가 생각한 코드들이 합쳐져 있습니다.

require(survival)
str(colon)
#status를 종속변수로 두려고 함. 0일 경우 생존, 1일 경우 사망, 재발 등을 뜻함.
colon1<-na.omit(colon) #NA값들을 제거함. 하지만 함부로 제거해서는 안되는 경우가 있음. 
#구조적인 데이터 결측의 가능성을 배제하면 안됨.
View(colon)
dim(colon)
dim(colon1)

result<-glm(status~rx+sex+age+obstruct+perfor+adhere+nodes+differ+extent+surg, family = binomial, data=colon1)
summary(result)

#회귀식을 그대로 받아들여서는 안됨.
#모델 선택에 대한 적용
#변수를 제거해가며 AIC확인하기
Reduced.model<-step(result)

summary(Reduced.model)
#보고를 할 때 회귀식을 이용하기보다 오즈레이쇼?를 중심으로 사용
require(moonBook)
#moonbook안에 extractOR이라는 함수를 이용
extractOR(Reduced.model) #95%신뢰구간에 있어서 하한치, 상한치, p-value값이 나옴.
#해석 : obstruct의 경우 OR가 1.25. 하한치와 상한치 값이 1을 포함하기 때문에 p-value는 의미 없음.
#nodes는 하한치와 상한치가 1을 포함하지 않음. p-value역시 굉장히 작은 값.

#이 결과를 받아들이기 전에 과산포(?)를 확인해야함.
#로지스틱 회귀분석 할 때, family는 늘 binomial을 쓰는 것은 아님. 과산포를 염두.
#과산포(Overdispersion)가 있는지의 유무. 예측모델에서 잔차가 얼마나 되는지에 따라 
#family옵션 중 (binomial,quasibinomial) 선택하여 사용해야함.
#과산포가 없다면 binomial. 과산포가 있다면 quasibinomial)를 선택.
#따라서 과산포를 확인하기 위해서는 family함수를 이용해 2개의 모델을 만들어야함.
#pchisq() 함수를 이용해서 과산포유무를 확인할 수 있음.

fit = glm(formula= status~rx+obstruct+adhere+nodes+extent+surg, family = binomial, data=colon1)
fit.od = glm(formula= status~rx+obstruct+adhere+nodes+extent+surg, family = quasibinomial, data=colon1)

#fit.od안에 있는 dispersion값과 fit모델에서 자유도 값을 가져오 
pchisq(summary(fit.od)$dispersion*fit$df.residual,fit$df.residual,lower=F)
#0.2803691 값이 0.05보다 크다면 과산포는 없다고 확신할 수 있습니다. 
#따라서 현재 여기서는 과산포가 없으므로 binomial 옵션을 사용하면 됨.

#예측모델을 함수 안에 바로 넣어보기
#moonbook 패키지에 있는 ORplot
ORplot(fit)
#값들이 1을 포함하지 않으면 의미가 있다는 것을 확인할 수 있음.

ORplot(fit, main="Plot for Odds ratios")

#SHOW.OR은 그래프 안에 값을 보여줄지 않을지. SHOW.CI는 옆에 값을 출력할지 안할지
ORplot(fit, type=2, show.OR = FALSE, show.CI=TRUE, pch=10, lwd=3, col=c("blue","red"), main="Plot of OR")


ORplot(fit, type=1, show.OR = FALSE, show.CI=TRUE, pch=10, lwd=3, col=c("blue","red"), main="Plot of OR")