[R] 상관 관계 분석, 선형 회귀 분석, 다중 회귀 분석, 로지스틱 회귀 분석

library 설치

install.packages("car")
install.packages("descr")
install.packages('agricolae')
install.packages("nparcomp")
install.packages("survival")
install.packages("Rcmdr")
require(nparcomp)
require(stats)
library(Rcmdr)
library(survival)
library(nparcomp)
library(agricolae)
library(descr)
library(car)

데이터를 받으면 바탕화면으로 빼고, rstudio로 연다. comma로 나뉘어져 있으면 sep만 ,로 바꿔주자.

dat<-read.table(file="C:/Users/cse/Desktop/~~~~.txt",header=TRUE,sep="\t",fill=TRUE)

결측값을 제거 해준다.

# NA 가 있는 행 제거
dat1<-dat1[complete.cases(dat1),]

# NA 값을 0으로 변경
dat1[is.na(dat1)]<-0

# NA 제거
data1<-na.omit(colon)

문제에서 그룹을 나누라고 했으면 그룹을 먼저 나누어준다. ifelse 문으로 나누어줄 수 있다.

dat$grp<-ifelse(dat$group=="ctrl",0,1)

분포를 보고 아웃 라이어를 제거해준다.

# 분포
hist(dat$A)
hist(dat$A, breaks = seq(0, 10000, by=100))

# 산점도
plot(cars)
plot(dist ~ speed, data=cars)
lines(cars)
lines(lowess(cars))

# freq
freq(dat$A)

# 아웃라이어 제거
res <- res[res$A>0 & res$A<10000,]
res <- res[(res$grp!="trt2"),]

모수 검정 vs 비모수 검정

데이터가 5000 이상이면 정규분포에 가까워진다는 중심극한정리.

# shapiro test
shapiro.test(data$weight) ## p-value가 0.05보다 커야 정규분포

데이터가 정규분포하면 모수 검정, 아니면 비모수 검정으로 가자.

모수 검정

• t-test : 범주형(2군) vs 연속형 평균 비교
  • A에 대해 두 그룹 간 유의미한 차이가 있는지?
• ANOVA : 범주형(3군) vs 연속형 평균 비교
  • A에 대해 세 그룹 간 유의미한 차이가 있는지?
  • 사후 검정 : TukeyHSD
• Pearson correlation : 연속형 vs 연속형 상관성
  • A와 B에 상관관계가 있는지?

# bartlett test
bartlett.test(weight~group, data=data) ## 3개 그룹 퍼진 정도(분산) 검사
## p-value가 0.05보다 커야 분산이 같다는 의미.

# ANOVA
out = aov(HT ~ factor(F_BR_3), data=a)
summary(out)
## F값이 0.05보다 작으면 서로 다른 그룹 간 평균 차이가 유의미하다.

TukeyHSD(out)
plot(TukeyHSD(out))

# Pearson correlation

cor.test(cars$speed, cars$dist) # cor 값이 1에 가까울 수록 양의 상관관계에 있다.

# cor 값을 제곱하면 (cor^2) 선형 회귀 분석의 설명력이 된다.

비모수검정

• Mann-Whitney U test : 범주형 (2군) vs 연속형 중위수 비교
  • A에 대해 두 그룹 간 유의미한 차이가 있는지?
• Wilcoxon-Signed rank test : 범주형(2군) vs 연속형 중위수 비교
  • A에 대해 두 그룹 간 유의미한 차이가 있는지?
• kruskal-wallis test : 범주형(3군) vs 연속형 중위수 비교
  • A에 대해 세 그룹 간 유의미한 차이가 있는지?
  • 사후검정 : mctp
• Spearman Correlation : 연속형 vs 연속형 상관성
  • A와 B에 상관관계가 있는지?

# wilcox test
wilcox.test(weight~gr, data=data1, exact=FALSE, conf.level=0.95)

# kruskal test

kruskal.test(weight ~ group, data=data)
result = mctp(weight ~ group, data=data)
summary(result)

## factor를 사용하지 않아도 괜찮다.

# Spearman Correlation

cor.test(cars$speed, cars$dist, method="spearman") # cor 값이 1에 가까울 수록 양의 상관관계에 있다.

선형 회귀 분석

영향을 주는 변수(x축, 독립 변수) : speed

영향을 받는 변수(y축, 종속 변수) : dist

(두 변수 모두 정규분포 해야 함)

우리가 알고 싶은 것? : speed에 따른 dist 변화

# simple regression, 둘 다 continuous일 때

res = lm(dist ~ speed, data=cars)
summay(res)

# simple regression, 영향을 주는 변수가 categorical (2 groups)일 때
res = lm(weight ~ factor(gr), data=data) ## 0과 1만 있다면 factor 안붙여도 됨
summary(res)

# simple regression, 영향을 주는 변수가 categorical (3 groups)일 때
res = lm(weight ~ factor(gr), data=data)
summary(res)
## weight = 5.03  -0.37 * gr1 + 0.49 * gr2
## 결과 분석 팁 : gr이 0일 때가 default라고 하자.
## gr이 1이라면 weight = 5.03 - 0.37
## gr이 2라면 weight = 5.03 + 0.49
## 즉, gr이 0일 때 비해 얼마나 차이가 나는지?

다중 회귀 분석

종속변수(y축) : mpg

독립변수(x축) : wt, vs, gear

res = lm(mpg ~ wt + vs + factor(gear), data=input())
summary(res)

step(res)

## AIC 오름차순으로 표가 나옴. AIC 값이 낮을수록 좋은 것.
## <none> 은 아무것도 안뺐을 때
## 앞에 붙는 +는 추가했을 때, -는 제거했을 때
## 현재보다 더 나은 것이 있으면 자동으로 함수를 실행.

# VIF Variance Inflation Factor : X들간 상관관계 분석
vif(res)

## gvif 값들이 10이하면 괜찮다고 생각. 10이 넘어가면 fomula에서 빼야 함.

# plot 그리기

vif_values <- vif(m)
barplot(vif_values, main="VIF Values", horiz=TRUE, col="steelblue")
abline(v=5, lwd=3, lty=2) # 위에 그린 barplot에 v=5인 점선을 그어준다.

Chi-Square 검정

두 범주형 변수간 관련성이 있는지 판별.

a4<-CrossTable(b$obstruct, b$status)
a4

# chi-square
## 관찰된 빈도가 기대되는 빈도와 의미있게 다른지의 여부를 검정하기 위해 사용되는 검정방법

chisq.test(b$obstruct, b$status)

## p-value가 0.05 이하면 두 변수간 관련성 O

Logistic Regression

종속변수가 categorical 변수일 때 적용. 정규분포를 따르지 않아도 됨.

result = glm(status~ sex + age + obstruct + perfor + adhere + nodes + 
               factor(differ) +
               factor(extent) + surg, family=binomial, data=data1)
summary(result)

step(result)
step(result, direction="backward")
step(result, direction="forward")

Odds Ratio

실패에 비해 성공할 확률의 비.

0일 때 비해 1일 확률의 비

ORtable = function(x, digits=2){
  suppressMessages(a<-confint(x))
  result=data.frame(exp(coef(x)),exp(a))
  result=round(result,digits)
  result=cbind(result, round(summary(x)$coefficient[,4],3))
  colnames(result)=c("OR", "2.5%", "97.5%", "p")
  result
}

ORtable(m)

result = glm(status ~ obstruct  , family=binomial, data=b)
summary(result)
ORtable(result)

# visualization of Odds Ratio
install.packages("moonBook")
library(moonBook)
odds_ratio = ORtable(m)
odds_ratio = odds_ratio[2:nrow(odds_ratio),]
HRplot(odds_ratio, type=2, show.CI=TRUE, cex=2)

해석 : 장이 폐쇄되지 않았을 때 비해서 장이 폐쇄된 경우는 odds ratio배(1.26배)만큼 status확률(재발 또는 사망 확률)이 증가한다. 1에 걸쳐있으면 유의성이 떨어진다고 본다.

해석2 : 암세포가 침습한 깊이가 4(adjacent organ, 장기와 인접)했을 때는 깊이가 1일때 비해 status확률이 4.83배 증가한다.

Logistic 회귀 부분 전체 코드

################################## categorical data analysis ##################################

# dependent var : status ( recur or death of colorectal cancer =1 )
# independent var : 
obstruct : obastruction 종양에 의한 장의 폐쇄

install.packages("survival")
library(survival)
str(colon)
View(colon)

library(descr)
b<-colon

#1. categorical variable : frequency
a1<-freq(b$status)
a1
View(a1)

a2<-freq(b$obstruct)
a2
View(a2)

a3<-table(b$obstruct, b$status)
a3
View(a3)

install.packages("Rcmdr")
library(Rcmdr)
rowPercents(a3)
colPercents(a3)

a4<-CrossTable(b$obstruct, b$status)
a4

#2. chi-sqaure test
chisq.test(b$obstruct, b$status)

# dependent var : status ( recur or death of colorectal cancer =1 )
# independent var : 
obstruct : obastruction 종양에 의한 장의 폐쇄
perfor : perforation 장의 천공
adhere : adherence 인접장기와의 유착
nodes : number of lymphatic gland 암세포가 확인된 림프절 수
differ : 암세포의 조직학적 분화정도(1=well, 2=moderate, 3=poor)
extent : 암세포가 침습한 깊이 (1=submucosa, 2=muscle, 3=serosa, 4=adjacent organ)
surg : time until registration after surgery (0=short, 1=long)

freq(b$obstruct)
freq(b$perfor)
freq(b$adhere)
freq(b$nodes)  # NA's         36
hist(b$nodes)
summary(b$nodes)
hist(b$nodes, breaks=seq(0,35, by=1))
freq(b$differ) # NA's         46
freq(b$extent)
freq(b$surg)

#2. chi-sqaure test
CrossTable(b$obstruct, b$status)
chisq.test(b$obstruct, b$status)

CrossTable(b$perfor, b$status)
chisq.test(b$perfor, b$status)

CrossTable(b$adhere, b$status)
chisq.test(b$adhere, b$status)

CrossTable(b$differ, b$status)
chisq.test(b$differ, b$status)

CrossTable(b$extent, b$status)
chisq.test(b$extent, b$status)

CrossTable(b$surg, b$status)
chisq.test(b$surg, b$status)

# continuous var -> categorical var
b$nodes_gr=ifelse(b$nodes>2,1,0)

CrossTable(b$nodes_gr, b$status)
chisq.test(b$nodes_gr, b$status)

###### logistic regression model  #########################################

# GLM (Genelized Linear Model)

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

dim(b)
data1<-na.omit(colon)
dim(data1)

# sample = data1
result = glm(status~ sex + age + obstruct + perfor + adhere + nodes + differ +
               extent + surg, family=binomial, data=data1)
summary(result)

# differ, extent => gr
result = glm(status~ sex + age + obstruct + perfor + adhere + nodes + 
               factor(differ) +
               factor(extent) + surg, family=binomial, data=data1)
summary(result)

step(result)
step(result, direction="backward")
step(result, direction="forward")
m=step(result, direction="backward")
summary(m)

m=step(result, direction="forward")
m=step(result, direction="both")

# Odds Ratio
ORtable = function(x, digits=2){
  suppressMessages(a<-confint(x))
  result=data.frame(exp(coef(x)),exp(a))
  result=round(result,digits)
  result=cbind(result, round(summary(x)$coefficient[,4],3))
  colnames(result)=c("OR", "2.5%", "97.5%", "p")
  result
}

ORtable(m)

# sex + age + obstruct + perfor + adhere + nodes + differ + extent + surg

result = glm(status ~ obstruct  , family=binomial, data=b)
summary(result)
ORtable(result)

result = glm(status ~ adhere , family=binomial, data=b)
summary(result)
ORtable(result)

result = glm(status ~ nodes , family=binomial, data=b)
summary(result)
ORtable(result)

result = glm(status ~ factor(extent) , family=binomial, data=b)
summary(result)
ORtable(result)

result = glm(status ~ surg , family=binomial, data=b)
summary(result)
ORtable(result)

# visualization of Odds Ratio
install.packages("moonBook")
library(moonBook)
odds_ratio = ORtable(m)
odds_ratio = odds_ratio[2:nrow(odds_ratio),]
HRplot(odds_ratio, type=2, show.CI=TRUE, cex=2)