R을 활용한 외부데이터 분석

이번 게시글에선 주어진 txt파일을 분석해보겠습니다.

3개의 txt파일에 대한 정보입니다.

 

## ex1.txt  
# 대상자 번호  pid 
# 당뇨측정시점  redcap_event_name  
# 당화혈색소  HbA1c 
# HbA1c(당화혈색소) ≥ 6.5% : 당뇨 정의 

## ex2.txt  
# 체중_대상자 번호 BW_myhealth_id 
# 체중_측정 일시 BW_date_measured 
# 체중_측정값 BW_value 

## ex3.txt  
# 대상자 번호 pid 
# 대상자 나이 age  
# 대상자 성별 sex

 

1. 데이터 불러오기

input1 = read.table("exam302_1.txt", header=T, sep="\t")
input2 = read.table("exam302_2.txt", header=T, sep="\t")
input3 = read.table("exam302_3.txt", header=T, sep="\t")

3개의 파일을 각각 다른 변수에 할당합니다.

 

2. 데이터 전처리

input1 = input1[,c("pid","hba1c")]
input2 = input2[,c("pid","BW_value")]
input3 = input3[,c("pid","age")]

분석에 필요한 변수만 추출합니다.

 

data = merge(input3,input1,by="pid")
data = merge(data,input2,by="pid")

대상자번호(pid) 기준으로 하나의 테이블로 합칩니다.

 

data = data[data$"BW_value"<300 & data$"BW_value">0,]

체중이 비정상정으로 크거나 0인 행을 제거합니다.

 

data = na.omit(data)

결측값을 제거합니다.

 

mainData = data[!duplicated(data$pid),]
mainData = aggregate(cbind(BW_value, hba1c) ~ pid+age, data, mean)

pid 기준, 중복되는 행들이 많습니다. 체중과 당화혈색소를 평균내서 중복을 제거합니다.

 

 

가설 1 : 당뇨인 군과 당뇨가 아닌 군의 평균체중은 차이가 있다.

1. 데이터 불러오기

ret1 = mainData
ret1$"hba1c" = ifelse(ret1$"hba1c">=6.5,1,0)

가설1의 결과를 담을 변수 ret1에 전처리된 데이터를 저장합니다.

또한, 당뇨여부를 판단해야하므로 hba1c 열을 당뇨여부를 표시하도록 변경합니다.

 

2. 정규분포성 검사

qqnorm(ret1$BW_value[ret1$BW_value=="0"])
qqnorm(ret1$BW_value[ret1$BW_value=="1"])

..

샘플이 5000개보다 많으므로, 정규분포를 가정합니다.

 

3. 분산 비교

var.test(BW_value~hba1c,ret1)

var.test 결과 p-value > 0.05 이므로, 분산이 같다는 결과를 얻습니다.

 

4. 평균 비교

t.test(BW_value~hba1c,ret1,var.equal=T)

t.test 결과 p-value > 0.05 이므로, 평균이 같다는 결과를 얻습니다.

 

5. 함의

당뇨 여부와 평균체중은 통계적으로 유의한 차이가 없다.

두 집단의 평균체중 차이는 약 0.2 정도이다.

 

 

가설 2 : 당뇨환자의 연령그룹이 (,40) [40,60) [60,) 라면, 그룹별 평균체중은 차이가 있다.

1.  데이터 불러오기

ret2 = mainData[mainData$"hba1c">=6.5,]

가설2의 결과를 담을 변수 ret2에 전처리된 데이터를 저장합니다.

 

2.  정규분포성 검사

샘플이 5000개보다 많으므로, 정규분포를 가정합니다.

 

3. 등분산 검사

bartlett.test(BW_value~factor(gp),ret2)

bartlett.test 결과 p-value < 0.05 이므로, 분산이 다르다는 결과를 얻었습니다.

 

4. 평균 비교

oneway.test(BW_value~factor(gp), ret2, var.equal=F)

oneway.test 결과 p-value < 0.05 이므로 평균이 다르다는 결과를 얻었습니다.

 

5. 사후검정

pairwise.t.test(ret2$BW_value, ret2$gp, p.adjust.method="bonferroni")

해당 테스트 결과, 2,3 그룹에서 통계적으로 유의한 차이가 보였습니다.

 

만약 등분산이었다면..

result = aov(BW_value ~ gp, ret2)
plot(TukeyHSD(result))

그래프에서, 0을 포함하지 않는다면 통계적으로 유의한 차이를 보인 것이다.