'KNOU' 카테고리의 글 목록 (2 Page)

제6.2장 비모수적 방법 - 2. 누적한계추정법

KNOU/2 보건 정보 데이터 분석 2017. 1. 26. 09:37

비모수적 방법 non-parametric method

6.2.2 누적한계추정법 product-limit method

생명표 방법: 기간을 1년, 6개월 등 특정 단위로 나눠 구분

누적한계추정법: 매 사건이 발생할 때 마다 해당 시점을 표기

누적한계추정법(product-limit method)은 생존함수를 추정하는 대표적인 방법 중 하나로 연구자의 이름을 따서 Kaplan-Meier추정법이라고도 한다. 모든 환자의 생존시간 또는 중도절단 시간이 각각 관찰되었다고 하자. 모든 자료의 생존시간 도는 중도절단 x1,…,투을 순서대로 배열한 것을 t1<t2<…<tn이라 하고, δi = 0, 그렇지 않은 경우 δi =1로 정의한다. 즉 중도 절단 되지 않고 사건이 발생한 경우 status = 1, 중도절단된 경우 stats = 0 로 표기한다.

예] 신장이식수술을 받은 15명 환자들의 호전기간(remission duration)이 아래와 같다고 하자. 여기서 +로 표기된 것은 중도절단된 자료를 가르킨다. 각 시점에서 생존함수를 누적한계추정법을 이용하여 추정해보자

표. 신장이식 환자들의 호전기간 (단위 : 일)

3.0 4.0+ 4.5 4.5 5.5 6.0 6.4 6.5 7.0 7.5 8.4+ 10.0+ 10.0 12.0 15.0

1. 데이터 불러오기

> library(survival)

> setwd('c:/Rwork')

> kidney<-read.table('kidney.txt',header=T)

> kidney #status =1 사건, 0 = 절단, 절단 표시 매우 중요

time status

1 3.0 1

2 4.0 0

3 4.5 1

4 4.5 1

5 5.5 1

6 6.0 1

7 6.4 1

8 6.5 1

9 7.0 1

10 7.5 1

11 8.4 0

12 10.0 0

13 10.0 1

14 12.0 1

15 15.0 1

> attach(kidney)

2. 누적한계추정치(Kaplan-Meier 추정치)

> fit1 = survfit(Surv(time,status)~1, data=kidney) #~ 우측에는 공변량

> summary(fit1)

Call: survfit(formula = Surv(time, status) ~ 1, data = kidney)

time n.risk n.event survival std.err lower 95% CI upper 95% CI

3.0 15 1 0.933 0.0644 0.8153 1.000

4.5 13 2 0.790 0.1081 0.6039 1.000

5.5 11 1 0.718 0.1198 0.5177 0.996

6.0 10 1 0.646 0.1275 0.4389 0.951

6.4 9 1 0.574 0.1320 0.3660 0.901

6.5 8 1 0.503 0.1336 0.2984 0.846

7.0 7 1 0.431 0.1324 0.2358 0.787

7.5 6 1 0.359 0.1283 0.1781 0.723

10.0 4 1 0.269 0.1237 0.1094 0.663

12.0 2 1 0.135 0.1135 0.0258 0.703

15.0 1 1 0.000 NaN NA NA

> fit1

Call: survfit(formula = Surv(time, status) ~ 1, data = kidney)

n events median 0.95LCL 0.95UCL

15 12 7 6 NA

Error: invalid multibyte character in parser at line 1

-> fit1 = survfit(Surv(time,status)~1, data=kidney) #~ 우측에는 공변량

6.5일 에서의 생존확률은 50.3%, 7.0일 이 지나면 생존율이 43.1%로 50%이하가 된다.

3. 사망시점의 사분위수 추정치와 그의 신뢰구간, 가령 사망자가 50%되는 시점은 언제인가?

> quantile(fit1,probs=c(0.25,0.5,0.75),conf.int=T) #신뢰구간까지

$quantile

25 50 75

5.5 7.0 12.0

$lower

25 50 75

4.5 6.0 7.0

$upper

25 50 75

7.5 NA NA

-> 생존확률이 75%인 경우, 즉 사망확률이 25%인 경우의 시점은 5.5일, 생존확률이 50%인 경우 7.0일, 생존확율이 25%인 경우의 시점은 12.0 일

4. 누적한계추정치의 95% 신뢰구간 그래프

> plot(fit1,xlab='time',ylab='Survival function',lwd=2)

> legend(0.2,0.3,c('KM estimate','95% CI'),lty=c(1,2))

-> 점선은 신뢰구간을 의미, 실선은 생존함수추정치.생존확율이 50%인 경우는 약 7.0일임을 알 수 있다.

출처: 보건정보 데이터분석 (이태림, 이재원, 김주한, 장대흥 공저)

'KNOU > 2 보건 정보 데이터 분석' 카테고리의 다른 글

제6.4장 모수적 방법 (0)	2017.01.26
제6.3장 비모수적 방법을 이용한 생존함수의 비교 (0)	2017.01.26
제6.2장 비모수적 방법 - 1. 생명표 방법 (0)	2017.01.26
제5.1장 공분산 분석 - 2 공변량이 둘 이상인 경우 (2)	2017.01.13
제5.1장 공분산 분석 - 1. 공변량이 하나인 경우 (0)	2017.01.06

Posted by 마르띤

,

제6.2장 비모수적 방법 - 1. 생명표 방법

KNOU/2 보건 정보 데이터 분석 2017. 1. 26. 09:36

비모수적 방법 non-parametric method

6.2.1생명표 방법 life table method

생명표 방법 개념 추가 설명

예] 협심증(angina pectoris)이 있는 2,418명의 남성들에 대한 생존자료이다. 이 생명표에서 경과기관에 따른 생존확률을 구해보자

진단 후 경과 기관(단위: 년)	사망자 수	중도절단 수
( 0 – 1 ]	456	0
( 1 – 2 ]	226	39
( 2 – 3 ]	152	22
( 3 – 4 ]	171	23
( 4 – 5 ]	135	24
( 5 – 6 ]	125	107
( 6 – 7 ]	83	133
( 7 – 8 ]	74	102
( 8 – 9 ]	51	68
( 9 – 10 ]	42	64
( 10 – 11 ]	43	45
( 11 – 12 ]	34	53
( 12 – 13 ]	18	33
( 13 – 14 ]	9	27
( 14 – 15 ]	6	23
( 15 –	0	30

1. 데이터 입력

> library(KMsurv)

> setwd('C:/Rwork')

> 협심증환자자료 = read.csv('협심증환자자료.csv')

> head(협심증환자자료)

time censor freq

1 0.5 1 456

2 0.5 0 0

3 1.5 1 226

4 1.5 0 39

5 2.5 1 152

6 2.5 0 22

> attach(협심증환자자료)

2. lifetab 함수 연구

> lifetab

function (tis, ninit, nlost, nevent)

(이하 생략)

-> tis: 시간, 길이는 17, ninit = 전체 데이터, 2418명, nlost = 중도절단, nevent 사건발생

3. nevent 자료 만들기

> which(censor==1) #사망 데이터

[1] 1 3 5 7 9 11 13 15 17 19 21 23 25 27 29 31

> 집단.사망=협심증환자자료 [which(censor==1),]

> head(집단.사망)

time censor freq

1 0.5 1 456

3 1.5 1 226

5 2.5 1 152

7 3.5 1 171

9 4.5 1 135

11 5.5 1 125

4. nlost 자료 만들기

> which(censor==0) #절단 데이터

[1] 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 32

> 집단.절단=협심증환자자료[which(censor==0),]

> head(집단.절단)

time censor freq

2 0.5 0 0

4 1.5 0 39

6 2.5 0 22

8 3.5 0 23

10 4.5 0 24

12 5.5 0 107

5. ninit전체데이터 자료 만들기

> 사망자수=집단.사망[,3]

> 사망자수

[1] 456 226 152 171 135 125 83 74 51 42 43 34 18 9 6 0

> 절단자수=집단.절단[,3]

> 절단자수

[1] 0 39 22 23 24 107 133 102 68 64 45 53 33 27 23 30

> 합=사망자수+절단자수

> 합

[1] 456 265 174 194 159 232 216 176 119 106 88 87 51 36 29 30

> sum(합) #2418

[1] 2418

6. tis 시간 변수 자료 만들기

> 년 = floor(집단.사망$time) #floor 내림값

> 년

[1] 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15

> lt=length(년)

> lt

[1] 16

> length(년)

[1] 17

> 년[lt+1] = NA

> 년[lt+1]

[1] NA

7. 생명표 만들기

> 생명표=lifetab(년,sum(합),절단자수,사망자수)

> 생명표

nsubs nlost nrisk nevent surv pdf hazard se.surv

0-1 2418 0 2418.0 456 1.0000000 0.18858561 0.20821918 0.000000000

1-2 1962 39 1942.5 226 0.8114144 0.09440394 0.12353102 0.007955134

2-3 1697 22 1686.0 152 0.7170105 0.06464151 0.09440994 0.009179397

3-4 1523 23 1511.5 171 0.6523689 0.07380423 0.11991585 0.009734736

4-5 1329 24 1317.0 135 0.5785647 0.05930618 0.10804322 0.010138361

5-6 1170 107 1116.5 125 0.5192585 0.05813463 0.11859583 0.010304216

6-7 938 133 871.5 83 0.4611239 0.04391656 0.10000000 0.010379949

7-8 722 102 671.0 74 0.4172073 0.04601094 0.11671924 0.010450930

8-9 546 68 512.0 51 0.3711964 0.03697464 0.10483042 0.010578887

9-10 427 64 395.0 42 0.3342218 0.03553750 0.11229947 0.010717477

10-11 321 45 298.5 43 0.2986843 0.04302654 0.15523466 0.010890741

11-12 233 53 206.5 34 0.2556577 0.04209376 0.17941953 0.011124244

12-13 146 33 129.5 18 0.2135639 0.02968456 0.14937759 0.011396799

13-14 95 27 81.5 9 0.1838794 0.02030570 0.11688312 0.011765989

14-15 59 23 47.5 6 0.1635737 0.02066194 0.13483146 0.012259921

15-NA 30 30 15.0 0 0.1429117 NA NA 0.013300258

se.pdf se.hazard

0-1 0.007955134 0.009697769

1-2 0.005975178 0.008201472

2-3 0.005069200 0.007649121

3-4 0.005428013 0.009153696

4-5 0.004945997 0.009285301

5-6 0.005033980 0.010588867

6-7 0.004690538 0.010962697

7-8 0.005175094 0.013545211

8-9 0.005024599 0.014659017

9-10 0.005307615 0.017300846

10-11 0.006269963 0.023601647

11-12 0.006847514 0.030646128

12-13 0.006682743 0.035110295

13-14 0.006514794 0.038894448

14-15 0.008035120 0.054919485

15-NA NA NA

nsubs	nlost	nrisk	nevent	surv
생존자수	중도절단수	유효인원수	사망자수	생존함수

pdf	hazard	se.surv	se.pdf	se.hazard
확률밀도 함수	위험함수	생존함수의 표준오차	확률밀도 함수의 표준오차	위험함수의 표준오차

-> 5번째 열 surv는 생존 함수, 7번째 열 hazard는 위험 함수. 협심증 환자의 약 19%가 1년 이내, 약 28%가 2년 이내에 사망하는 것으로 추정된다. 따라서 협심증 환자가 2년 이상 생존할 확률은 약 72%가 된다는 것을 알 수 있다. 또한 5년 이상 생존할 확률은 약 52%이다.

#전체 환자 중 생존확율이 70% 되는 지점

> names(생명표)

[1] "nsubs" "nlost" "nrisk" "nevent" "surv" "pdf"

[7] "hazard" "se.surv" "se.pdf" "se.hazard"

> which.min(abs(생명표$surv-0.7))

[1] 3

> 생명표[3,]

nsubs nlost nrisk nevent surv pdf hazard se.surv

2-3 1697 22 1686 152 0.7170105 0.06464151 0.09440994 0.009179397

se.pdf se.hazard

2-3 0.0050692 0.007649121

-> abs:절대값, which.min 최소값, 또는 생명표에서 3번째 행을 보면 survival이 70%에 가장 가까이 접근했음을 알 수 있다. 2년 이상 생존할 확률은 72%

#전체 환자 중 생존확율이 50% 되는 지점

> which.min(abs(생명표$surv-0.5))

[1] 6

> 생명표[6,]

nsubs nlost nrisk nevent surv pdf hazard se.surv

5-6 1170 107 1116.5 125 0.5192585 0.05813463 0.1185958 0.01030422

se.pdf se.hazard

5-6 0.00503398 0.01058887

-> 5년이 지나면 surv 즉 생존할 확률이 약 52%가 됨을 알 수 있다.

8. 데이터 시각화 – 생존함수 추정치 그래프

> plot(년[1:lt], 생명표[,5],type='s',xlab='year',ylab='Survival function',ylim=c(0,1),lwd=2)

> abline(h=0.5)

> abline(v=5)

> plot(년[1:lt], 생명표[,5],type='o',xlab='year',ylab='Survival function',ylim=c(0,1),lwd=2)

> abline(h=0.5)

> abline(v=5)

9. 데이터 시각화 – 위험함수 추정치 그래프

> names(생명표)

[1] "nsubs" "nlost" "nrisk" "nevent" "surv" "pdf"

[7] "hazard" "se.surv" "se.pdf" "se.hazard"

> mean(생명표$hazard)

[1] NA

> 생명표$hazard

[1] 0.20821918 0.12353102 0.09440994 0.11991585 0.10804322 0.11859583

[7] 0.10000000 0.11671924 0.10483042 0.11229947 0.15523466 0.17941953

[13] 0.14937759 0.11688312 0.13483146 NA

> mean(생명표$hazard,na.rm=T) #NA 값 제외한 평균

[1] 0.1294874

> plot(년[1:lt], 생명표[,7],type='s',xlab='년',ylab='Hazard function',ylim=c(0,0.25),lwd=2)

> abline(h=0.1294874)

> plot(년[1:lt], 생명표[,7],type='o',xlab='년',ylab='Hazard function',ylim=c(0,0.25),lwd=2)>

> abline(h=0.1294874)

출처: 보건정보 데이터분석 (이태림, 이재원, 김주한, 장대흥 공저)

'KNOU > 2 보건 정보 데이터 분석' 카테고리의 다른 글

제6.3장 비모수적 방법을 이용한 생존함수의 비교 (0)	2017.01.26
제6.2장 비모수적 방법 - 2. 누적한계추정법 (0)	2017.01.26
제5.1장 공분산 분석 - 2 공변량이 둘 이상인 경우 (2)	2017.01.13
제5.1장 공분산 분석 - 1. 공변량이 하나인 경우 (0)	2017.01.06
제4장 범주형 자료의 분석 - 4.3 로짓분석 (0)	2017.01.04

Posted by 마르띤

,

제5.1장 공분산 분석 - 2 공변량이 둘 이상인 경우

KNOU/2 보건 정보 데이터 분석 2017. 1. 13. 10:06

[예] 어떤 약물에 대한 체내 배출연구에서 얻은 자료이다. 연구자는 약의 형태에 따라 체내로부터 배출되는 약물의 양이 달라지는지를 알고자 한다. 그런데 배출되는 약물의 양은 약의 형태뿐만 아니라 배출된 약물을 측정한 시간과 각 개체의 항신진대사 점수에도 영향을 받을 것으로 생각한다. 이러한 경우에는 측정시간과 항신진대사 점수를 2개의 공변량으로 하여 이들을 제어한 약의 형태에 대한 효과를 공분산분석을 통해 알 수 있다.

관측번호	약의형태(trt)	항신진대사점수(x1)	소요시간(x2)	약물량(y)
1	1	37	61	11.3208
2	2	37	37	12.9151
3	3	45	53	18.8947
…	…	…	…	…

공변량이 2개 이상인 경우에는 공변량이 하나인 경우의 모형을 그대로 확장해서 각 모수의 추정과 검정을 할 수 있다.

<공분산분석을 위한 두 가지 가정>

1) 각 처리 안에서 반응변수Y에 미치는 공변량x의 효과가 모두 동일해야 한다. 즉, 회귀계수가 모든 약의 형태에 대해서 동일해야 하며, 교호작용이 없어야 한다.

2) 공변량x 효과가 0이 아니다. 효과가 0이라면 분산분석을 하면 된다.

2. 공분산분석에서의 검정

1) H0: β1 = 0 항신진대사의 효과가 없음

2) H0: β2 = 0 소요시간의 효과가 없음

이상 2개의 공변량 효과를 제어한 후 배출된 약물량의 모평균이 약의 형태type에 따라 차이가 있는가를 검정할 수 있는데, 이를 검정하기 위한 귀무가설은 아래와 같다.

3) H0 : α1 = α2 = ... = αI

1. library 호출 및 데이터 불러오기

> library(lsmeans)

> setwd('C:/Rwork')

> drug = read.csv('약물배출량자료.csv')

> head(drug)

관측번호 약형태 항신진대사점수 소요시간 약물량

1 1 1 37 61 11.3208

2 2 2 37 37 12.9151

3 3 3 45 53 18.8947

4 4 4 41 41 14.6739

5 5 5 57 41 8.6493

6 6 6 49 33 9.5238

> colname<-c('obs','type','x1','x2','y')

> colnames(drug)<-colname

> head(drug)

obs type x1 x2 y

1 1 1 37 61 11.3208

2 2 2 37 37 12.9151

3 3 3 45 53 18.8947

4 4 4 41 41 14.6739

5 5 5 57 41 8.6493

6 6 6 49 33 9.5238

> attach(drug)

2. 회귀계수의 동일성 검정(교호작용 존재 확인)

공분산 분석을 위한 중요한 가정으로 교호작용이 존재하는지 점검한다. 교호작용이 존재하면 공분산분석이 아닌 분산분석을 실시한다.

> model1 = aov(y~factor(type) + factor(type)*x1 + factor(type)*x2 , data = drug)

> summary(model1)

Df Sum Sq Mean Sq F value Pr(>F)

factor(type) 5 250.3 50.1 2.099 0.113

x1 1 696.6 696.6 29.206 3.9e-05 ***

x2 1 54.4 54.4 2.282 0.148

factor(type):x1 5 160.8 32.2 1.349 0.289

factor(type):x2 5 42.3 8.5 0.355 0.872

Residuals 18 429.3 23.9

---

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘

-> 공변량이 2개 이므로 우변에 * 형태로 입력

귀무가설 H0 : 교호작용이 존재하지 않는다.

대립가설 H1 : 교호작용이 존재한다.

p-value: x1 항신진대사점수 – 0.289

x2 소요시간 – 0.872

의사결정: 유의수준 5% 하에서 귀무가설을 기각할 수 없다.

결론: 두 공변량(항신진대사점수, 소요시간)과 약의 형태 사이에 교호작용이 존재하지 않으므로 공분산분석을 할 수 있다. ( = 각 처리 안에서 반응변수Y에 미치는 공변량x의 효과가 모두 동일하다.)

3. 이원공분산분석 (Two-way ANOVA)

> model2 = aov(y~factor(type) + x1 + x2, data=drug)

> summary(model2)

Df Sum Sq Mean Sq F value Pr(>F)

factor(type) 5 250.3 50.1 2.216 0.0808 .

x1 1 696.6 696.6 30.837 6.13e-06 ***

x2 1 54.4 54.4 2.409 0.1319

Residuals 28 632.5 22.6

---

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

-> 2개의 공변량을 보정하지 않았을 때의 F-값과 p값은 각각 2.216, 0.0808으로 유의수준 5%하에서 유의하지 않다.

귀무가설 H0 : α1 = α2 = ... = αI

대립가설 H1 : 최소한 하나 이상의 약은 효과가 있다.

p-value: 0.0808

의사결정: 유의수준 5% 하에서 귀무가설을 기각할 수 없다.

결론: 약의 형태에 따라 배출된 약물량의 차이가 없다.

3-1. 공변량 효과 제어 시 치료법의 효과 검정

> model3 = lm(y~factor(type) + x1 + x2, data=drug)

> summary(model3)

Call:

lm(formula = y ~ factor(type) + x1 + x2, data = drug)

Residuals:

Min 1Q Median 3Q Max

-7.222 -2.634 -0.379 1.475 9.646

Coefficients:

Estimate Std. Error t value Pr(>|t|)

(Intercept) 36.6370 8.9591 4.089 0.000331 ***

factor(type)2 0.8965 2.7823 0.322 0.749677

factor(type)3 7.9097 2.7684 2.857 0.007973 **

factor(type)4 3.0722 2.8247 1.088 0.286025

factor(type)5 9.5434 2.8534 3.345 0.002355 **

factor(type)6 5.8389 2.8391 2.057 0.049149 *

x1 -0.7606 0.1375 -5.531 6.51e-06 ***

x2 0.1647 0.1061 1.552 0.131868

---

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 4.753 on 28 degrees of freedom

Multiple R-squared: 0.6129, Adjusted R-squared: 0.5161

F-statistic: 6.332 on 7 and 28 DF, p-value: 0.0001641

귀무가설 H0 : β1 = β2 = 0 공변량(항신진대사, 소요시간)의 효과가 없음

H0 : α1 = α2 = ... = αI

대립가설 not H0

p-value : 0.0001641

의사결정: 귀무가설을 강하게 기각한다. 5%의 유의수준 하에서 매우 유의하다.

결론: 우리가 세운 모형의 자료에 적합하다는 것을 알 수 있다. 배출된 약물량의 모평균이 약 형태type의 효과와 공변량들의 효과가 없다라고 말할 수 없다.

4. 공변량 효과 제어시 치료법의 효과 검정 - 모형 제곱합

1)제1종 제곱합(Type I SS): SS(type)

> model2 = aov(y~factor(type) + x1 + x2, data=drug)

> summary(model2)

Df Sum Sq Mean Sq F value Pr(>F)

factor(type) 5 250.3 50.1 2.216 0.0808 .

x1 1 696.6 696.6 30.837 6.13e-06 ***

x2 1 54.4 54.4 2.409 0.1319

Residuals 28 632.5 22.6

---

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

OR

> summary(aov(y~factor(type)+x1+x2,data=drug))

Df Sum Sq Mean Sq F value Pr(>F)

factor(type) 5 250.3 50.1 2.216 0.0808 .

x1 1 696.6 696.6 30.837 6.13e-06 ***

x2 1 54.4 54.4 2.409 0.1319

Residuals 28 632.5 22.6

---

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’

-> 결과해석: 약의 형태type가 기여한 제1종 제공합의 p-value는 0.0808로 유의수준 5%하에서 유의하지 않다.

2)제3종 제곱합(Type III SS): SS(type | x1, x2)

> summary(aov(y~x1+x2+factor(type),data=drug))

Df Sum Sq Mean Sq F value Pr(>F)

x1 1 516.5 516.5 22.867 5.03e-05 ***

x2 1 62.8 62.8 2.779 0.1067

factor(type) 5 422.0 84.4 3.736 0.0102 *

Residuals 28 632.5 22.6

---

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘

-> 결과해석: x1,2 given type, 공변량x1,x2가 기여한 상태에서 순수하게 약의 형태type가 기여한 제3종 제공합의 p-value는 0.0102로 유의수준 5%하에서 유의하다.

5. 각 수준 별 추정치를 알아보자.

> model3 = lm(y~factor(type) + x1 + x2, data=drug)

> summary(model3)

Call:

lm(formula = y ~ factor(type) + x1 + x2, data = drug)

Residuals:

Min 1Q Median 3Q Max

-7.222 -2.634 -0.379 1.475 9.646

Coefficients:

Estimate Std. Error t value Pr(>|t|)

(Intercept) 36.6370 8.9591 4.089 0.000331 ***

factor(type)2 0.8965 2.7823 0.322 0.749677

factor(type)3 7.9097 2.7684 2.857 0.007973 **

factor(type)4 3.0722 2.8247 1.088 0.286025

factor(type)5 9.5434 2.8534 3.345 0.002355 **

factor(type)6 5.8389 2.8391 2.057 0.049149 *

x1 -0.7606 0.1375 -5.531 6.51e-06 ***

x2 0.1647 0.1061 1.552 0.131868

---

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 4.753 on 28 degrees of freedom

Multiple R-squared: 0.6129, Adjusted R-squared: 0.5161

F-statistic: 6.332 on 7 and 28 DF, p-value: 0.0001641

-> 출력결과 유의하게 나타나는 것은 3,5,번째 약의 형태의 추정치들이 첫 번째 약의 형태와의 차이를 타나내는데 p-값이 유의수준 5% 하에서 유의하다. 항신진대사 점수(x1)는 p 값이 <0.0001로 유의수준 5%하에서 매우 유의하지만 약의 배출 소요시간(x2)은 p 값이 0.131868로 유위수준 5%하에서 유의하지 않은 것으로 나타났다. 위의 출력 결과를 바탕으로 약의 형태별 공분산모형식을 쓰면 다음과 같다.

ŷ_1j = β^₀ + α^₁ + β^1x₁₁_j₊β^2x₂₁_j = 36.637 + 0 - 0.761 x₁₁_j + 0.165x₂₁_j

ŷ_2j = β^₀ + α^₂ + β^1x₁₂_j₊β^2x₂₁_j = 36.637 + 0.897- 0.761 x₁₂_j + 0.165x₂₂_j

ŷ_3j = β^₀ + α^₃ + β^1x₁₃_j₊β^2x₂₁_j = 36.637 + 7.908 - 0.761 x₁₃_j + 0.165x₂₃_j

ŷ_4j = β^₀ + α^₄ + β^1x₁₄_j₊β^2x₂₁_j = 36.637 + 3.072 - 0.761 x₁₄_j + 0.165x₂₄_j

ŷ_5j = β^₀ + α^₅ + β^1x₁₅_j₊β^2x₂₁_j = 36.637 +9.544 - 0.761 x₁₅_j + 0.165x₂₅_j

ŷ6_j = β^₀ + α^₆ + β^1x₁₆_j₊β^2x₂₆_j = 36.637 + 5.839 - 0.761 x₁₆_j + 0.165x₂₆_j

위의 식에 나타난 각 회귀식의 추정결과 모든 약의 형태에 대해서 항신진대사 점수(x1)와 약의 배출 소요시간(x2)의 회귀계수는 동일하다. 약의 배출 소요시간을 제어한 항신진대사 점수의 효과를 보면 항신진대사 점수가 1단위 높아짐에 따라 약의 배출량은 0.761만큼 감소하고 통계적으로 유의한 효과가 있고, 항신진대사 점수를 제어한 약의 배출 소요시간은 1단위 증가할수록 0.165만큼 증가하지만 그 효과는 유의하지 않다.

6. LSMEANS(Adjusted means) 계산

> lsmeans(model3,~type)

type lsmean SE df lower.CL upper.CL

1 5.674385 1.991087 28 1.595828 9.752942

2 6.570912 1.946518 28 2.583650 10.558174

3 13.584122 1.969292 28 9.550210 17.618034

4 8.746633 1.955440 28 4.741095 12.752171

5 15.217831 1.987611 28 11.146394 19.289268

6 11.513284 1.980761 28 7.455879 15.570690

Confidence level used: 0.95

-> 각 약의 형태에 대한 LSMEAN(보정된 평균)값이 출력되어 있다. 보정된 평균은 항신진대사 점수x1와 소요시간x2의 효과를 제어했을 때 반응변수Y인 배출되는 약물량의 평균값이다. 앞의 식에서 각 약의 형태의 회귀식의 공변량값에 전체 평균값을 넣어서 약의 형태별 보정된 평균(adjusted mean)을 아래 식과 같이 계산할 수 있다.

y bar_ad1 = 36.387 + 0 – 0.761*51.22 + 0.165*48.56 = 5.67

y bar_ad2 = 36.387 + 0.897 – 0.761*51.22 + 0.165*48.56 = 6.57

y bar_ad3 = 36.387 + 7.908 – 0.761*51.22 + 0.165*48.56 = 13.59

y bar_ad4 = 36.387 + 3.072 – 0.761*51.22 + 0.165*48.56 = 8.75

y bar_ad5 = 36.387 + 9.544 – 0.761*51.22 + 0.165*48.56 = 15.22

y bar_ad6 = 36.387 + 5.839 – 0.761*51.22 + 0.165*48.56 = 11.51

보정된 평균값을 보면 첫 번째와 두 번째 약의 형태의 보정된 평균값이 다소 작고 세 번째, 다섯 번째, 여섯 번째 약의 형태에 대한 보정된 평균값이 크다는 것을 알 수 있다.

6. 처리 간 다중비교

> model3.lsm = lsmeans(model3,pairwise ~ type, glhargs = list())

> print(model3.lsm,omit=1)

$lsmeans

type lsmean SE df lower.CL upper.CL

1 5.674385 1.991087 28 1.595828 9.752942

2 6.570912 1.946518 28 2.583650 10.558174

3 13.584122 1.969292 28 9.550210 17.618034

4 8.746633 1.955440 28 4.741095 12.752171

5 15.217831 1.987611 28 11.146394 19.289268

6 11.513284 1.980761 28 7.455879 15.570690

Confidence level used: 0.95

$contrasts

contrast estimate SE df t.ratio p.value

1 - 2 -0.8965269 2.782316 28 -0.322 0.9995

1 - 3 -7.9097368 2.768397 28 -2.857 0.0771

1 - 4 -3.0722481 2.824671 28 -1.088 0.8822

1 - 5 -9.5434459 2.853395 28 -3.345 0.0257

1 - 6 -5.8388992 2.839070 28 -2.057 0.3378

2 - 3 -7.0132099 2.783087 28 -2.520 0.1525

2 - 4 -2.1757212 2.753630 28 -0.790 0.9669

2 - 5 -8.6469191 2.802553 28 -3.085 0.0468

2 - 6 -4.9423723 2.758561 28 -1.792 0.4869

3 - 4 4.8374887 2.801787 28 1.727 0.5266

3 - 5 -1.6337091 2.782316 28 -0.587 0.9911

3 - 6 2.0708376 2.840329 28 0.729 0.9766

4 - 5 -6.4711978 2.783344 28 -2.325 0.2178

4 - 6 -2.7666511 2.753890 28 -1.005 0.9125

5 - 6 3.7045468 2.831753 28 1.308 0.7781

P value adjustment: tukey method for comparing a family of 6 estimates

위의 결과를 보면 약의 형태 5의 경우 보정된 평균이 형태 1과 2와 유의하게 차이가 나는 것을 볼 수 있다.

> names(model3.lsm)

[1] "lsmeans" "contrasts"

> plot(model3.lsm[[1]])

> plot(model3.lsm[[2]])

->Tukey(HSD) 검정 결과를 신뢰구간으로 보면, 1-5, 2-5가 유의함을 알 수 있다.

Dunnett vs Tukey

다중분석 하기 전 순서형인 type 변수를 명목형 변수로 변경

> str(drug)

'data.frame': 36 obs. of 5 variables:

$ obs : int 1 2 3 4 5 6 7 8 9 10 ...

$ type: int 1 2 3 4 5 6 1 2 3 4 ...

$ x1 : int 37 37 45 41 57 49 49 53 53 53 ...

$ x2 : int 61 37 53 41 41 33 49 53 45 53 ...

$ y : num 11.32 12.92 18.89 14.67 8.65 ...

> drug$type<-as.factor(drug$type)

> summary(drug)

obs type x1 x2 y

Min. : 1.00 1:6 Min. :37.00 Min. :33.00 Min. : 0.0017

1st Qu.: 9.75 2:6 1st Qu.:49.00 1st Qu.:45.00 1st Qu.: 5.9561

Median :18.50 3:6 Median :53.00 Median :49.00 Median : 8.4527

Mean :18.50 4:6 Mean :51.22 Mean :48.56 Mean :10.2179

3rd Qu.:27.25 5:6 3rd Qu.:53.00 3rd Qu.:53.00 3rd Qu.:13.9175

Max. :36.00 6:6 Max. :61.00 Max. :65.00 Max. :28.1828

> model4<-lm(y~type+x1+x2,data=drug)

다중 분석 시작

> library(multcomp)

> tukey<-glht(model4,linfct=mcp(type='Tukey'))

> summary(tukey)

Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts

Fit: lm(formula = y ~ type + x1 + x2, data = drug)

Linear Hypotheses:

Estimate Std. Error t value Pr(>|t|)

2 - 1 == 0 0.8965 2.7823 0.322 0.9995

3 - 1 == 0 7.9097 2.7684 2.857 0.0771 .

4 - 1 == 0 3.0722 2.8247 1.088 0.8821

5 - 1 == 0 9.5434 2.8534 3.345 0.0257 *

6 - 1 == 0 5.8389 2.8391 2.057 0.3374

3 - 2 == 0 7.0132 2.7831 2.520 0.1524

4 - 2 == 0 2.1757 2.7536 0.790 0.9669

5 - 2 == 0 8.6469 2.8026 3.085 0.0469 *

6 - 2 == 0 4.9424 2.7586 1.792 0.4867

4 - 3 == 0 -4.8375 2.8018 -1.727 0.5265

5 - 3 == 0 1.6337 2.7823 0.587 0.9911

6 - 3 == 0 -2.0708 2.8403 -0.729 0.9766

5 - 4 == 0 6.4712 2.7833 2.325 0.2178

6 - 4 == 0 2.7667 2.7539 1.005 0.9124

6 - 5 == 0 -3.7045 2.8318 -1.308 0.7780

---

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Adjusted p values reported -- single-step method)

> plot(tukey)

->Tukey(HSD) 검정 결과를 신뢰구간으로 보면, 1-5, 2-5가 유의함을 알 수 있다.

> dunnett <- glht(model4,linfct=mcp(type='Dunnett'))

> summary(dunnett)

Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Dunnett Contrasts

Fit: lm(formula = y ~ type + x1 + x2, data = drug)

Linear Hypotheses:

Estimate Std. Error t value Pr(>|t|)

2 - 1 == 0 0.8965 2.7823 0.322 0.9974

3 - 1 == 0 7.9097 2.7684 2.857 0.0324 *

4 - 1 == 0 3.0722 2.8247 1.088 0.7130

5 - 1 == 0 9.5434 2.8534 3.345 0.0103 *

6 - 1 == 0 5.8389 2.8391 2.057 0.1733

---

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Adjusted p values reported -- single-step method)

> plot(dunnett)

7. 모형 적합성 검토

> par(mfrow=c(2,2))

> plot(model3)

-> 등분산성, 정규성 가정에는 큰 문제가 없음을 알 수 있다.

출처: 보건정보데이터 분석(이태림, 이재원, 김주한, 장대흥 공저), R을 이용한 누구나 하는 통계분석(안재형 저)

'KNOU > 2 보건 정보 데이터 분석' 카테고리의 다른 글

제6.2장 비모수적 방법 - 2. 누적한계추정법 (0)	2017.01.26
제6.2장 비모수적 방법 - 1. 생명표 방법 (0)	2017.01.26
제5.1장 공분산 분석 - 1. 공변량이 하나인 경우 (0)	2017.01.06
제4장 범주형 자료의 분석 - 4.3 로짓분석 (0)	2017.01.04
제4장 범주형 자료의 분석 - 4.2.4 대응자료 및 사례 - 대조군 검정 (0)	2016.12.29

Posted by 마르띤

,

제5.1장 공분산 분석 - 1. 공변량이 하나인 경우

KNOU/2 보건 정보 데이터 분석 2017. 1. 6. 14:32

[예] 당뇨병에 걸린 20명의 환자에 대해 혈당을 낮추는 서로 다른 다섯가지 치료법의 효능을 비교하고자 환자 20명을 랜덤하게 5룹으로 나누어 각각의 치료법을 적용하여 한달 후의 혈당량 수치를 측정하였다.

그러나 한 달 후의 혈당량 수치가 초기 혈당량 수치에 영향을 받을 것으로 생각하여 초기 혈당량

수치도함께 측정하였다

관측번호	치료법(trt)	x(초기수치)	y(한달 후 수치)
1	A	27.2	32.6
2	B	22	36.6
…
…

- 독립변수: 치료법(trt)

- 공변량: x(초기 혈당 수치)

- 반응변수: y(한달 후 혈당 수치)

- 반응변수y의 모평균에 영향을 끼칠 수 있는 또 다른 변수가 존재 -> 공변량의 영향을 고려해야 함.

1. 문제제기 - 공분산 분석의 필요성

치료법(trt)에 따른 혈당이 낮아지는 효과를 알아보려 하였으나 초기수치(x)가 낮은 경우 한달 후 수치(y)도 낮아질 것으로 예상된다면, 정확한 치료법(trt)의 효과를 알 수가 없다. 이 때 반응변수 y에 영향을 끼칠 수 있는 또 다른 변수 초기수치(x)에 대해 공변량(Covariate)라 하고, 초기수치(x)를 보정한 상태(adjusted)에서 한달 후 혈당량 수치의(y) 보정된 모평균에 차이가 있는지 보는 것을 공분산분석(Analysis of Covariance: ANCOVA)라 한다. 공분산 분석은 회귀분석과 분산분석의 결합으로 각 처리안에서 공변량을 설명변수 x로 하여 회귀분석을 실시하며, 이렇게 공변량을 고려하면 이를 고려하지 않은 분산분석보다 추정의 정도(precision)을 높일 수 있다.

2. 공분산분석을 위한 두 가지 가정

1) 각 처리 안에서 반응변수Y에 미치는 공변량x의 효과가 모두 동일해야 한다. 즉 교호작용이 없어야 한다.

2) 공변량x 효과가 0이 아니다. 효과가 0이라면 분산분석을 하면 된다.

3. 공분산 모형

yij = β0 + αi + βXij + εij

- Yij : i번째 처리에서 j번째 개체의 반응값

- Xij : i번째 처리에서 j번째 개체의 공변량 값

- αi : 처리의 효과

- β : 모든 처리에 공통으로 작용하는 공변량의 효과

- ε : 등분산을 갖는 정규분포를 따른다고 가정

4. 각 모수의 검정

1) H01: β = 0

- 귀무가설은 처리효과를 제어한 상태에서 반응변수Y에 미치는 공변량효과가 없다는 가정을 검정.

- 만약 처리와 공변량 사이에 교호작용이 존재하면 처리간에 회귀계수가 동일하지 않다는 것을 의미하고, 귀무가설이 기각되지 않으면 분산분석 시행

2) H02 : α1 = α2 = ... = αI

- 공변량 효과를 제어한 상태에서 처리 간 반응변수의 차이가 있는지를 검정

-> 통계 모형을 통해 진짜 알려고 하는 내용

이상 내용에 대해 R 분석을 한 내용은 아래와 같다.

1. library 호출 및 데이터 입력

> library(HH)

> library(lsmeans)

> glucose = read.csv('혈당량자료.csv',header=T)

> head(glucose)

관측번호 치료법 초기혈당량 치료후혈당량

1 1 A 27.2 32.6

2 2 A 22.0 36.6

3 3 A 33.0 37.7

4 4 A 26.8 31.0

5 5 B 28.6 33.8

6 6 B 26.8 31.7

> colname<-c('obs','trt','x','y')

> colnames(glucose)<-colname

> head(glucose)

obs trt x y

1 1 A 27.2 32.6

2 2 A 22.0 36.6

3 3 A 33.0 37.7

4 4 A 26.8 31.0

5 5 B 28.6 33.8

6 6 B 26.8 31.7

> attach(glucose)

2. 회귀계수의 동일성 검정 (교호작용 존재 확인)

공분산분석을 하기 전에 먼저 살펴보아야 할 가정 중의 하나가 바로 처리 간 회귀계수 β의 동일성이다.즉, 처리마다 공변량 효과가 동일해야 함을 의미하는데 만약 처리와 공변량 사이에 교호작용이 존재하면 처리 간에 회귀 계수가 동일하지 않다는 것을 의미하고, 이 경우 공분산 분석을 하는 것은 바람직하지않다.

따라서 먼저 둘 사이에 교호작용이 존재하는가의 여부를 살펴본 후 공분산 분석을 해야 한다.

첫번째 귀무가설은 처리효과를 제어한 반응변수Y에 미치는 공변량 효과를 검정하기 위한 것이다.

귀무가설 H01: β = 0

대립가설 H11: β ≠ 0

두 번째 귀무가설은 공변량 효과를 제어한 상태에서 처리 간 반응변수의 차이가 있는지를 검정.

귀무가설 H02 : α1 = α2 = … = αI

대립가설 H12 : not H02

위 두 가지 귀무가설이 기각되지 않으면 공분산분석을 하지 않고 분산분석을 해야 한다.

> model1 = aov(y~factor(trt)*x,data=glucose)

> summary(model1)

Df Sum Sq Mean Sq F value Pr(>F)

factor(trt) 4 198.41 49.60 15.868 0.000248 ***

x 1 92.53 92.53 29.601 0.000285 ***

factor(trt):x 4 36.48 9.12 2.917 0.077290 .

Residuals 10 31.26 3.13

---

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

-> aov함수: 공분산분석의 가정 중 하나인 처리 간 회귀계수의 동일성을 확인하기 위해 처리와 공변량 사이의 교호작용의 유무를 검정. aov함수의 우변에 factor(trt) * 공변량 X를 입력. 치료법 trt앞에 factor를 입력하는 것은 치료법trt를 의미하는 ABCDE는 각각 명목형 변수이기 때문.

귀무가설 H01: 교호작용이 존재하지 않는다.

대립가설 H11: 교호작용이 존재한다.

p-value : 0.077290

의사결정: 유의수준 5% 하에서 귀무가설을 기각할 수 없다.

결론: 교호작용이 존재하지 않으므로 공분산분석을 할 수 있다. ( = 각 처리 안에서 반응변수Y에 미치는 공변량x의 효과가 모두 동일하다.)

Ancova(HH Library)함수를 이용하면, 치료법 간 공변량의 효과가 동일하다는 가정을 할 수 있는 xyplot을 그릴 수 있다.

> ancova(y~trt*x,data=glucose)

Analysis of Variance Table

Response: y

Df Sum Sq Mean Sq F value Pr(>F)

trt 4 198.407 49.602 15.8683 0.0002484 ***

x 1 92.528 92.528 29.6012 0.0002846 ***

trt:x 4 36.476 9.119 2.9173 0.0772904 .

Residuals 10 31.258 3.126

---

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1

-> HH라이브러리 ancova 함수는 aov함수와 동일한 기능, trellis 그림인 xyplot시각화 기능 제공.함수 출력 결과는 위 aov함수와 동일하며, F값은 2.92, p-value는 0.0773으로 유의수준 5% 하에서 유의하지 않으므로 귀무가설을 기각할 수 없다.따라서 공변량과 처리 사이에 교호작용이 존재하지 않으므로 공분산분석을 할 수 있다는 것을 알 수 있다.

3. 일원 공분산분석 (One-way ANCOVA) - 공변량 효과 제어 시 치료법의 효과 검정

처리 간 공변량의 효과가 동일하다는 가정을 확인한 후 공분산 분석을 출력

> model2 = lm(y~factor(trt)+x,data=glucose)

> summary(model2)

Call:

lm(formula = y ~ factor(trt) + x, data = glucose)

Residuals:

Min 1Q Median 3Q Max

-3.1360 -1.0024 -0.2827 0.7257 6.0806

Coefficients:

Estimate Std. Error t value Pr(>|t|)

(Intercept) 13.9437 4.8219 2.892 0.011834 *

factor(trt)B -2.7685 1.5554 -1.780 0.096793 .

factor(trt)C -1.6660 1.6186 -1.029 0.320776

factor(trt)D -1.6284 1.5618 -1.043 0.314787

factor(trt)E -4.5903 1.9115 -2.401 0.030788 *

x 0.7534 0.1723 4.373 0.000637 ***

---

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.2 on 14 degrees of freedom

Multiple R-squared: 0.8112, Adjusted R-squared: 0.7437

F-statistic: 12.03 on 5 and 14 DF, p-value: 0.0001164

-> lm 함수: 공분산분석에는 lm함수를 사용. ~를 중심으로 좌변에는 반응변수y, 우변에는 치료법trt와 공변량 x를 입력하고 그 사이는 *가 아닌 +를 사용.

귀무가설 H01: β = 0

H02 : α1 = α2 = … = αI

대립가설 H11: β ≠ 0

H12 : not H02

p-value : 0.0001164

의사결정: 유의수준 5% 하에서 귀무가설을 매우 강하게 기각할 수 있다.

결론: 우리가 세운 모형의 자료에 적합하다는 것을 알 수 있다. 낮아진 혈당량의 모평균이 처trt의 효과와 공변량들의 효과가 없다라고 말할 수 없다.유의수준 5%하에서 혈당량의 초기수치(x)가 즉 모든 공변량이 0이라는 귀무가설을 기각할만한 증거가 충분하다.

-> 결과 해석:

- 치료법A를 0으로 둘 때 치료법E p-value가 0.030788으로 유의하게 차이가 난다는 것을 알 수 있다.

1종 제곱합과 회귀계수 동일성 확인 위한 xyplot 그래프는 아래와 같다.

> ancova(y~trt+x,data=glucose)

Analysis of Variance Table

Response: y

Df Sum Sq Mean Sq F value Pr(>F)

trt 4 198.407 49.602 10.252 0.0004301 ***

x 1 92.528 92.528 19.125 0.0006369 ***

Residuals 14 67.734 4.838

---

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’

-> ancova(y~trt+x,data=glucose), 공분산분석이므로 *가 아닌 +를 입력.

-> 그래프 결과해석: 5개 그래프의 적합회귀식 절편은 다르지만 기울기는 동일하므로 공변량 효과가 처리마다 다르지 않음을 알 수 있다. “2. 공분산분석을 위한 두 가지 가정” 중 첫 번째 내용 “각 처리 안에서 반응변수Y에 미치는 공변량x의 효과가 모두 동일해야 한다. 즉 교호작용이 없어야 한다.”을 만족시킨다.

4. 공변량 효과 제어시 치료법의 효과 검정 - 모형 제곱합

1) 제1종 제곱합 Type I SS

SS(trt,X) = SS(trt) + SS(X | trt)

처리가 기여한 부분 + 처리의 기여 후 순수 공변량 x이 기여한 부분

2) 제3종 제곱합 Type III SS

SS(trt | X) + SS(X | trt)

공변량x가 기여한 상태에서 처리가 기여한 부분 + 처리의 기여 후 순수 공변량 x이 기여한 부분

x given trt, 초기혈당수치 x가 고려된 후 치료법간 차이trt의 차이를 확인한다는 분석의 목적. 즉 x 기여 외 순수 trt 기여를 아는 것이 목표이기 때문에, 우리가 관심가지는 분야 역시 제3종 제곱함 SS(trt | X) 부분

제1종 제곱합

> ancova(y~trt+x,data=glucose)

Analysis of Variance Table

Response: y

Df Sum Sq Mean Sq F value Pr(>F)

trt 4 198.407 49.602 10.252 0.0004301 ***

x 1 92.528 92.528 19.125 0.0006369 ***

Residuals 14 67.734 4.838

---

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’

제3종 제곱합

> summary(aov(y~x+factor(trt),data=glucose))

Df Sum Sq Mean Sq F value Pr(>F)

x 1 256.75 256.75 53.067 4e-06 ***

factor(trt) 4 34.19 8.55 1.767 0.192

Residuals 14 67.73 4.84

---

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1

제1종 제곱합은

- F-값10.252

- p-value: 0.0004301 < 0.05

- 의사결정: ss(trt) 제1종 제곱합은 유의하게 나와 치료법이 5% 유의수준 하에서 유의하다.

제3종 제곱합

- F-값1.767

- p-value: 0.192 > 0.05

- 의사결정: ss(trt | x) x given trt, 공변량x이 기여한 상태에서 처리trt가 기여한 순수한 부분을 확인하는 제3종 제곱합은 치료법trt이 5% 유의수준 하에서 유의하지 않다. 공분산분석은 공변량x의 효과를 제어했을 때 치료법trt에 따라 혈압의 보정평균이 차이가 나는지가 관심사이기 때문에 치료법의 유의성 결과는 제3종 제곱합에 나타난 결과를 보아야 한다.

초기수치

초기수치에 대해서는 제1종 제곱합과 제3종 제곱합이 같게 나오는데 이것은 제1종 제곱합에서는 SS(X | trt)이고 고유기여분을 나타내는 제3종 제곱합에서도 SS(X | trt)이기 때문이다. 초기 수치에 대한 p-value는 0.0006으로 유의수준 5%하에서 매우 유의함을 알 수 있다.

귀무가설 H02 : α1 = α2 = … = αI

p-value: 0.0006

의사결정: 초기수치의 효과가 0이라는 귀무가설을 기각하게 되어 앞의 공분산분석을 하기 위해 만족해야 하는 첫 번째 가정이 충족되는 것을 알 수 있다.

4. 각 요인의 수준별 추정치 분석

> model2 = lm(y~factor(trt)+x,data=glucose)

> summary(model2)

Call:

lm(formula = y ~ factor(trt) + x, data = glucose)

Residuals:

Min 1Q Median 3Q Max

-3.1360 -1.0024 -0.2827 0.7257 6.0806

Coefficients:

Estimate Std. Error t value Pr(>|t|)

(Intercept) 13.9437 4.8219 2.892 0.011834 *

factor(trt)B -2.7685 1.5554 -1.780 0.096793 .

factor(trt)C -1.6660 1.6186 -1.029 0.320776

factor(trt)D -1.6284 1.5618 -1.043 0.314787

factor(trt)E -4.5903 1.9115 -2.401 0.030788 *

x 0.7534 0.1723 4.373 0.000637 ***

---

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.2 on 14 degrees of freedom

Multiple R-squared: 0.8112, Adjusted R-squared: 0.7437

F-statistic: 12.03 on 5 and 14 DF, p-value: 0.0001164

위의 각 추정 계수를 바탕으로 공분산 모형식을 쓰면 아래 식과 같다.

ŷ_1j = β^₀ + α^₁ + β^x_1j = 13.9437 + 0 + 0.7534x_1j

ŷ_2j = β^₀ + α^₂ + β^x_2j = 13.9437 -2.7685 + 0.7534x_2j

ŷ_3j = β^₀ + α^₃ + β^x_3j = 13.9437 -1.6660 + 0.7534x_3j

ŷ_4j = β^₀ + α^₄ + β^x_4j = 13.9437 -1.6284 + 0.7534x_4j

ŷ_5j = β^₀ + α^₅ + β^x_5j = 13.9437 -4.5903 + 0.7534x_5j

각 처리의 회귀식을 살펴보면 공변량 초기혈당량의 회귀계수 0.7534는 모두 처리에 대해 공통이고 반응변수y의 절편에서만 값이 차이가 난다는 것을 알 수 있다. 첫 번째 치료법 A의 처리효과를 0으로놓았기 때문에 Inetercept 부분의 회귀계수 13.9437은 첫번째 치료법 A의 상수항임을 알 수 있고, 치료법B부터 치료법 E에 대한 추정치는 각각 첫 번째

치료법 A의 추정치와의 차이가 난다.

치료법A(trt A)와 치료법 E(trt E)의 차이에 대한 p-value가 0.0308로 유의하게 나와 치료법 A가 치료법E와 유의하게 차이가 난다는 것을 알 수 있다.

혈당량의 초기 수치(x)에 대한 모수 추정치는 0.753이고 t-값과 p-값이 각각 4.373과 0.000637으로유의수준 5%하에서 혈당량의 초기 수치(x)가 0이라는 귀무가설을 기각할 만한 증거가 충분하다.

즉 공변량 효과가 있다. 이는 공분산결과 분석을 해도 된다.

5. LSMEAN(adjusted mean) 분석

> lsmeans(model2, ~ trt)

trt lsmean SE df lower.CL upper.CL

A 32.97565 1.151991 14 30.50487 35.44642

B 30.20716 1.148211 14 27.74449 32.66982

C 31.30960 1.104799 14 28.94004 33.67916

D 31.34724 1.117954 14 28.94946 33.74501

E 28.38536 1.341631 14 25.50785 31.26287

Confidence level used: 0.95

Warning message:

In model.frame.default(trms, grid, na.action = na.pass, xlev = xlev) :

variable 'trt' is not a factor

보정된(adjust) 평균은 각 처리마다 추정된 회귀식에서 공변량값에 각 처리 평균 대신 공변량의 전체평균을

사용하였을 때 기대되는 반응변수의 평균이다. 즉, 공변량 효과를 보정한 상태에서의 반응변수의 평균으로, 모든 처리에서 공변량 평균이 같다고 했을 때의 평균이다. 공분산분석에서는 공변량 효과를

보정한 보정 평균 간에 차이가 있는지가 가장 큰 관심사이다. 앞에서 구한 각 처리 회귀식의 공변량에

공변량의 전체 평균값을 넣어서 처리별 보정된 평균을 다음 식과 같이 계산할 수 있다.

y bar_ad1 = β^₀ + α^₁ + β^xbar = 13.9437 + 0 + 0.7534 X 25.26 = 32.975

y bar_ad2 = β^₀ + α^₂ + β^xbar = 13.9437 -2.7685 + 0.7534 X 25.26 =30.207

y bar_ad3 = β^₀ + α^₃ + β^xbar = 13.9437 -1.6660 + 0.7534 X 25.26 =31.309

y bar_ad4 = β^₀ + α^₄ + β^xbar = 13.9437 -1.6284 + 0.7534 X 25.26 =31.347

y bar_ad5 = β^₀ + α^₅ + β^xbar = 13.9437 -4.5903 + 0.7534 X 25.26 =28.385

각 처리마다 계산된 보정된 평균값의 차이는 공변량의 값이 전체 평균으로 같기 때문에 평균값의 차이는 절편의 차이가 됨을 알 수 있다. 각 처리마다 추정된 회귀계수 값들은 첫 번째 처리(trt A)와 추정치 차이인 동시에 첫 번째 처리와의 보정된 평균 차이이기도 하다.

6. Tukey(HSD) 검정 (처리 간 다중 비교)

> model2.lsm = lsmeans(model2,pairwise ~ trt,glhargs=list())

> print(model2.lsm, omit = 2)

$lsmeans

trt lsmean SE df lower.CL upper.CL

A 32.97565 1.151991 14 30.50487 35.44642

B 30.20716 1.148211 14 27.74449 32.66982

C 31.30960 1.104799 14 28.94004 33.67916

D 31.34724 1.117954 14 28.94946 33.74501

E 28.38536 1.341631 14 25.50785 31.26287

Confidence level used: 0.95

$contrasts

contrast estimate SE df t.ratio p.value

A - B 2.76849175 1.555390 14 1.780 0.4216

A - C 1.66604727 1.618557 14 1.029 0.8378

A - D 1.62840923 1.561818 14 1.043 0.8316

A - E 4.59029035 1.911531 14 2.401 0.1718

B - C -1.10244448 1.615029 14 -0.683 0.9570

B - D -1.14008252 1.560695 14 -0.730 0.9457

B - E 1.82179860 1.904048 14 0.957 0.8696

C - D -0.03763804 1.585115 14 -0.024 1.0000

C - E 2.92424307 1.690871 14 1.729 0.4485

D - E 2.96188112 1.832554 14 1.616 0.5115

P value adjustment: tukey method for comparing a family of 5 estimates

> plot(model2.lsm[[2]])

앞서 본 모델에서는 A와 E 치료법 사이에 유의한 차이가 있었지만, LSD 검정(Tukey(HSD)) 검정에서는어느 치료법도 혈당 수치에 있어 유의한 차이가 없다. 왜일까? Tukey검정방법은 상대적으로 보수적이기 때문에 정말 큰 차이가 날 경우에만 유의한 차이를 보인다.

Dunnett vs Tukey

모든 평균이 같다는 귀무가설이 기각되었다는 말은 그룹 중 최소한 하나는 0이 아니다라는 말이다. 어느 쌍의 차이로 귀무가설이 기각되었는지 조사하기 위해 다중비교를 한다. 분산분석에서 많이 쓰이는 다중비교 방법은 Dunnett와 Tukey이다. Tukey는 가능한 모든 조합의 쌍을, Dunnett는 하나의 대조군(reference)을 나머지 비교군(treatment)들과 비교한다.

6-1) Tukey

> library(multcomp)

> model3 = lm(y~trt+x,data=glucose)

> tukey = glht(model3,linfct=mcp(trt='Tukey'))

> summary(tukey)

Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts

Fit: lm(formula = y ~ trt + x, data = glucose)

Linear Hypotheses:

Estimate Std. Error t value Pr(>|t|)

B - A == 0 -2.76849 1.55539 -1.780 0.419

C - A == 0 -1.66605 1.61856 -1.029 0.836

D - A == 0 -1.62841 1.56182 -1.043 0.830

E - A == 0 -4.59029 1.91153 -2.401 0.170

C - B == 0 1.10244 1.61503 0.683 0.956

D - B == 0 1.14008 1.56069 0.730 0.945

E - B == 0 -1.82180 1.90405 -0.957 0.868

D - C == 0 0.03764 1.58512 0.024 1.000

E - C == 0 -2.92424 1.69087 -1.729 0.446

E - D == 0 -2.96188 1.83255 -1.616 0.509

(Adjusted p values reported -- single-step method)

결과해석: 모든 치료방법 간 차이의 p-value가 0.05보다 크기 때문에 유의하지 않다.

> plot(tukey)

결과 해석: 치료 방법의 차이 신뢰구간이 0을 포함하고 있으므로 서로 유의하지 않다. 즉, 어느 치료법도 혈당수치에 있어 유의한 차이가 없다.

6-2) Dunnett

> dunnett=glht(model3,linfct=mcp(trt='Dunnett'))

> summary(dunnett)

Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Dunnett Contrasts

Fit: lm(formula = y ~ trt + x, data = glucose)

Linear Hypotheses:

Estimate Std. Error t value Pr(>|t|)

B - A == 0 -2.768 1.555 -1.780 0.2700

C - A == 0 -1.666 1.619 -1.029 0.7006

D - A == 0 -1.628 1.562 -1.043 0.6919

E - A == 0 -4.590 1.912 -2.401 0.0953 .

---

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Adjusted p values reported -- single-step method)

결과해석: 모든 치료방법 간 차이의 p-value가 0.05보다 크기 때문에 유의하지 않다.

> plot(dunnett)

결과 해석: 치료 방법의 차이 신뢰구간이 0을 포함하고 있으므로 서로 유의하지 않다. 즉, 어느 치료법도 혈당수치에 있어 유의한 차이가 없다.

7. 잔차 분석

> par(mfrow=c(2,2))

> plot(모형2)

특이값 2번이 존재하지만, 정규분포를 따르고 등분산 가정은 큰 문제는 없다.

출처: 보건정보데이터 분석(이태림, 이재원, 김주한, 장대흥 공저), R을 이용한 누구나 하는 통계분석(안재형 저)

'KNOU > 2 보건 정보 데이터 분석' 카테고리의 다른 글

제6.2장 비모수적 방법 - 1. 생명표 방법 (0)	2017.01.26
제5.1장 공분산 분석 - 2 공변량이 둘 이상인 경우 (2)	2017.01.13
제4장 범주형 자료의 분석 - 4.3 로짓분석 (0)	2017.01.04
제4장 범주형 자료의 분석 - 4.2.4 대응자료 및 사례 - 대조군 검정 (0)	2016.12.29
제2장 보건정보 데이터의 기초분석 (1)	2016.12.27

Posted by 마르띤

,

제4장 범주형 자료의 분석 - 4.3 로짓분석

KNOU/2 보건 정보 데이터 분석 2017. 1. 4. 10:03

로짓분석이란?

생존여부나 교통사고 발성 여부와 같이 반응변수(Y)는 범주형(명목형(남자 1, 여자2), 순서형(초졸1, 중졸2, 고졸3, 대졸4))이고 이에 대한 설명변수(X)는 범주형과 이산형(방의 개수: 2개)dl 혼합된 경우 관련성 여부를 규명하기 위해서는 모형으로 로짓logit 모형을 적용할 수 있다.

오즈: 기본이 되는 변수로 오즈odds를 들 수 있는데 이것은 두 확률의 비를 의미하는 것이다.

odds=p/(1-p)

로짓: 흡연 산모가 미숙아를 출산할 확률이 비흡연 산모를 1로 했을 때 2.5배라고 발표하는 것이 그 예이다. 로짓은 오즈에 자연로그를 취한 형태를 의미한다.

logit=ln(p/1-p)

로짓 모형: 범주형 자료 분석에서 실험자들의 설명변수(X)에 대한 통제가 가능한 경우 반응변수(Y)에 대한 로짓모형을 적용할 수 있다.

logit=ln(p/1-p) = β₀ + β₁x

예) B.J.T Morgan은 폴란드 바르샤바의 3918명 여성들을 대상으로 초사한 초경자료를 가지고 범주형 자료의 회귀분석(을 소개하였다). 연령과 월경 사이에 어떠한 관계가 있는지 알아보자.

> 초경연령자료 = read.table('c:/Rwork/바르샤바초경연령자료.csv',sep=',',header=T)

> head(초경연령자료)

평균연령 초경경험자 군의총수

1 9.21 0 376

2 10.21 0 200

3 10.58 0 93

4 10.83 2 120

5 11.08 2 90

6 11.33 5 88

> attach(초경연령자료)

> 확률 = 초경경험자/군의총수

> plot(평균연령,확률)

#로짓분석

> logit <- glm(확률~평균연령,data=초경연령자료,family='binomial')

> logit

Call: glm(formula = 확률 ~ 평균연령, family = "binomial", data = 초경연령자료)

Coefficients:

(Intercept) 평균연령

-20.907 1.608

Degrees of Freedom: 23 Total (i.e. Null); 22 Residual

Null Deviance: 19.06

Residual Deviance: 0.2214 AIC: 11.38

> summary(logit)

Call:

glm(formula = 확률 ~ 평균연령, family = "binomial", data = 초경연령자료)

Deviance Residuals:

Min 1Q Median 3Q Max

-0.20043 -0.08458 -0.05277 0.06104 0.13296

Coefficients:

Estimate Std. Error z value Pr(>|z|)

(Intercept) -20.9067 8.1314 -2.571 0.0101 *

평균연령 1.6077 0.6244 2.575 0.0100 *

---

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

Null deviance: 19.06046 on 23 degrees of freedom

Residual deviance: 0.22137 on 22 degrees of freedom

AIC: 11.383

Number of Fisher Scoring iterations: 6

결과: 로짓모형을 추정하면 아래와 같다. 즉, X가 1단위 증가할 때 변화하는 로그 오즈(odds)의 비율은 1.6077이다.

logit(p) = -20.9067 + 1.6077 X

> exp(coef(logit))

(Intercept) 평균연령

8.323842e-10 4.991105e+00

결과해석: 나이가 한 살 증가할 때 초경 경험자의 오즈 odds는 4.991 (exp(1.6077))배 더 많다고 할 수 있다.

그룹의 50%가 초경을 하는 연령을 구하기 위해 p=0.5일 때의 연령 X를 구하면

logit(0.5) =	log(	0.5	) =	log1 = 0
		1-0.5

이고 이것이 바로 유효 중앙값 ED₅₀(Effective Dose of 50%)이 되어 50%가 초경을 경험한 소녀들의 나이가 된다. 이를 계산하면 아래 수식을 통해 x=13, 즉 약 13세임을 알 수 있다.

0 = -20.9067 + 1.6077 X

이식은 양의 용량이 늘어날수록 반응률이 높아지는 반응곡선, 즉 용량-반응곡선(dose response curve)에도 응용된다. 예를 들면 100명의 고열 환자에게 해열제를 투여했을 때 환자 50%의 열이 떨어지는 효과를 보여주는 약의 용량인 ED₅₀을 구하는데 적용할 수 있다.

예) M.J.R Healy는 비타민 E의 용량에 따른 임신한 쥐의 숫자를 아래와 같이 발표하였다. 투입된 비타민 E의 용량과 임신에 대한 관계를 알아보자.

용량(mg)	수	임신
3.75	5	0
5	10	2
6.25	10	4
7.5	10	8
10	11	10
15	11	11

> 용량 = c(3.75,5,6.25,7.5,10,15)

> 수= c(5,10,10,10,11,11)

> 임신 = c(0,2,4,8,10,11)

> 쥐임신<-cbind(용량,수,임신)

> 쥐임신

용량 수 임신

[1,] 3.75 5 0

[2,] 5.00 10 2

[3,] 6.25 10 4

[4,] 7.50 10 8

[5,] 10.00 11 10

[6,] 15.00 11 11

> 확률2 = 쥐임신[,3] / 쥐임신[,2]

> 확률2

[1] 0.0000000 0.2000000 0.4000000 0.8000000 0.9090909 1.0000000

#로짓분석

> logit2 <- glm(확률2~log10(쥐임신[,1]),family='binomial')

> logit2

Call: glm(formula = 확률2 ~ log10(쥐임신[, 1]), family = "binomial")

Coefficients:

(Intercept) log10(쥐임신[, 1])

-12.42 15.35

Degrees of Freedom: 5 Total (i.e. Null); 4 Residual

Null Deviance: 4.297

Residual Deviance: 0.136 AIC: 6.315

> summary(logit2)

Call:

glm(formula = 확률2 ~ log10(쥐임신[, 1]), family = "binomial")

Deviance Residuals:

1 2 3 4 5 6

-0.23099 0.11888 -0.09852 0.15618 -0.16518 0.08464

Coefficients:

Estimate Std. Error z value Pr(>|z|)

(Intercept) -12.42 10.39 -1.195 0.232

log10(쥐임신[, 1]) 15.35 12.76 1.203 0.229

(Dispersion parameter for binomial family taken to be 1)

Null deviance: 4.29706 on 5 degrees of freedom

Residual deviance: 0.13603 on 4 degrees of freedom

AIC: 6.3152

Number of Fisher Scoring iterations: 6

결과: glm(formula = 확률2 ~ log10(쥐임신[, 1]), family = "binomial")를 보면 모형을 적합하는데 용량(X) 대신 log ₁₀(X)를 사용하였다. 위의 결과를 통해 다음 추정식을 회귀할 수 있다.

ln(	p_x	) =	-12.42 + 15.35 log₁₀x
	1 - p_x

따라서 비타민 E 용량이 10배 증가함에 따라 변화하는 임신 로그 오즈의 비율은15.35임을 알 수 있다.

> exp(coef(logit2))

(Intercept) log10(쥐임신[, 1])

4.033060e-06 4.631322e+06

비타민 용량 E가 10배 증가할 때 임신의 오즈는 4.631(=exp(15.35))배 더 많다고 할 수 있다.

전체 그룹 중 50%가 임신하는 비타민 용량 E는 아래와 같이 구할 수 있다.

logit(0.5) =	log(	0.5	) =	log1 = 0
		1-0.5

수식을 이용하면 ED₅₀(Effective Dose of 50%)은 0 = -12.42 + 15.35 log₁₀x이므로 ED₅₀(Effective Dose of 50%)의 추정량은 6.44mg이 된다.

> 10^(12.42/15.35)

[1] 6.443481

출처: 보건정보데이터분석(이태림, 이재원, 김주한, 장대흥 공저)

'KNOU > 2 보건 정보 데이터 분석' 카테고리의 다른 글

제5.1장 공분산 분석 - 2 공변량이 둘 이상인 경우 (2)	2017.01.13
제5.1장 공분산 분석 - 1. 공변량이 하나인 경우 (0)	2017.01.06
제4장 범주형 자료의 분석 - 4.2.4 대응자료 및 사례 - 대조군 검정 (0)	2016.12.29
제2장 보건정보 데이터의 기초분석 (1)	2016.12.27
제4장 범주형 자료의 분석 - 4.2.2 독립성 검정 (카이제곱 검정) (0)	2016.11.08

Posted by 마르띤

,

제4장 범주형 자료의 분석 - 4.2.4 대응자료 및 사례 - 대조군 검정

KNOU/2 보건 정보 데이터 분석 2016. 12. 29. 10:09

4.2.4 대응자료 및 사례 – 대조군 검정

유권자	이전	이후
1	1	0
2	0	0
3	1	1
. . .	. . .	. . .
100	0	1

선거유세의 효과를 파악하는 방법으로 유권자 100명을 확률표본으로 추출하여 유세 이전의 지지(1=여당, 0 = 야당)와 유세 이후의 지지(1=여당, 0=야당) 데이터를 만들었다. 유세 이전과 이후 유권자의 지지변화가 있었는지 알아보자

데이터를 table로 만들면 아래와 같다.

구분		유세이후
구분		X2=1	X2=0	합
유세이전	X1=1	63	4	67
	X1=0	21	12	33
	합	84	16	100

> 선거유세효과<-matrix(c(63,4,21,12),2,2,byrow=T,dimnames=list(유세이전=c('유세전 여당','유세전 야당'),유세이후=c('유세후 여당','유세후 야당')))

> 선거유세효과

유세이후

유세이전 유세후 여당 유세후 야당

유세전 여당 63 4

유세전 야당 21 12

> mcnemar.test(선거유세효과,correct=F) #McNemar 검정, {(21-4)/(21+4)}^2=11.56, 11.56보다 같거나 더 극단적인 값이 나올 확률 0.0006739

McNemar's Chi-squared test

data: 선거유세효과

McNemar's chi-squared = 11.56, df = 1, p-value = 0.0006739

> mcnemar.test(선거유세효과) #연속성 수정

McNemar's Chi-squared test with continuity correction

data: 선거유세효과

McNemar's chi-squared = 10.24, df = 1, p-value = 0.001374

> library(exact2x2)

> exact2x2(선거유세효과,paired=T)

Exact McNemar test (with central confidence intervals)

data: 선거유세효과

b = 4, c = 21, p-value = 0.0009105

alternative hypothesis: true odds ratio is not equal to 1

95 percent confidence interval:

0.04753664 0.56452522

sample estimates:

odds ratio

0.1904762

귀무가설: H0 유세 전과 후의 지지율이 같다

대립가설: H1 유세 전과 후의 지지율이 다르다

검정 결과: McNemar 검정 결과 p-값 0.0009105

결론: 귀무가설을 기각, 유세 이전과 이후에 유권자의 지지에 변화가 일어났다.

McNemar 검정은 변화의 유의성 검정 이외에 역학 연구에서 흔히 이루어지는 사례-대조군 연구에도 사용된다. 역학연구에서는 후향연구의 한 형태로 사례-대조군 연구를 실시하고 경우에 따라서 사례에 대한 대조를 1:1 대응으로 찾는 경우가 있다. 다음의 호지킨병과 편도적출(tonsillectomy)자료는 1:1 대응의 사례-대조군 연구 자료이고 McNemar 검정의 좋은 예이다.

[예] 1972년에 편도적출과 호지킨병의 관계를 규명하는 연구보고가 있었다. 미국 국립 암 연구소에서 치료를 받은 174명의 호지킨병 환자들을 사례군으로 하고 환자들의 형제자매 472명을 대조군으로 구성하였다. 그리고 이 자료를 기초로 한 추후의 연구보고에서는 사례와 대조를 1:1 대응으로 하기 위해서 사례 한 명에 대응되는 대조는 나이 차이가 5년 이내이고 같은 성별을 가진 형제자매 중 사례와 나이가 가장 가까운 사람으로 정하였다. 이렇게 하여 85쌍의 사례 – 대조군 자료가 관찰되었다. 호지킨병과 편도적출이 독립적인가를 검정하는 문제를 생각해보자.

> 호지킨병 = matrix(c(26,15,7,37),nrow=2,byrow=T)

> 호지킨병

[,1] [,2]

[1,] 26 15

[2,] 7 37

> dimnames(호지킨병) = list(사례군=c('편도추출 유','편도추출 무'), 대조군=c('편도추출 유','편도추출 무'))

> 호지킨병

대조군

사례군 편도추출 유 편도추출 무

편도추출 유 26 15

편도추출 무 7 37

> 분할표_호지킨병 = addmargins(호지킨병)

> 분할표_호지킨병

대조군

사례군 편도추출 유 편도추출 무 Sum

편도추출 유 26 15 41

편도추출 무 7 37 44

Sum 33 52 85

> mcnemar.test(호지킨병,correct=F) #McNemar 검정

McNemar's Chi-squared test

data: 호지킨병

McNemar's chi-squared = 2.9091, df = 1, p-value = 0.08808

> mcnemar.test(호지킨병) #McNemar 검정(연속성 수정)

McNemar's Chi-squared test with continuity correction

data: 호지킨병

McNemar's chi-squared = 2.2273, df = 1, p-value = 0.1356

> library(exact2x2)

> exact2x2(호지킨병,paired=T)

Exact McNemar test (with central confidence intervals)

data: 호지킨병

b = 15, c = 7, p-value = 0.1338

alternative hypothesis: true odds ratio is not equal to 1

95 percent confidence interval:

0.8224084 6.2125863

sample estimates:

odds ratio

2.142857

귀무가설: H0 호지킨병과 편도적출이 독립적이다

대립가설: H1 호지킨병과 편도적출이 독립적이지 않다

검정결과: McNemar 검정결과 p-value 0.1338

결론: a=0.05에서 H0를 기각시키기에 충분한 증거를 제시하지 않는다고 결론

출처: 보건정보데이터 분석(이태림, 이재원, 김주한, 장대흥 공저)

'KNOU > 2 보건 정보 데이터 분석' 카테고리의 다른 글

제5.1장 공분산 분석 - 1. 공변량이 하나인 경우 (0)	2017.01.06
제4장 범주형 자료의 분석 - 4.3 로짓분석 (0)	2017.01.04
제2장 보건정보 데이터의 기초분석 (1)	2016.12.27
제4장 범주형 자료의 분석 - 4.2.2 독립성 검정 (카이제곱 검정) (0)	2016.11.08
제4장 범주형 자료의 분석 - 4.2 범주형 자료의 검정(카이제곱 검정) (0)	2016.11.04

Posted by 마르띤

,

제2장 보건정보 데이터의 기초분석

KNOU/2 보건 정보 데이터 분석 2016. 12. 27. 10:11

p.46 제2장 보건 정보 데이터의 기초 분석

2.3 자료의 기술 및 요약

> setwd('c:/Rwork/')

> 담즙과포화비율자료 = read.table('담즙과포화비율.txt',header=T)

> attach(담즙과포화비율자료)

> head(담즙과포화비율자료,2)

성별 담즙의과포화비율

1 남자 40

2 남자 88

> plot(담즙의과포화비율,type='p',xlab='자료',ylab='담즙과포화비율',main='담즙과포화비율')

> par(new=T)

> plot(담즙의과포화비율,type='h',xlab='자료',ylab='담즙과포화비율',main='담즙과포화비율')

> length(담즙의과포화비율)

[1] 60

> length(담즙의과포화비율) #길이

[1] 60

> sort(담즙의과포화비율,decreasing=T) #내림차순

[1] 146 142 137 128 127 123 123 120 118 116 112 111 110 110 107 106 106 98 91 90 89 88 88

[24] 88 87 87 86 86 86 84 84 82 80 80 80 79 78 77 76 76 75 74 73 73 69 67

[47] 66 66 65 65 58 58 57 56 55 52 52 47 40 35

> sum(담즙의과포화비율) #합

[1] 5185

> cumsum(담즙의과포화비율) #누적합

[1] 40 128 238 295 381 518 596 707 795 875 961 1041 1088 1194 1259 1333 1399 1478

[19] 1536 1659 1746 1834 1924 1980 2053 2165 2275 2393 2445 2551 2618 2683 2735 2819 2905 2940

[37] 3056 3132 3187 3260 3349 3476 3563 3705 3782 3858 3916 4007 4114 4212 4340 4424 4570 4645

[55] 4765 4845 4927 5050 5116 5185

> mean(담즙의과포화비율);median(담즙의과포화비율) #평균과 중앙값

[1] 86.41667

[1] 84

> mean(담즙의과포화비율,trim=1/10) #10% 절삭

[1] 85.4375

> var(담즙의과포화비율);sd(담즙의과포화비율) #분산과 표준편차

[1] 657.1624

[1] 25.63518

> fivenum(담즙의과포화비율) #다섯숫자요약

[1] 35.0 68.0 84.0 106.5 146.0

> quantile(담즙의과포화비율)

0% 25% 50% 75% 100%

35.00 68.50 84.00 106.25 146.00

> IQR(담즙의과포화비율) #3사분위수-1사분위수

[1] 37.75

> quantile(담즙의과포화비율)[4]-quantile(담즙의과포화비율)[2]

75%

37.75

> mad(담즙의과포화비율) #Median Absolut Deviation 각 데이터에서 중앙값을 뺀 후 절대값을 취한 값들의 중앙값

[1] 26.6868

> max(담즙의과포화비율);min(담즙의과포화비율)

[1] 146

[1] 35

> range(담즙의과포화비율)

[1] 35 146

> R=max(담즙의과포화비율)-min(담즙의과포화비율)

> R

[1] 111

2. 4 표와 그래프를 이용한 자료의 요약

> head(담즙과포화비율자료,2)

성별 담즙의과포화비율

1 남자 40

2 남자 88

> 계급 = cut(담즙의과포화비율, breaks=c(20,40,60,80,100,120,140,160))

> head(계급)

[1] (20,40] (80,100] (100,120] (40,60] (80,100] (120,140]

Levels: (20,40] (40,60] (60,80] (80,100] (100,120] (120,140] (140,160]

> table(계급)

계급

(20,40] (40,60] (60,80] (80,100] (100,120] (120,140] (140,160]

2 8 18 15 10 5 2

> par(mfrow=c(1,1))

> hist(담즙의과포화비율,breaks=c(20,40,60,80,100,120,140,160),main='히스토그램')

> rug(담즙의과포화비율)

> stem(담즙의과포화비율) #나무 줄기 그림

The decimal point is 1 digit(s) to the right of the |

2 | 5

4 | 072256788

6 | 556679334566789

8 | 000244666778889018

10 | 667001268

12 | 033787

14 | 26

> stem(담즙의과포화비율,scale=2) #줄기의 마디 2배로 늘리기

The decimal point is 1 digit(s) to the right of the |

3 | 5

4 | 07

5 | 2256788

6 | 556679

7 | 334566789

8 | 000244666778889

9 | 018

10 | 667

11 | 001268

12 | 03378

13 | 7

14 | 26

> library(vioplot)

> vioplot(담즙의과포화비율,names='담즙의과포화비율',col='yellow')

> n=length(담즙의과포화비율)

> plot(sort(담즙의과포화비율),(1:n)/n,type='s',ylim=c(0,1),main='Ogive of bile supersaturation',ylab='ECDF',xlab='담즙과포화비율') #ogive 오자이브 그래프

> rug(담즙의과포화비율)

출처: 보건정보데이터 분석(이태림, 이재원, 김주한, 장대흥 공저)

'KNOU > 2 보건 정보 데이터 분석' 카테고리의 다른 글

제4장 범주형 자료의 분석 - 4.3 로짓분석 (0)	2017.01.04
제4장 범주형 자료의 분석 - 4.2.4 대응자료 및 사례 - 대조군 검정 (0)	2016.12.29
제4장 범주형 자료의 분석 - 4.2.2 독립성 검정 (카이제곱 검정) (0)	2016.11.08
제4장 범주형 자료의 분석 - 4.2 범주형 자료의 검정(카이제곱 검정) (0)	2016.11.04
제4장 범주형 자료의 분석 - 4.1 범주형 자료와 분할표 (0)	2016.11.04

Posted by 마르띤

,

제8장 연관성분석

KNOU/2 데이터마이닝 2016. 12. 26. 14:06

1. 연관분석의 의의

연관분석(association analysis)은 연관 규칙을 통해 하나의 거래나 사건에 포함되어 있는 둘 이상의 품목 간 상호 연관성을 발견해 내는 것. 예를 들어 와인을 구매한 고객 중 15%가 치즈를 구매한다, 기저귀를 산 남성 고객은 맥주를 산다 등. 이러한 분석을 통하여 효율적인 매장 진열, 패키지 품목의 개발, 교차판매전략 등에 응용할 수 있다.

2. 연관 분석의 측도

(1) 지지율: A -> B의 지지율은 전체 거래 중 A와 B가 동시에 포함된 거래의 수. 전체 거래 수 5건 중 빵과 버터를 동시에 구매한 거래수가 3건이라면 지지율은 3/5 = 60%. 지지율이 1에 가까울수록 연관도가 높다고 할 수 있다.

지지율의 단점은 표본수가 적은 경우 통계적 유의성을 증명하기 어렵고, A->B 연관규칙과 B->A의 연관규칙 차이를 알 수 없으므로 다른 평가지표가 필요하다.

(2) 신뢰도: 지지율의 경우 기준이 되는 사건(구매)이 전체 집합인 데 비하여 신뢰도는 기준이 되는 사건을 특정 품목을 구매한 것에 한정한다. A->B 신뢰도는 A를 구매하였을 때, 품목 B를 추가로 구매할 확률. 예를 들어 장바구니 수가 150개였고, 빵과 우유가 함께 들어가 있는 장바구니가 30개였다면, 빵-> 우유의 지지율은 30/150 = 20%. 그런데 빵이 들어가 있는 장바구니만 추렸더니 100개였고, 이중 우유가 들어있는 것이 30개였다면, 빵-> 우유의 신뢰도는 30/100 = 30%. 신뢰도가 1에 가까울수록 연관도가 높다고 할 수 있다.

하지만 빵->우유의 신뢰도가 30%였으나, 빵이 없는 장바구니 50개 중 우유가 있는 경우가 15개라면 이 경우에도 신뢰도가 30%라, 또 다른 지표가 필요하다.

(3)향상도: 향상도는 A->B의 신뢰도를 독립 가정하에서의 신뢰도로 나눈 것을 의미한다. 우유를 살 확률 p(A)는 0.6, 콜라를 살 확률 p(B)는 0.5, 동시 구매할 지지율이 0.4, 신뢰도가 0.67일 경우, 향샹도는 신뢰도 / p(B) 즉 0.67/0.5 = 1.34가 된다. 이는 우유 구매 데이터를 안 상태에서 콜라 마케팅을 진행 할 경우, 우유 구매 데이터를 모를 때 보다 그 효과가 34% 증가한다고 알 수 있다.

l A->B 지지율 : 품목 A를 구매한 후 B를 구매한 거래 수 / 전체 거래수

l A->B 신뢰도 : 품목 A를 구매한 후 B를 구매한 거래 수 / A를 포함한 전체 거래 수

l A->B 향상도 : A->B 신뢰도 / B를 포함한 거래의 비중

3. 연관 분석의 장단점

l 장점: 쉽게 이해할 수 있고, 적용하기도 편하다. 목적변수 없는 직관적인 분석(자율학습 unsupervised learning에 해당)이라는 점도 쉽게 활용할 수 있는 장점. 데이터마이닝에 앞서 탐색 도구로도 유용

l 단점: 품목수의 증가에 따라 계산량이 증가. 유사한 품목은 한 범주로 일반화해야 함. 연속형 변수를 사용해서는 규칙을 찾기 힘들고, 거래가 드문 품목에 대한 정보를 찾기 어렵다.

4. R 실습

1) 데이터 입력

> setwd('c:/Rwork')

> read.table('trains3.txt')

V1

1 癤퓅ilk,bread,butter

2 milk,butter,coke

3 bread,butter,coke

4 milk,coke,ramen

5 bread,butter,ramen

> tr1=read.transactions('trains3.txt',format='basket',sep=',')

> tr1

transactions in sparse format with

6 transactions (rows) and

6 items (columns)

> as(tr1,'data.frame')

items

1 {bread,butter,癤퓅ilk}

2 {butter,coke,milk}

3 {bread,butter,coke}

4 {coke,milk,ramen}

5 {bread,butter,ramen}

6 {}

as 함수를 통해 tr1이라는 이름으로 저장되어 있는 거래 object를 출력해 볼 수 있다.

2) apriori 함수로 연관규칙을 산출.

> rules1=apriori(tr1,parameter=list(supp=0.4,conf=0.4))

Apriori

Parameter specification:

confidence minval smax arem aval originalSupport maxtime support minlen maxlen target ext

0.4 0.1 1 none FALSE TRUE 5 0.4 1 10 rules FALSE

Algorithmic control:

filter tree heap memopt load sort verbose

0.1 TRUE TRUE FALSE TRUE 2 TRUE

Absolute minimum support count: 2

set item appearances ...[0 item(s)] done [0.00s].

set transactions ...[6 item(s), 6 transaction(s)] done [0.00s].

sorting and recoding items ... [3 item(s)] done [0.00s].

creating transaction tree ... done [0.00s].

checking subsets of size 1 2 done [0.00s].

writing ... [5 rule(s)] done [0.00s].

creating S4 object ... done [0.00s].

> inspect(rules1)

lhs rhs support confidence lift

[1] {} => {bread} 0.5000000 0.5000000 1.0

[2] {} => {coke} 0.5000000 0.5000000 1.0

[3] {} => {butter} 0.6666667 0.6666667 1.0

[4] {bread} => {butter} 0.5000000 1.0000000 1.5

[5] {butter} => {bread} 0.5000000 0.7500000 1.5

함수 설명: rules1=apriori(tr1,parameter=list(supp=0.4,conf=0.4)) default값으로는 지지율support는 0.1, 신뢰도confidence는 0.8. bread

결과 해석: bread -> butter의 지지율은 0.5, 신뢰도는 1.0(100%), 향상율은 1.5(150%). 이는 빵을 구매한 고객은 반드시 버터를 사고 이 사실을 알고 버터에 대한 마케팅을 할 경우 모를 때 보다 그 효과가 25% 증가한다는 것을 의미한다.

3) 연관규칙의 시각화

> library(arulesViz)

> plot(rules1)

가로축은 지지율, 세로축은 신뢰도, 우측 막대기는 향상도

> plot(rules1,'grouped')

원의 크기는 지지율, 색상의 진하기는 향상도. LHS는 조건으로 이름 앞의 숫자는 각 규칙의 수를 의미

> plot(rules1,'graph')

원의 크기는 각 규칙의 지지율을, 색상의 진하기는 향상도를 의미. Butter와 bread가 상대적으로 중심에 위치하고 있고, 콜라의 경우 떨어져 있다.

> plot(rules1,'paracoord')

품목 간 연간관계를 병렬적으로 확인할 수 있다. 화살표의 두께는 지지율을, 화살표의 색상의 진하기는 향상도를 나타낸다.

출처: 데이터마이닝(장영재, 김현중, 조형준 공저)

소감: 데이터 양이 많을 경우 입력이 조금 힘들 것 같지만, A를 사는 고객은 반드시 B를 사니 A를 사지 않은 고객 대상 타켓 마케팅이 가능해보여 적용할 수 있는 분야가 많아 보인다. 그래프도 조금 더 공부는 해야겠지만, 시각화된 내용이 매우 직관적이라 활용도가 많아 보인다.

'KNOU > 2 데이터마이닝' 카테고리의 다른 글

제6장 모형비교평가 - 연속형 목표변수 (0)	2017.02.24
제5장 신경망모형 - 분류 (5)	2016.12.20
제5장 신경망모형 - 회귀 (0)	2016.11.14
제4장 앙상블 모형 - 분류앙상블모형 - 랜덤 포레스트 (0)	2016.11.09
제4장 앙상블 모형 - 분류앙상블모형 - 부스팅 (0)	2016.11.07

Posted by 마르띤

,

제5장 신경망모형 - 분류

KNOU/2 데이터마이닝 2016. 12. 20. 10:12

독일 신용평가 데이터를 활용한 신경망 모형. 목표변수 y는 good / bad의 범주형 데이터로 모든 변수를 수치화 한 후 신경망 모형을

1) 데이터 입력

> set.seed(1000)

> library(neuralnet)

> library(dummy)

> setwd('c:/Rwork')

> german = read.table('germandata.txt',header = T)

추가 공부: dummy화란?

2) 데이터 및 타입 변경

> dvar=c(4,9,10,15,17) #명목변수 지정 purpose(a43,a40..), personal, debtors, housing, job

> german2 = dummy(x=german[,dvar]) #명목변수를 더미변수화

> head(german2,1)

purpose_A40 purpose_A41 purpose_A410 purpose_A42 purpose_A43 purpose_A44 purpose_A45

1 0 0 0 0 1 0 0

purpose_A46 purpose_A48 purpose_A49 personal_A91 personal_A92 personal_A93 personal_A94

1 0 0 0 0 0 1 0

debtors_A101 debtors_A102 debtors_A103 housing_A151 housing_A152 housing_A153 job_A171

1 1 0 0 0 1 0 0

job_A172 job_A173 job_A174

1 0 1 0

> german2 = german2[,-c(10,14,17,20,24)] #더미변수생성

> head(german,1)

check duration history purpose credit savings employment installment personal debtors

1 A11 6 A34 A43 1169 A65 A75 4 A93 A101

Residence property age others housing numcredits job residpeople telephone foreign y

1 4 A121 67 A143 A152 2 A173 1 A192 A201 good

> german2 = cbind(german[,-dvar],german2) #변수 결함

> str(german2)
'data.frame':   1000 obs. of 40 variables:
$ check       : Factor w/ 4 levels "A11","A12","A13",..: 1 2 4 1 1 4 4 2 4 2 ...
$ duration    : int 6 48 12 42 24 36 24 36 12 30 ...
$ history     : Factor w/ 5 levels "A30","A31","A32",..: 5 3 5 3 4 3 3 3 3 5 ...
$ credit      : int 1169 5951 2096 7882 4870 9055 2835 6948 3059 5234 ...
$ savings     : Factor w/ 5 levels "A61","A62","A63",..: 5 1 1 1 1 5 3 1 4 1 ...
$ employment : Factor w/ 5 levels "A71","A72","A73",..: 5 3 4 4 3 3 5 3 4 1 ...
$ installment : int 4 2 2 2 3 2 3 2 2 4 ...
$ residence   : int 4 2 3 4 4 4 4 2 4 2 ...
$ property    : Factor w/ 4 levels "A121","A122",..: 1 1 1 2 4 4 2 3 1 3 ...
$ age         : int 67 22 49 45 53 35 53 35 61 28 ...
$ others      : Factor w/ 3 levels "A141","A142",..: 3 3 3 3 3 3 3 3 3 3 ...
$ numcredits : int 2 1 1 1 2 1 1 1 1 2 ...
$ residpeople : int 1 1 2 2 2 2 1 1 1 1 ...
$ telephone   : Factor w/ 2 levels "A191","A192": 2 1 1 1 1 2 1 2 1 1 ...
$ foreign     : Factor w/ 2 levels "A201","A202": 1 1 1 1 1 1 1 1 1 1 ...
$ y           : Factor w/ 2 levels "bad","good": 2 1 2 2 1 2 2 2 2 1 ...
$ purpose_A40 : Factor w/ 2 levels "0","1": 1 1 1 1 2 1 1 1 1 2 ...
$ purpose_A41 : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 2 1 1 ...
$ purpose_A410: Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
$ purpose_A42 : Factor w/ 2 levels "0","1": 1 1 1 2 1 1 2 1 1 1 ...
$ purpose_A43 : Factor w/ 2 levels "0","1": 2 2 1 1 1 1 1 1 2 1 ...
$ purpose_A44 : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
$ purpose_A45 : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
$ purpose_A46 : Factor w/ 2 levels "0","1": 1 1 2 1 1 2 1 1 1 1 ...
$ purpose_A48 : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
$ purpose_A49 : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
$ personal_A91: Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 2 1 ...
$ personal_A92: Factor w/ 2 levels "0","1": 1 2 1 1 1 1 1 1 1 1 ...
$ personal_A93: Factor w/ 2 levels "0","1": 2 1 2 2 2 2 2 2 1 1 ...
$ personal_A94: Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 2 ...
$ debtors_A101: Factor w/ 2 levels "0","1": 2 2 2 1 2 2 2 2 2 2 ...
$ debtors_A102: Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
$ debtors_A103: Factor w/ 2 levels "0","1": 1 1 1 2 1 1 1 1 1 1 ...
$ housing_A151: Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 2 1 1 ...
$ housing_A152: Factor w/ 2 levels "0","1": 2 2 2 1 1 1 2 1 2 2 ...
$ housing_A153: Factor w/ 2 levels "0","1": 1 1 1 2 2 2 1 1 1 1 ...
$ job_A171    : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
$ job_A172    : Factor w/ 2 levels "0","1": 1 1 2 1 1 2 1 1 2 1 ...
$ job_A173    : Factor w/ 2 levels "0","1": 2 2 1 2 2 1 2 1 1 1 ...
$ job_A174    : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 2 1 2 ...

> nrow(german2);ncol(german2)
[1] 1000
[1] 40

> for(i in 1:ncol(german2)) if(!is.numeric(german2[,i])) german2[,i] = as.numeric(german2[,i])#여타 순서가 있는 범주형 변수의 수치형 변수화

> german2$y = ifelse(german$y == 'good',1,0) #목표변수 변환

> head(german$y)

[1] good bad good good bad good

Levels: bad good

> head(german2$y)

[1] 1 0 1 1 0 1

## 중요 : 신경망에서는 범주형 데이터를 수치화하여 적용한다. ##

3) 75% 랜덤 추출

> i = sample(1:nrow(german2),round(0.75*nrow(german2))) #75%랜덤 추출

> length(i)
[1] 750

4) 변수의 표준화 과정

> max2 = apply(german2, 2, max)

> min2 = apply(german2, 2, min)

> gdat = scale(german2, center = min2, scale = max2 - min2) # 변수조정(0,1,dummy variable은 변화 없음)

> gdat = as.data.frame(gdat) #데이터 프레임 형태로 변경, 데이터 준비 끝!

> str(gdat)

'data.frame': 1000 obs. of 35 variables:

$ check : num 0 0.333 1 0 0 ...

$ duration : num 0.0294 0.6471 0.1176 0.5588 0.2941 ...

$ history : num 1 0.5 1 0.5 0.75 0.5 0.5 0.5 0.5 1 ...

$ credit : num 0.0506 0.3137 0.1016 0.4199 0.2542 ...

$ savings : num 1 0 0 0 0 1 0.5 0 0.75 0 ...

$ employment : num 1 0.5 0.75 0.75 0.5 0.5 1 0.5 0.75 0 ...

$ installment : num 1 0.333 0.333 0.333 0.667 ...

$ residence : num 1 0.333 0.667 1 1 ...

$ property : num 0 0 0 0.333 1 ...

$ age : num 0.8571 0.0536 0.5357 0.4643 0.6071 ...

$ others : num 1 1 1 1 1 1 1 1 1 1 ...

$ numcredits : num 0.333 0 0 0 0.333 ...

$ residpeople : num 0 0 1 1 1 1 0 0 0 0 ...

$ telephone : num 1 0 0 0 0 1 0 1 0 0 ...

$ foreign : num 0 0 0 0 0 0 0 0 0 0 ...

$ y : num 1 0 1 1 0 1 1 1 1 0 ...

$ purpose_A40 : num 0 0 0 0 1 0 0 0 0 1 ...

$ purpose_A41 : num 0 0 0 0 0 0 0 1 0 0 ...

$ purpose_A410: num 0 0 0 0 0 0 0 0 0 0 ...

$ purpose_A42 : num 0 0 0 1 0 0 1 0 0 0 ...

$ purpose_A43 : num 1 1 0 0 0 0 0 0 1 0 ...

$ purpose_A44 : num 0 0 0 0 0 0 0 0 0 0 ...

$ purpose_A45 : num 0 0 0 0 0 0 0 0 0 0 ...

$ purpose_A46 : num 0 0 1 0 0 1 0 0 0 0 ...

$ purpose_A48 : num 0 0 0 0 0 0 0 0 0 0 ...

$ personal_A91: num 0 0 0 0 0 0 0 0 1 0 ...

$ personal_A92: num 0 1 0 0 0 0 0 0 0 0 ...

$ personal_A93: num 1 0 1 1 1 1 1 1 0 0 ...

$ debtors_A101: num 1 1 1 0 1 1 1 1 1 1 ...

$ debtors_A102: num 0 0 0 0 0 0 0 0 0 0 ...

$ housing_A151: num 0 0 0 0 0 0 0 1 0 0 ...

$ housing_A152: num 1 1 1 0 0 0 1 0 1 1 ...

$ job_A171 : num 0 0 0 0 0 0 0 0 0 0 ...

$ job_A172 : num 0 0 1 0 0 1 0 0 1 0 ...

$ job_A173 : num 1 1 0 1 1 0 1 0 0 0 ...

5) 신경망 모델 구축 및 신경망 그래프 그리기

> train = gdat[i,] #학습샘플과 테스트 샘플 추출

> test = gdat[-i,]

> gn = names(german2)

> gn

[1] "check" "duration" "history" "credit" "savings" "employment"

[7] "installment" "residence" "property" "age" "others" "numcredits"

[13] "residpeople" "telephone" "foreign" "y" "purpose_A40" "purpose_A41"

[19] "purpose_A410" "purpose_A42" "purpose_A43" "purpose_A44" "purpose_A45" "purpose_A46"

[25] "purpose_A48" "personal_A91" "personal_A92" "personal_A93" "debtors_A101" "debtors_A102"

[31] "housing_A151" "housing_A152" "job_A171" "job_A172" "job_A173"

> f = as.formula(paste('y~',paste(gn[!gn %in% 'y'],collapse = '+')))

> f

y ~ check + duration + history + credit + savings + employment +

installment + residence + property + age + others + numcredits +

residpeople + telephone + foreign + purpose_A40 + purpose_A41 +

purpose_A410 + purpose_A42 + purpose_A43 + purpose_A44 +

purpose_A45 + purpose_A46 + purpose_A48 + personal_A91 +

personal_A92 + personal_A93 + debtors_A101 + debtors_A102 +

housing_A151 + housing_A152 + job_A171 + job_A172 + job_A173

> nn1 = neuralnet(f,data=train,hidden=c(3,2),linear.output=F) #은닉층은 2개, 첫번째 노드는 3개, 두번째 노드는 2개. 분류의 경우 linear.output = F

> summary(nn1)

Length Class Mode

call 5 -none- call

response 750 -none- numeric

covariate 25500 -none- numeric

model.list 2 -none- list

err.fct 1 -none- function

act.fct 1 -none- function

linear.output 1 -none- logical

data 35 data.frame list

net.result 1 -none- list

weights 1 -none- list

startweights 1 -none- list

generalized.weights 1 -none- list

result.matrix 119 -none- numeric

> plot(nn1)

은닉층이 2개인 신경망 모형의 그래프가 완성된다. 이 그래프에 대한 해석을 좀 더 공부하고 싶은데 아직은 잘 모르겠음.

6) 모형 추정:

> dim(german2)[1]
[1] 1000

> dim(german2)[2]
[1] 35

> colnames(test)[16]
[1] "y"

> pred.nn0 = compute(nn1,train[,c(1:15,17:dim(german2)[2])]) #학습데이터의 실제값과 예측값 비교, 16번째 열의 값은 y

함수 설명:

> pred.nn0 = compute(nn1,train[,c(1:15,17:dim(german2)[2])]) 16번째 변수가 y 목표변수. compute는 작성된 신경망모형을 이용하여 새로운 예에 적용하여 결과를 도출. nn1는 새롭게 예측에 적용할 자료, train[,c(1:15,17:dim(german2)[2])]는 신경망모형으로적합한 결과 오브젝트

7) 학습샘플의 실제값과 예측값을 비교해보자.

> head(cbind(german2[1,16],round(pred.nn0$net.result,10)))

[,1] [,2]

931 1 0.7786470581

546 1 0.0000005387

56 1 0.8208161458

883 1 0.9999999722

11 1 0.0000004232

354 1 0.0000046419

#왼쪽이 실제값, 오른쪽이 예측값. 분류의 문제이므로 값은 0과 1사이. 0.5를 cut off값을 둘 수 있다. 또는 0.3미만 폐기, 0.7이하 보류, 0.7 초과만 사용한느 cutoff도 가능. 왼쪽이 실제 값, 오른쪽이 학습된 데이터. 4번째 행은 실제 1의 값을 1에 가깝게 예측하였고, 5번째 행은 실제 1이지만 0에 가깝게 예측한 사례. german2[,16]은 16번째 열 즉 y값임을 알겠는데 german2[1,16]은 뭘까… 궁금

8) 예측 정확도 평가

> pred.nn1 = compute(nn1,test[,c(1:15,17:dim(german2)[2])]) #test data를 바탕으로 판단해보자

> pred.nn2 = ifelse(pred.nn1$net.result>0.5,1,0) #0.5를 경계로 1과 0 분류

> head(cbind(german2[-i,16],pred.nn2)) #테스트 샘플의 실제값과 예측값

[,1] [,2]

3 1 1

12 1 0

13 1 1

14 1 0

15 0 1

26 1 1

> sum(german2[-i,16]!=pred.nn2) / length(german2[-i,16])

[1] 0.404

> sum(german2[-i,16]!=pred.nn2)

[1] 101

> length(german2[-i,16])

[1] 250

#테스트 샘플의 오분류율, 16번째 값은 목표변수, sum(german2[-i,16]!=pred.nn2)는 pred.nn2와 같지 않은 값을 전체길이 length(german2[-i,16])로 나눔. i를 빼고 16번째 컬럼을 사용

> library(nnet)

> nnet1 =nnet(f,data=train,size = 3, linout = F)

# weights: 109

initial value 181.570914

iter 10 value 113.679457

iter 20 value 96.943318

iter 30 value 82.803121

iter 40 value 73.239858

iter 50 value 70.807278

iter 60 value 69.865795

iter 70 value 69.476434

iter 80 value 69.158350

iter 90 value 69.039026

iter 100 value 68.929899

final value 68.929899

stopped after 100 iterations

> pred.nnet1 = predict(nnet1,test[,c(1:15,17:dim(german2)[2])])

> pred.nnet2 = ifelse(pred.nnet1>0.5,1,0)

> head(cbind(german2[-i,16],pred.nnet2)) #테스트 샘플의 실제값과 예측값

[,1] [,2]

3 1 1

12 1 0

13 1 1

14 1 1

15 0 0

26 1 1

> sum(german2[-i,16]!=pred.nnet2) / length(german2[-i,16]) #테스트 샘플 예측의 오분류율

[1] 0.408

이 부분은 교재에 별도 설명이 없어서 추가 공부가 필요함.

소감: 알파고 딥마이닝으로 인해 관심을 가지게 된 신경망 모형. 이론 공부도 해보고 R도 따라해보니 약 50%정도 이해된 상태. 궁금한 점은 아래와 같음.

1. 명복 변수 중 더미화 하지 않은 것들도 있음.

-> 교수님 답변: 해당 신경망 모형에서는 숫자로 입력되어 있는 범주형 변수들은 수치변수로 그대로 사용하고 표준화만 하였습니다. 순서가 있는 범주형 변수라고 판단한 변수였습니다. 해당 변수들을 제외하고 나머지 변수들은 변환이 필요하여 dummy 함수를 사용하였습니다.

2. 위 신경망 plot 그래프가 무엇을 의미하는지 더 자세히 해석할 능력이 필요함.

-> 교수님 답변: 신경망 모형의 해석은 그림을 보고 해석하기가 상당히 힘듭니다. 워낙 복잡한 함수의 결합이기 때문입니다. 다만, 화살표 위의 수치(절대값)를 보고 연결강도가 강한지 여부를 판단할 수 있습니다. 신경망 모형의 태생적인 한계점인 것 같습니다.

참고로 교재에서 사용하였던 neuralnet 패키지는 plot을 제공합니다. 과거 강의에서 사용했던 다른 패키지는 직접적으로 plot을 산출할 수는 없었습니다.

3. 위 신경망 모형을 통해서 오분류율은 40.4%로 나왔는데 너무 높은건 아닌지 생각됨.

-> 교수님 답변: 오분율은 상대적인 것이긴 하지만 높은 수준으로 보입니다. 교재에서는 예측 판별 기준으로 0.5라는 값을 사용했는데 실제 실무에서는 이 값을 적절하게 변경시키면서 오분류율을 낮추어주는 것이 좋을 것입니다. 예를 들어 원 데이터가 1이 0에 비해 많이 분포되어 있다면, 0.5보다 큰 값을 기준치로 삼는 편이 좋습니다.

4. 은닉층의 개수는 어떻게 설정하면 좋을지?

5. 가장 궁금한 것은 이 모델을 어떻게 실무에 적용할지 잘 모르겠음. 위 독일 신용평가 데이터로 신경망 모형을 만들고 오분류율도 체크하고, 실제값과 예측값도 비교하였는데, 이것이 의미하는 것들. 가령, 변수의 중요성, 그래서 신용도가 좋은 경우는 어느 경우이고, 또 다른 고객 데이터 셋이 있는데 이 새로운 셋에서 어떤 고객, 변수, 값들을 추출해야 우리가 원하는 우수한 고객을 알아낼 수 있는지, 실제 비즈니스 적용 포인트가 가장 궁금함.

-> 교수님 답변: 실무에서 신경망 모형을 잘 해석하여 활용하기는 어려운 것이 사실입니다. 다만, 상대적으로 예측력이 높다는 점을 살려 새로운 데이터의 모든 변수들을 활용하여 1이나 0값을 예측하거나 목표변수의 값 자체를 산출하는 목적으로는 유용성이 높다고 봅니다. 변수의 해석이나 변수의 선택보다는 예측이나 분류 자체의 목적으로 사용하기에 적합하다고 보시면 좋습니다.

6. 전반적으로 의사결정나무, 앙상블 모형에 이어 데이터 마이닝 분야에 깊은 관심을 가지게 되는 좋은 계기.

출처: 데이터마이닝(장영재, 김현중, 조형준 공저)

'KNOU > 2 데이터마이닝' 카테고리의 다른 글

제6장 모형비교평가 - 연속형 목표변수 (0)	2017.02.24
제8장 연관성분석 (0)	2016.12.26
제5장 신경망모형 - 회귀 (0)	2016.11.14
제4장 앙상블 모형 - 분류앙상블모형 - 랜덤 포레스트 (0)	2016.11.09
제4장 앙상블 모형 - 분류앙상블모형 - 부스팅 (0)	2016.11.07

Posted by 마르띤

,

제5장 신경망모형 - 회귀

KNOU/2 데이터마이닝 2016. 11. 14. 14:11

보스턴 하우징을 이용한 신경망 모형. 목표변수를 주택가격의 중간값인 medv(연속형 변수)로 하고, 나머지 변수를 입력변수로 하는 신경망 모형

1) 데이터 입력

> set.seed(100)

> library(MASS)

> library(neuralnet)

> bdat = Boston

> str(bdat)

'data.frame': 506 obs. of 15 variables:

$ crim : num 0.00632 0.02731 0.02729 0.03237 0.06905 ...

$ zn : num 18 0 0 0 0 0 12.5 12.5 12.5 12.5 ...

$ indus : num 2.31 7.07 7.07 2.18 2.18 2.18 7.87 7.87 7.87 7.87 ...

$ chas : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...

$ nox : num 0.538 0.469 0.469 0.458 0.458 0.458 0.524 0.524 0.524 0.524 ...

$ rm : num 6.58 6.42 7.18 7 7.15 ...

$ age : num 65.2 78.9 61.1 45.8 54.2 58.7 66.6 96.1 100 85.9 ...

$ dis : num 4.09 4.97 4.97 6.06 6.06 ...

$ rad : Factor w/ 9 levels "1","2","3","4",..: 1 2 2 3 3 3 5 5 5 5 ...

$ tax : num 296 242 242 222 222 222 311 311 311 311 ...

$ ptratio : num 15.3 17.8 17.8 18.7 18.7 18.7 15.2 15.2 15.2 15.2 ...

$ black : num 397 397 393 395 397 ...

$ lstat : num 4.98 9.14 4.03 2.94 5.33 ...

$ medv : num 24 21.6 34.7 33.4 36.2 28.7 22.9 27.1 16.5 18.9 ...

$ medv.hat: num 23.7 23.7 32.4 32.4 32.4 ...

주의: chas, rad를 수치형 변수로 바꾸고, 앙상블 모형에서 사용하였던 medv.hat 변수를 삭제해야 한다.

신경망에서는 범주형 데이터를 수치화하여 적용한다. 대체로 순서가 의미 있는 범부형 데이터는 수치 전환을 한 뒤 표준화하여 사용하게 된다.

2) 데이터 및 타입 변경

> bdat<-bdat[,-15]

> bdat$chas = as.numeric(bdat$chas)

> bdat$rad = as.numeric(bdat$rad)

> class(bdat$chas);class(bdat$rad)

[1] "numeric"

3) 50% 랜덤 추출

> i = sample(1:nrow(bdat), round(0.5*nrow(bdat)))

> max1 = apply(bdat, 2, max)

> min1 = apply(bdat, 2, min)

결과 해석 : 50% 랜덤 추출하여 훈련데이터(train)로 저장하고 나머지는 검증 데이터(test)로 저장

4) 변수의 표준화 과정

> sdat = scale(bdat,center=min1,scale = max1 - min1)

> sdat = as.data.frame(sdat)

> head(sdat,3)

crim zn indus chas nox rm age dis rad

1 0.0000000000000 0.18 0.06781524927 0 0.3148148148 0.5775052692 0.6416065911 0.2692031391 0.000

2 0.0002359225392 0.00 0.24230205279 0 0.1728395062 0.5479977007 0.7826982492 0.3489619802 0.125

3 0.0002356977440 0.00 0.24230205279 0 0.1728395062 0.6943858977 0.5993820803 0.3489619802 0.125

tax ptratio black lstat medv

1 0.2080152672 0.2872340426 1.0000000000 0.08967991170 0.4222222222

2 0.1049618321 0.5531914894 1.0000000000 0.20447019868 0.3688888889

3 0.1049618321 0.5531914894 0.9897372535 0.06346578366 0.6600000000

함수 설명:

> sdat = scale(bdat,center=min1,scale = max1 - min1) 신경망에서 사용하는 수치형 변수를 0과 1사이의 값으로 바꾸어 주기 위한 단계. 위의 [0-1] 변환은 (x-min(x) /(max(x)-min(x)) 의 수식으로도 계산 가능.

5) 신경망 모델 구축 및 신경망 그래프 그리기

> train = sdat[i,] #학습, 훈련샘플 training

> test = sdat[-i,] #테스트샘플 test

> n = names(train) ;n

[1] "crim" "zn" "indus" "chas" "nox" "rm" "age" "dis" "rad" "tax"

[11] "ptratio" "black" "lstat" "medv"

> form = as.formula(paste('medv~',paste(n[!n %in% 'medv'],collapse = '+')))

> form

medv ~ crim + zn + indus + chas + nox + rm + age + dis + rad + tax + ptratio + black + lstat

> nn1 = neuralnet(form,data=train,hidden = c(5,3),linear.output = T)

> summary(nn1) #작성된 신경망모형의 오브젝트의 구조를 정리

Length Class Mode

call 5 -none- call

response 253 -none- numeric

covariate 3289 -none- numeric

model.list 2 -none- list

err.fct 1 -none- function

act.fct 1 -none- function

linear.output 1 -none- logical

data 14 data.frame list

net.result 1 -none- list

weights 1 -none- list

startweights 1 -none- list

generalized.weights 1 -none- list

result.matrix 95 -none- numeric

함수 설명

> form = as.formula(paste('medv~',paste(n[!n %in% 'medv'],collapse = '+'))) #종속변수는 mdev, 나머지는 독립변수이다 틸다 ~. 와같은 개념. 종속변수를 제외한 n 안에 있는 모든 변수명을 다 집어넣고 더하여라. 변수명을 직접 입력해도 된다.

실제 form라고 입력하면 medv ~ crim + zn + indus + chas + nox + rm + age + dis + rad + tax + ptratio + black + lstat

라는 공식을 볼 수 있다.

> nn1 = neuralnet(form,data=train,hidden = c(5,3),linear.output = T) #은닉층은 2개이고 처음은 5개 나머지는 3개, 회귀의 문제이므로 linear.output = T

질문: 은닉층의 수는 어떻게 결정할까? 통상 은닉층의 마디수는 입력층의 마디수 두배를 넘지 않도록 해야한다. 위에서는 c(5,3)으로 결정하였는데 이 수치가 최적의 수치인지는 어떻게 알까? 더 공부가 필요하다.

<참고> 딥러닝(deep learning)

기계학습 기법 중 하나로 신경망모형으로부터 비롯된 딥러닝(deep learning)은 기본적으로 은닉층이 많이 쌓여 가면서 복잡하고 깊은 구조로 발전하면서 deep 이라는 이름이 붙여졌다. 입력변수와 출력변수 간 복잡한 관계를 가중치를 통해 조정할 수 있는 구조와 매커니즘을 갖고 있다.

> plot(nn1)

①로 표시된 것이 상수항에 해당하고 가중치는 각각의 화살 표 위에 출력된다. 해석 공부도 더 필요하다.

처음 은닉층은 5개, 두번째 은닉층은 3개를 가지고 있는 신경망 모형이 완성된다.

6) 모형 추정: 위의 그래프는 50%의 training data만을 사용하였으니 나머지 50% test data를 쓰자.

> pred.nn0 = compute(nn1,test[,1:13])

> summary(pred.nn0)
               Length Class Mode
    neurons      3    -none- list
    net.result 253    -none- numeric

> names(pred.nn0)

[1] "neurons" "net.result"

> pred0 = pred.nn0$net.result*(max(bdat$medv)-min(bdat$medv))+min(bdat$medv)

> head(pred0)
                    [,1]
    1 26.55282080
    2 25.15912380
    3 34.48852794
    4 31.56317828
    5 32.47594868
    6 25.32302689

함수 설명:

> pred.nn0 = compute(nn1,test[,1:13])

14번째 변수가 medv 목표변수.compute는 작성된 신경망모형을 이용하여 새로운 예에 적용하여 결과를 도출. nn1는 새롭게 예측에 적용할 자료, test[,1:13]는 신경망모형으로 적합한 결과 오브젝트

> pred0 = pred.nn0$net.result*(max(bdat$medv)-min(bdat$medv))+min(bdat$medv)

목표변수를 조정전 변수값으로 환원하는 과정. 원래 가지고 있는 값으로 환원. 역변환 과정

7) 학습샘플의 실제값과 예측값 (적합값)을 비교해보자.

아래 함수 head(cbind(bdat[i,14],pred0),50)에서 cbind(bdat[i,14],는 실제값(training data)과 예측값(neural network)으로 구성된 행렬을 구성하는 함수

> head(cbind(bdat[i,14],pred0),50)

[,1] [,2]

1 15.6 28.12482177

2 19.2 21.36405177

3 29.1 34.24035993

4 18.4 35.88436410

6 24.0 25.77273199

9 22.2 14.86779840

12 11.9 20.82595913

14 26.4 20.03825755

18 24.4 17.19077183

19 23.9 18.69934412

21 20.3 14.85391889

23 9.6 15.77921679

24 13.3 14.95074253

25 33.3 16.19467498

26 15.0 15.53094329

27 19.3 16.50978025

28 27.5 15.65043905

30 22.6 19.71434020

33 29.4 13.39051466

34 19.5 15.45310394

36 33.8 22.75202201

41 31.2 37.96410157

42 21.2 36.32787047

43 19.9 27.29495234

44 42.3 26.74797523

46 24.8 20.55343697

47 50.0 19.83522653

48 20.8 16.66204238

50 48.8 16.90901569

51 23.0 19.54979162

54 41.7 22.43724314

55 17.1 16.69040814

57 25.0 22.91956778

58 15.2 29.22726384

62 8.1 18.82030494

63 8.8 24.13500281

66 19.7 25.15181243

67 28.6 20.16650282

68 20.2 20.27218139

69 18.7 17.48036500

70 17.0 19.87776814

71 12.1 25.21432307

73 25.0 23.16314867

75 12.3 23.82067236

77 23.9 20.42068536

79 37.6 20.09872166

80 22.7 20.89416546

81 5.0 28.30078783

84 28.4 22.10685814

87 14.0 21.30514814

결과 해석: 26번째 행의 값처럼 유사한 결과값도 있지만 81번 결과 값처럼 서로 차이가 나는 경우도 있다.

8) 예측 정확도 평가: 모형 추정에 사용되지 않은 test data를 가지고 예측을 해보자.

> pred.nn1 = compute(nn1,test[,1:13]) #목표 변수 제외, 새로운 변수로 간주하고 nn1에 적용

> pred1 = pred.nn1$net.result*(max(bdat$medv)-min(bdat$medv))+min(bdat$medv) #목표변수를 조정전 변수값으로 환원. 역변환

> head(cbind(bdat[-i,14],pred1),50) #좌측은 test 변수의 실제값, 우측은 예측값(적합값), 14번째 값은 목표변수, pred1은 예측값(적합값)

[,1] [,2]

1 24.0 28.12482177

2 21.6 21.36405177

3 34.7 34.24035993

4 33.4 35.88436410

6 28.7 25.77273199

9 16.5 14.86779840

12 18.9 20.82595913

14 20.4 20.03825755

18 17.5 17.19077183

19 20.2 18.69934412

21 13.6 14.85391889

23 15.2 15.77921679

24 14.5 14.95074253

25 15.6 16.19467498

26 13.9 15.53094329

27 16.6 16.50978025

28 14.8 15.65043905

30 21.0 19.71434020

33 13.2 13.39051466

34 13.1 15.45310394

36 18.9 22.75202201

41 34.9 37.96410157

42 26.6 36.32787047

43 25.3 27.29495234

44 24.7 26.74797523

46 19.3 20.55343697

47 20.0 19.83522653

48 16.6 16.66204238

50 19.4 16.90901569

51 19.7 19.54979162

54 23.4 22.43724314

55 18.9 16.69040814

57 24.7 22.91956778

58 31.6 29.22726384

62 16.0 18.82030494

63 22.2 24.13500281

66 23.5 25.15181243

67 19.4 20.16650282

68 22.0 20.27218139

69 17.4 17.48036500

70 20.9 19.87776814

71 24.2 25.21432307

73 22.8 23.16314867

75 24.1 23.82067236

77 20.0 20.42068536

79 21.2 20.09872166

80 20.3 20.89416546

81 28.0 28.30078783

84 22.9 22.10685814

87 22.5 21.30514814

결과해석: 전반적으로 실제값과 예측값이 서로 유사함을 알 수 있다.

> head(cbind(bdat[-i,14],pred1)[,1],10)
   1    2    3    4    5    6    7    8    9   10
24.0 21.6 34.7 33.4 36.2 28.7 22.9 27.1 16.5 18.9
> head(cbind(bdat[-i,14],pred1),10)
   [,1]        [,2]
1 24.0 26.55282080
2 21.6 25.15912380
3 34.7 34.48852794
4 33.4 31.56317828
5 36.2 32.47594868
6 28.7 25.32302689
7 22.9 19.92609580
8 27.1 20.01203029
9 16.5 19.68915018
10 18.9 19.33552651
> obs <-cbind(bdat[-i,14],pred1)[,1]
> pdt <-cbind(bdat[-i,14],pred1)[,2]

> out = cbind(obs,pdt)
> head(out)
   obs         pdt
1 24.0 26.55282080
2 21.6 25.15912380
3 34.7 34.48852794
4 33.4 31.56317828
5 36.2 32.47594868
6 28.7 25.32302689
> which.min(abs(out$obs - out$pdt))
Error in out$obs : $ operator is invalid for atomic vectors
> out = as.data.frame(out)
> which.min(abs(out$obs - out$pdt))
[1] 43
> out[43,]
    obs         pdt
43 25.3 25.31048382
> which.max(abs(out$obs - out$pdt))
[1] 197
> out[197,]
    obs         pdt
372 50 18.13253008

9) PMSE: 예측된 값이 얼마나 잘 예측되었을까?

> PMSE = sum((bdat[-i,14] - pred1)^2) / nrow(test)

> PMSE

[1] 18.70948041

함수 설명:

> PMSE = sum((bdat[-i,14] - pred1)^2) / nrow(test) #예측 평균 제곱 오차,

결과 해석: 예측평균제곱오차는 18.7

출처: 데이터마이닝(장영재, 김현중, 조형준 공저)

'KNOU > 2 데이터마이닝' 카테고리의 다른 글

제8장 연관성분석 (0)	2016.12.26
제5장 신경망모형 - 분류 (5)	2016.12.20
제4장 앙상블 모형 - 분류앙상블모형 - 랜덤 포레스트 (0)	2016.11.09
제4장 앙상블 모형 - 분류앙상블모형 - 부스팅 (0)	2016.11.07
제4장 앙상블 모형 - 분류앙상블모형 - 배깅 (0)	2016.11.03

Posted by 마르띤

,

데이터마이너를 꿈꾸며

'KNOU'에 해당되는 글 52건

제6.2장 비모수적 방법 - 2. 누적한계추정법

'KNOU > 2 보건 정보 데이터 분석' 카테고리의 다른 글

제6.2장 비모수적 방법 - 1. 생명표 방법

'KNOU > 2 보건 정보 데이터 분석' 카테고리의 다른 글

제5.1장 공분산 분석 - 2 공변량이 둘 이상인 경우

'KNOU > 2 보건 정보 데이터 분석' 카테고리의 다른 글

제5.1장 공분산 분석 - 1. 공변량이 하나인 경우

'KNOU > 2 보건 정보 데이터 분석' 카테고리의 다른 글

제4장 범주형 자료의 분석 - 4.3 로짓분석

'KNOU > 2 보건 정보 데이터 분석' 카테고리의 다른 글

제4장 범주형 자료의 분석 - 4.2.4 대응자료 및 사례 - 대조군 검정

'KNOU > 2 보건 정보 데이터 분석' 카테고리의 다른 글

제2장 보건정보 데이터의 기초분석

'KNOU > 2 보건 정보 데이터 분석' 카테고리의 다른 글

제8장 연관성분석

'KNOU > 2 데이터마이닝' 카테고리의 다른 글

제5장 신경망모형 - 분류

'KNOU > 2 데이터마이닝' 카테고리의 다른 글

제5장 신경망모형 - 회귀

'KNOU > 2 데이터마이닝' 카테고리의 다른 글

링크

카테고리

최근에 올라온 글

최근에 받은 트랙백

글 보관함

티스토리툴바