기타문서

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

Mosu(정종인) 2022. 12. 12. 17:12
반응형

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)
반응형
저작자표시 (새창열림)