데이터마이너를 꿈꾸며

제3장 연속형 자료의 분석

3.1 두 집단의 평균비교

3.1.1 독립표본의 평균비교 two sample test

예제) 흡연자 집단과 비흡연자 집단 간 폐 파괴지수를 측정하였다. 높은 수치는 폐의 손상이 크다는 것을 뜻한다. 흡연자와 비흡연자의 폐 파괴지수의 평균이 같다고 할 수 있는가? (각 그룹에서의 관측치들은 정규분포를 따르는 모집단으로부터 독립적으로 얻어진 것이며 두 그룹에서의 모분산은 같다고 가정하자. )

귀무가설 H0: 흡연자와 비흡연자의 폐 파괴지수 평균은 같다.

대립가설 H1: 흡연자와 비흡연자의 폐 파괴지수 평균은 다르다. (양측 검정)

#1. 자료 입력

> smoke=c(16.6,13.9,11.3,26.5,17.4,15.3,15.8,12.3,18.6,12,24.1,16.5,21.8,16.3,23.4,18.8)

> nonsmoke=c(18.1,6,10.8,11,7.7,17.9,8.5,13,18.9)

> sapply(list(smoke,nonsmoke),mean)

[1] 17.53750 12.43333

> sapply(list(smoke,nonsmoke),sd)

[1] 4.475247 4.849227

#2. 정규성 검정

> qqnorm(smoke,main='smoke')

> qqline(smoke,col='blue')

> shapiro.test(smoke)

Shapiro-Wilk normality test

data:  smoke

W = 0.94511, p-value = 0.4163

결과 해석: shapiro.test의 결과에 따라 p value = 0.4163 > 0.05 이므로 귀무가설 기각 못한다, 즉 정규분포를 따른다

> qqnorm(nonsmoke,main = 'nonsmoke')

> qqline(nonsmoke,col='red')

> shapiro.test(nonsmoke)

Shapiro-Wilk normality test

data:  nonsmoke

W = 0.90366, p-value = 0.274

#boxplot과 vioplot

> boxplot(smoke,nonsmoke,col='yellow',names=c('smoke','nonsmoke'))

> library(vioplot)

> vioplot(smoke,nonsmoke,col='yellow',names=c('smoke','nonsmoke'))

결과 해석:  두 집단에 차이가 있음을 알 수 있다.

#3. 두 모분산 비교 (양측검정)

#대립 가설의 형태: alternative = c('two.sided','less','greater')

> var.test(smoke,nonsmoke)

F test to compare two variances

data:  smoke and nonsmoke

F = 0.8517, num df = 15, denom df = 8, p-value = 0.7498

alternative hypothesis: true ratio of variances is not equal to 1

95 percent confidence interval:

  0.2076714 2.7243799

sample estimates:

  ratio of variances

0.8517046

결과 해석: p value 가 0.7498로 분산이 같다는 귀무가설을 기각할 수 없다. 즉, 등분산 가정

#4.두 모분산 비교 (양측검정) - 등분산 가정

> t.test(smoke,nonsmoke,var.equal = T)

Two Sample t-test

data:  smoke and nonsmoke

t = 2.658, df = 23, p-value = 0.01405

alternative hypothesis: true difference in means is not equal to 0

95 percent confidence interval:

  1.131680 9.076653

sample estimates:

  mean of x mean of y

17.53750  12.43333

결과해석:

귀무가설 H0: 흡연자와 비흡연자의 폐 파괴지수 평균은 같다.

대립가설 H1: 흡연자와 비흡연자의 폐 파괴지수 평균은 다르다. (양측 검정)

결정: p value는 p value는 0.01405로 두 모 평균이 같다는 귀무가설을 기각한다. 즉, 두 모평균이 서로 다르다.

#5.두 모분산 비교 (양측검정) - 이분산 가정

> t.test(smoke,nonsmoke)

Welch Two Sample t-test

data:  smoke and nonsmoke

t = 2.5964, df = 15.593, p-value = 0.01978

alternative hypothesis: true difference in means is not equal to 0

95 percent confidence interval:

  0.9279143 9.2804190

sample estimates:

  mean of x mean of y

17.53750  12.43333

결과해석: 등분산 가정과 큰 차이는 없다.

3.1.2 짝지은 표본의 평균비교 paired sample test

예) 환자 15명에게 혈압강하제를 12주 투입 후 혈압을 비교하였다. 새로운 약은 효과적인가?

귀무가설 h0: u1-u2 = 0

대립가설 h1: u1 > u2

1) 데이터 입력

> before=c(90,56,49,64,65,88,62,91,74,93,55,71,54,64,54)

> after=c(72,55,56,57,62,79,55,72,73,74,58,59,58,71,61)

> diff = before - after

2) 정규성 차이: shapiro – wilk test

> qqnorm(diff)

> qqline(diff,col='red')

> shapiro.test(diff)

Shapiro-Wilk normality test

data:  diff

W = 0.90982, p-value = 0.1345

결과 해석: shapiro test 결과 p value는 0.1345로써 정규분포를 이루고 있다고 할 수 있다.

3) Paired sample test

> mean(diff) ; sd(diff)

[1] 4.533333

[1] 9.425396

> t.test(before, after, paired = T, alternative = 'greater') #μ 복용 전 > μ 복용 후

Paired t-test

data:  before and after

t = 1.8628, df = 14, p-value = 0.0418

alternative hypothesis: true difference in means is greater than 0

95 percent confidence interval:

  0.2469617       Inf

sample estimates:

  mean of the differences

4.533333

귀무가설 H0: 약 복용전과 복용 후의 혈압 수치는 같다. μ 복용 전 = μ 복용 후

대립가설 H1: 약 복용전 대비 복용 후의 혈압 수치기 다 낮다. μ 복용 전 > μ 복용 후

결론: 단측 검정에 대한 p value가 0.0418로서 유의수준 5%에서 두 그룹의 혈압 차이가 없다는 귀무가설을 기각할 만한 충분한 증거가 있으므로 새로운 약이 혈압을 내린다고 볼 수 있다.

출처: 보건 정보 데이터 분석(이태림 저자)

목표변수가 집단을 의미하는 범주형 의사결정나무 -> 분류나무모형

목표변수가 연속형 변수인 의사결정나무 -> 회귀나무모형

예제) 목표변수가 범주형 good, bad인 독일 신용평가 데이터를 이용하여 cart 방법을 이용한 의사결정나무 구축

1) 데이터 불러오기

> setwd('c:/Rwork')

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

> str(german)

'data.frame':        1000 obs. of  21 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 ...

$ purpose    : Factor w/ 10 levels "A40","A41","A410",..: 5 5 8 4 1 8 4 2 5 1 ...

$ 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 ...

$ personal   : Factor w/ 4 levels "A91","A92","A93",..: 3 2 3 3 3 3 3 3 1 4 ...

$ debtors    : Factor w/ 3 levels "A101","A102",..: 1 1 1 3 1 1 1 1 1 1 ...

$ 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 ...

$ housing    : Factor w/ 3 levels "A151","A152",..: 2 2 2 3 3 3 2 1 2 2 ...

$ numcredits : int  2 1 1 1 2 1 1 1 1 2 ...

$ job        : Factor w/ 4 levels "A171","A172",..: 3 3 2 3 3 2 3 4 2 4 ...

$ 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 ...

> german$numcredits<-factor(german$numcredits)

> german$residence<-factor(german$residence)

> german$residpeople<-factor(german$residpeople)

> class(german$numcredits);class(german$residence);class(german$residpeople)

[1] "factor"

[1] "factor"

[1] "factor"

2) cart 방법 적용

> library(rpart)

> my.control<-rpart.control(xval=10,cp=0,minsplit=5)

> fit.german<-rpart(y~.,data=german,method='class',control=my.control)

> fit.german #최초의 나무. 가지치기를 하지 않은 최대 크기의 나무 보기

n= 1000

node), split, n, loss, yval, (yprob)

* denotes terminal node

1) root 1000 300 good (0.300000000 0.700000000) 

2) check=A11,A12 543 240 good (0.441988950 0.558011050) 

.

. (너무 커서 중략)

.

    253) credit>=1273 10   0 good (0.000000000 1.000000000) *

  127) check=A14 122   1 good (0.008196721 0.991803279) *

함수 설명

 1. rpart.control:

- xval=10: 교타 타당성의 fold 개수, 디폴트는 10

- cp=0: 오분류율이 cp값 이상으로 향상되지 않으면 더 이상 분할하지 않고 나무구조 생성을 멈춘다. cp값이 0이면 오분류값이 최소, 디폴트는 0.01

- minsplit=5: 한 노드를 분할하기 위해 필요한 데이터의 개수. 이 값보다 적은 수의 관측치가 있는 노드는 분할하지 않는다. 디폴트는 20

2. r.part

 - method=class: 나무 모형을 지정한다. anova는 회귀나무, poisson 포아송 회귀나무, class는 분류나무 exp는 생존나무. 디폴트는 class

 - na.action=na.rpart: 목표변수가 결측치이면 전체 관측치를 삭제. 입력변수가 결측치인 경우에는 삭제하지 않는다.

결과 해석

중간노드를 분할하는 최소 자료의 수를 5개로 지정하였고, cp값은 0으로 하여 나무모형의 오분류값이 최소가 될 때 까지 분할을 진행하였다. 또한 10-fold 교차타당성을 수행하여 최적의 cp값을 찾도록 하였다. 나무가 너무나 큰 관계로 중간 부분을 생략하였고, 용이한 모형 분석을 위해 가지치기를 해보자.

3) 큰 나무를 줄이기 위한 가지치기 작업

> printcp(fit.german)

Classification tree:

  rpart(formula = y ~ ., data = german, method = "class", control = my.control)

Variables actually used in tree construction:

  [1] age         check       credit      debtors     duration    employment  history   

[8] housing     installment job         numcredits  others      personal    property  

[15] purpose     residence   savings   

Root node error: 300/1000 = 0.3

n= 1000

               CP nsplit    rel error  xerror       xstd

1  0.0516667      0   1.00000 1.00000 0.048305

2  0.0466667      3   0.84000 0.94667 0.047533

3  0.0183333      4   0.79333 0.86333 0.046178

4  0.0166667      6   0.75667 0.87000 0.046294

5  0.0155556      8   0.72333 0.88667 0.046577

6  0.0116667     11   0.67667 0.88000 0.046464

7  0.0100000    13   0.65333 0.85667 0.046062

8  0.0083333     16   0.62333 0.87000 0.046294

9  0.0066667     18   0.60667 0.87333 0.046351

10 0.0060000     38   0.44333 0.92000 0.047120

11 0.0050000     43   0.41333 0.91000 0.046960

12 0.0044444     55   0.35333 0.92000 0.047120

13 0.0033333     59   0.33333 0.92000 0.047120

14 0.0029167     83   0.25000 0.97000 0.047879

15 0.0022222     93   0.22000 0.97667 0.047976

16 0.0016667     96   0.21333 0.97667 0.047976

17 0.0000000    104   0.20000 1.01333 0.048486

결과 해석

10-fold 교차타당성 방법에 의한 오분율(xerror)이 최소가 되는 값은 0.85667이며 이때의 cp값은 0.01임을 알 수 있다. 이 때 분리의 횟수가 13회(nsplit=13)인 나무를 의미한다.

또는 아래와 같은 방법으로도 최소 오분류값(xerror)를 찾을 수 있다.

> names(fit.german)

[1] "frame"               "where"               "call"                "terms"             

[5] "cptable"             "method"              "parms"               "control"           

[9] "functions"           "numresp"             "splits"              "csplit"            

[13] "variable.importance" "y"                   "ordered"           

> fit.german$cptable[,'xerror']

1         2         3         4         5         6         7         8         9

1.0000000 0.9466667 0.8633333 0.8700000 0.8866667 0.8800000 0.8566667 0.8700000 0.8733333

10        11        12        13        14        15        16        17

0.9200000 0.9100000 0.9200000 0.9200000 0.9700000 0.9766667 0.9766667 1.0133333

> which.min(fit.german$cptable[,'xerror'])

7

7

> fit.german$cptable[7,]

CP       nsplit     rel error       xerror        xstd

0.01000000 13.00000000  0.65333333  0.85666667  0.04606167

> fit.german$cptable[7]

[1] 0.01

> fit.german$cptable[which.min(fit.german$cptable[,'xerror'])]

[1] 0.01

> min.cp<-fit.german$cptable[which.min(fit.german$cptable[,'xerror'])]

> min.cp

[1] 0.01

> fit.prune.german<-prune(fit.german,cp=min.cp)

4) 오분류율이 최소인 cp값(=0.011)을 찾았으니 이 값을 기준으로 가지치기를 시행하자.

> fit.prune.german<-prune(fit.german,cp=0.01)

> fit.prune.german

결과 해석

node), split, n, loss, yval, (yprob) 기준으로 첫번째 결과를 분석하면 다음과 같다.

노드, 분할점, 개수, …공부 필요

16) duration>=47.5 36   5 bad (0.8611111 0.1388889) *

duration 변수 중 47.5보다 큰 경우, 전체 36(n)개를 bad(yval)로 분류하였고 그 중 5개(loss)가 good이다. 그리하여 bad로 분류되는 것은 31/36 = 0.8611111로 표기하게 되고, 5개의 loss는 5/36 = 1388889 로 그 확률을 볼 수 있다. 아래 plot에서는 bad 31/5로 표기

376) property=A123,A124 20   3 bad (0.8500000 0.1500000) *

377) property=A121,A122 45  14 good (0.3111111 0.6888889) *

property가 a123(car), a124(unknown / no property)의 경우 전체 20개를 bad로 분류하였고 3개의 loss 즉 good (3/20 = 0.15)로 분류하였다. 아래 plot에서는 bad 17/3로 표기 

property가 a121(real estate), a122(building society savings agreement)인 경우에는 전체 45개를 good으로 분류하였고 14개의 loss 즉 bad로 분류 (14/45=0.3111111), 아래 plot에서는 good 14/31로 표기

<< 17.6.18(일)>> 해석 부분 내용 추가

duration > = 22.5인 경우, 전체 고객은 237명이고, 이 중 신용도가 나쁜 사람의 비율은 56.5%이고 좋은 사람의 비율은43.5%로 103명이다. 따라서 duration > 22.5 그룹은 bad로 분류된다.

가지치기를 한 모형을 그림으로 나타내는 함수는 아래와 같다.

> plot(fit.prune.german,uniform = T,compress=T,margin=0.1)

> text(fit.prune.german,use.n=T,col='blue',cex=0.7)

왼쪽 가지의 가장 아랫부분의 분할점인 ‘purpose=acdeghj’는 purpose 변수의 범주값 중에서 알파벳 순서로, 1(=a), 3(=c), 4(=d), 5(=e), 7(=g), 8(=h), 10(=j)번째 범주값을 의미하며, fit.prune.german에서 각각 A40,A410,A42,A43,A45,A46,A49 임을 알 수 있다.

34) purpose=A40,A410,A42,A43,A45,A46,A49 137  52 bad (0.6204380 0.3795620) *

<< 17.6.18(일)>> 해석 부분 내용 추가

가장 우측의 duration > = 11.5가 아닌 경우, 신용다가 나쁜 / 좋은 사람의 비율은 9명 / . 4명이고, 신용도가 좋은 good으로 분류된다.

5) 나무수를 더 줄여보자.

> printcp(fit.german)

Classification tree:

  rpart(formula = y ~ ., data = german, method = "class", control = my.control)

Variables actually used in tree construction:

  [1] age         check       credit      debtors     duration    employment  history   

[8] housing     installment job         numcredits  others      personal    property  

[15] purpose     residence   savings   

Root node error: 300/1000 = 0.3

n= 1000

                CP nsplit   rel error  xerror      xstd

1  0.0516667      0   1.00000 1.00000 0.048305

2  0.0466667      3   0.84000 0.94667 0.047533

3  0.0183333      4   0.79333 0.86333 0.046178

4  0.0166667      6   0.75667 0.87000 0.046294

5  0.0155556      8   0.72333 0.88667 0.046577

6  0.0116667     11   0.67667 0.88000 0.046464

7  0.0100000     13   0.65333 0.85667 0.046062

8  0.0083333     16   0.62333 0.87000 0.046294

9  0.0066667     18   0.60667 0.87333 0.046351

10 0.0060000     38   0.44333 0.92000 0.047120

11 0.0050000     43   0.41333 0.91000 0.046960

12 0.0044444     55   0.35333 0.92000 0.047120

13 0.0033333     59   0.33333 0.92000 0.047120

14 0.0029167     83   0.25000 0.97000 0.047879

15 0.0022222     93   0.22000 0.97667 0.047976

16 0.0016667     96   0.21333 0.97667 0.047976

17 0.0000000    104   0.20000 1.01333 0.048486

5번째 단계이며 분리의 횟수가 8회(nsplit=8)인 나무는 교차타당성 오분류율이 0.88667로 최소는 아니지만 7번째 단계의 분리의 횟수 13회 나무 가지의 최소 오분류율 0.85667과는 크게 차이가 나지 않는다. 그리고 최소 오분류율 표준편차의 1배 범위(0.88667  <  0.85667 + 0.046062)에 있다. 이런 경우에는 5번째 단계이며 분리의 횟수가 8인 나무를 선택하는 경우도 있다.

5번째 단계이며 분리 횟수가 8인 cp값 0.0155556의 반올림 값 0.016 적용하여 다시 가지치기

> fit.prune.german<-prune(fit.german,cp=0.016)

> fit.prune.german

> plot(fit.prune.german,uniform=T,compress=T,margin=0.1)

> text(fit.prune.german,use.n=T,col='blue',cex=0.7)

6) 목표변수의 분류예측치를 구하고 그 정확도에 대해서 평가해 보자

 > fit.prune.german<-prune(fit.german,cp=0.01) 

 > pred.german=predict(fit.prune.german,newdata=german,type='class')

 > tab=table(german$y,pred.german,dnn=c('Actual','Predicted'))

 > tab

           Predicted

  Actual  bad  good

    bad  180   120

   good   76   624

함수 설명

predict(fit.prune.german,newdata=german,type='class'), type = class는 분류나무의 집단값 예측결과, 회귀나무라면 type = vector라고 해야 한다.

결과 해석

실제 good인데 good으로 예측한 것이 624개, 실제 bad인데 bad로 예측한 것이 180

따라서 오분류율은 {1000 – (624+180)} / 1000 = 19.6%

R코드를 이용하면 1-sum(diag(tab)) / sum(tab)

7) 마지막으로 독일신용평가데이터를 훈련데이터와 검증 데이터로 분할하여 분류나무를 평가해보자.

> set.seed(1234)

> i=sample(1:nrow(german),round(nrow(german)*0.7)) #70% for training훈련 data, 30% for test검증

> german.train=german[i,]

> german.test=german[-i,]

> fit.german<-rpart(y~.,data=german.train,method='class',control=my.control)

> printcp(fit.german)

Classification tree:

  rpart(formula = y ~ ., data = german.train, method = "class",

        control = my.control)

Variables actually used in tree construction:

  [1] age         check       credit      debtors     duration    employment  history   

[8] housing     installment job         numcredits  others      personal    property  

[15] purpose     residence   savings     telephone 

Root node error: 201/700 = 0.28714

n= 700

CP  nsplit   rel error  xerror     xstd

1  0.05721393      0   1.00000 1.00000 0.059553

2  0.03482587      2   0.88557 1.00498 0.059641

3  0.02985075      5   0.78109 1.00000 0.059553

4  0.01990050    6   0.75124 0.95025 0.058631

5  0.01741294      8   0.71144 0.96020 0.058822

6  0.01492537     10   0.67662 1.00000 0.059553

7  0.01243781     14   0.61692 1.00000 0.059553

8  0.00995025     17   0.57711 1.00995 0.059728

9  0.00746269     35   0.39303 1.03980 0.060238

10 0.00621891     46   0.30846 1.06965 0.060722

11 0.00497512     50   0.28358 1.04975 0.060402

12 0.00331675     58   0.24378 1.09950 0.061181

13 0.00248756     61   0.23383 1.11940 0.061474

14 0.00124378     69   0.21393 1.14925 0.061894

15 0.00099502     73   0.20896 1.14925 0.061894

16 0.00000000     78   0.20398 1.14925 0.061894

> fit.prune.german<-prune(fit.german,cp=0.02)

> fit.prune.german

> p.german.test=predict(fit.prune.german,newdata=german.test,type='class')

> tab=table(german.test$y,p.german.test,dnn=c('Actual','Predicted'))

> tab

            Predicted

Actual  bad good

      bad   34   65

      good  14  187

> 1-sum(diag(tab))/sum(tab) #오분류율

[1] 0.2633333 

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

6.1 회귀분석

예) 페인트의 불순도는 페인트를 얼마나 빨리 저어주느냐에 따라 달라진다. 아래표는 휘젓는 장치의 회전율과 불순도를 측정한 데이터이다.

> setwd('c:/Rwork')

> paint=read.csv('paint.csv')

> paint

x    y

1  20  8.4

2  22  9.5

3  24 11.8

4  26 10.4

5  28 13.3

6  30 14.8

7  32 13.2

8  34 14.7

9  36 16.4

10 38 16.5

11 40 18.9

12 42 18.5

> out=lm(y~x,data=paint)

> plot(y~x,data=paint)

> abline(out)

> out

Call:

  lm(formula = y ~ x, data = paint)

Coefficients:

  (Intercept)            x 

-0.2893       0.4566 

> summary(out)

Call:

  lm(formula = y ~ x, data = paint)

Residuals:

  Min      1Q  Median      3Q     Max

-1.1834 -0.5432 -0.3233  0.8333  1.3900

Coefficients:

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

(Intercept) -0.28928    1.22079  -0.237    0.817   

x            0.45664    0.03844  11.880 3.21e-07 ***

  ---

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

Residual standard error: 0.9193 on 10 degrees of freedom

Multiple R-squared:  0.9338,           Adjusted R-squared:  0.9272

F-statistic: 141.1 on 1 and 10 DF,  p-value: 3.211e-07

회귀식은 Y = -0.28928 + 0.45664X

기울기의 0.45664의 p-value는 3.21e-07로 유의.

결정계수R²는 0.9338로, 총 변동 중 93.38%가 회귀모형에 의해 설명되고 있다. 이는 두 변수 사이의 피어슨 상관계수의 제곱이다.

> with(paint,cor(x,y))^2

[1] 0.933832

 F값(df1=1,df2=10)은 141.1로 p-value는 3.21e-07로 유의.

cor.test()로 구한 Pearson 상관계수의 t(df=10)=11.88, p-value 값 3.211e-07은 위의 p-value값과 일치한다. T 값인 11.88의 제곱  11.88^2 값은 141.1344로 모형의 적합성을 나타내는 F값 141.1과 일치하고 당연히 p-value도 3.21e-07로 동일하다.

> with(paint,cor.test(x,y))

Pearson's product-moment correlation

data:  x and y

t = 11.88, df = 10, p-value = 3.211e-07

alternative hypothesis: true correlation is not equal to 0

95 percent confidence interval:

0.8810937 0.9907768

sample estimates:

cor

0.9663498

out으로 저장된 결과를 plot()을 이용하여 회귀진단에 필요한 그래프를 얻을 수 있다.

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

> plot(out)

해석:표준화된 잔차의 Normal Q-Q plot이 직선에 가깝고 다른 그래프들도 별다른 추세를 보이지 않는다.

resid()을 이용하면 회귀분석 결과물로부터 잔차를 어을 수 있다. 이렇게 받은 잔차에 아래의 명령어를 이용하여 normal Q-Q Plot을 그릴 수 있다.

> qqnorm(resid(out))

> qqline(resid(out))

정확한 p-value를 알고 싶다면 shapiro.test()를 이용하여 잔차가 정규분포를 따르는지 검정한다.

> shapiro.test(resid(out))

Shapiro-Wilk normality test

data:  resid(out)

W = 0.9196, p-value = 0.2826

H0 정규분포를 따른다

H1 정규분포를 따르지 않는다

p-value는 0.2826로 잔차는 정규분포 가정을 만족시킨다. 분산분석을 통하여도 회귀의 유의성을 검정할 수 있다.

> anova(out)

Analysis of Variance Table

Response: y

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

x          1 119.275 119.275  141.13 3.211e-07 ***

Residuals  10   8.451   0.845                     

---

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

6.2 공분산분석

실험이 잘 통제되어 종속변수Y의 변동을 설명하는 데 그룹변수 이외에 다른 변인이 없다면 Two-sample t-test나 일원분산분석을 하면 된다. 그러나 현실적으로 통제는 쉽지 않기 때문에 이 경우에는 여러가지 다른 변수들을 통제해주어 조사하고자 하는 변수만의 효과를 조사해야 한다. 공분산분석은 분산분석에 연속형 변수를 추가한 것이다. 궁극적인 목적인 각 그룹 간 평균들의 차이가 있는 검정하는 것으로 분산분석과 동일하나, 통제가 안 되는 연속형 변수(covariate)를 추가하여 오차를 줄이고 검정력을 높이는 것이 차이점이다.

lm(y~공변량변수 + 그룹변수)

> machine<-read.csv('machine.csv')

> head(machine)

machine  y  x

1      m1 36 20

2      m1 41 25

3      m1 39 24

4      m1 42 25

5      m1 49 32

6      m2 40 22

> levels(machine$machine)

[1] "m1" "m2" "m3"

해석: m1기계1이 대조군이고 m2와 m3은 비교군이다.

> out=lm(y~x+machine,data=machine)

> anova(out)

Analysis of Variance Table

Response: y

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

x          1 305.130 305.130 119.9330 2.96e-07 ***

machine    2  13.284   6.642   2.6106   0.1181   

Residuals 11  27.986   2.544                     

---

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

 해석: p-value 값 0.1181로 기계 간의 차이가 있다고 결론을 내릴 수 없다. 

summary()로 자세한 결과물을 보자.

> summary(out)

Call:

  lm(formula = y ~ x + machine, data = machine)

Residuals:

  Min      1Q  Median      3Q     Max

-2.0160 -0.9586 -0.3841  0.9518  2.8920

Coefficients:

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

(Intercept)   17.360      2.961   5.862 0.000109 ***

x              0.954      0.114   8.365 4.26e-06 ***

machinem2      1.037      1.013   1.024 0.328012   

machinem3     -1.584      1.107  -1.431 0.180292   

---

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

Residual standard error: 1.595 on 11 degrees of freedom

Multiple R-squared:  0.9192,        Adjusted R-squared:  0.8972

F-statistic: 41.72 on 3 and 11 DF,  p-value: 2.665e-06

해석: m2를 reference인 m1와 비교한 machinem2의 p-value값이 0.328012이고 m3을 m1과 비교한 machinem3의 p-value는 0.180292이다. 그러나 이 p-value들을 조정없이 사용하여 두 machine 모두 m1와 유의한 차이가 없다고 결론을 내리지는 않는다.

1개의 대조군(m1)을 2개의 비교군과 비교하므로 Dunnett의 방법으로 p-value를 조정한다.

> library(multcomp)

> dunnett=glht(out,linfct=mcp(machine='Dunnett'))

> summary(dunnett)

Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Dunnett Contrasts

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

Linear Hypotheses:

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

m2 - m1 == 0    1.037      1.013   1.024    0.518

m3 - m1 == 0   -1.584      1.107  -1.431    0.304

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

> plot(dunnett)

 해석: 두 p-value 모두가 원래 p-value보다 커졌으며, 여전히 유의한 차이가 있다고 결론을 지을 수 없다. Plot 그래프로 확인할 수도 있다. 신뢰구간이 0을 포함하고 있으므로 유의하다고 볼 수 없다.

> with(data=machine,tapply(y,machine,mean))

m1   m2   m3

41.4 43.2 36.0

> with(data=machine,tapply(y,machine,sd))

m1       m2       m3

4.827007 3.701351 3.807887

해석: 각 기계에서 생산된 섬유 제품의 강도 평균은 41.4, 43.2, 36.0으로 얻어져 혹시 차이가 있다면 기계3의 강도가 떨어지는 듯 하다.

출처: 실험 계획과 응용, R로 하는 통계 분석

이원산분석은 그룹변수가 하나인 일원산분석의 확장으로 두 개의 그룹변수를 가진다. 이 분석의 특징은 두 그룹변수들의 효과뿐만 아니라 두 그룹 변수들이 서로 어떤 영향을 미치는지 교호작용을 볼 수 있다는 것이다

1) 교호 작용이 있는 경우: lm(y~변수a*변수b, data=)

2) 교호 작용이 없는 경우: summary(lm(y~변수a+변수b, data = )

예시) 각기 다른 사료 a1~a4와 돼지품종 b1~3이 체중 증가에 미치는 영향을 조사

> setwd('c:/Rwork')

> pigs<-read.csv('pigs.csv')

> head(pigs)

result breed feed

1     64    b1   a1

2     66    b1   a1

3     70    b1   a1

4     65    b1   a2

5     63    b1   a2

6     58    b1   a2

> with(pigs,tapply(result,breed,mean))

b1       b2       b3

60.25000 61.33333 54.25000

> with(pigs,tapply(result,feed,mean))

a1       a2       a3       a4

67.44444 56.66667 58.55556 51.77778

> with(pigs,tapply(result,list(breed,feed),mean))

a1       a2       a3       a4

b1 66.66667 62.00000 64.00000 48.33333

b2 72.33333 50.66667 65.33333 57.00000

b3 63.33333 57.33333 46.33333 50.00000

> boxplot(result~breed+feed,col='blue',data=pigs)

기술통계를 위해 tapply 함수를 활용하여 품종, 사료, 품종x사료 조합의 평균을 구하였고, boxplot 함수를 통해 각 변수 조합의 result 분포 차이를 나타냈다.

교호 작용 검사

> with(pigs,interaction.plot(x.factor=feed,trace.factor=breed,response=result,fun=mean,type='b',main='Interaction Plot',pch=c(1,2,16)))

#교호 작용이 있는 경우

> out=lm(result~breed*feed,data=pigs)

> out

> anova(out)

Analysis of Variance Table

Response: result

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

breed       2  349.39  174.69  2.7926 0.081214 .

feed        3 1156.56  385.52  6.1628 0.002939 **

breed:feed   6  771.28  128.55  2.0549 0.097118 .

Residuals  24 1501.33   62.56                   

---

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

품종(breed)는 유의하지 않지만 p-value가 0.05에 근접하고, 사료(feed)는 유의하다. 품종에 따라 사료가 주는 효과의 크기가 조금씩 다를 수 있지만, Interaction인 breed:feed의 p-value는 0.09로 유의하지 않으므로 교호작용이 있다는 결론을 내릴 수 없다. 유의성 검정 결과를 보면 사료의 종류에 따라 체중증가에 유의한 차이가 있으나, 돼지 품종에 따라서는 체중증가에 차이가 있다는 결론은 내릴 수 없다. 품종에 대한 유의확률은 0.081로 약간의 영향력을 행사하는 듯하다.

#교호 작용이 없는 경우

 to be updated..

출처: 실험 계획과 응용, R로 하는 통계 분석

제3장 일원배치 완전확률화 계획법,

예제) 직물의 긁힘에 대한 저항력을 측정하기 위해 원단 A1,2,3,4의 제품 마모도를 비교하조가 한다.

> setwd('c:/Rwork')

> fabric = read.csv('fabric.csv')

> head(fabric)

result group

1   1.93    a1

2   2.38    a1

3   2.20    a1

4   2.25    a1

5   2.55    a2

6   2.72    a2

> boxplot(result~group,data=fabric)

> with(fabric,tapply(result,group,mean))

a1   a2   a3   a4

2.19 2.68 2.42 2.31

> round(with(fabric,tapply(result,group,mean)))

a1 a2 a3 a4

2  3  2  2

> with(fabric,tapply(result,group,sd))

a1         a2         a3         a4

0.18920888 0.08906926 0.18036999 0.05715476

> out=lm(result~group,data=fabric)

> out

Call:

  lm(formula = result ~ group, data = fabric)

Coefficients:

  (Intercept)      groupa2      groupa3      groupa4 

2.19         0.49         0.23         0.12 

Lm의 결과물을 보면 intercept 추정치가 a1의 평균과 같으므로 a1이 reference로 사용되었음을알 수 있다. A1와 a2 평균이 각각 2.19, 2.68으로 그 차이가 groupa2 임을 알 수 있다.

> summary(out)

Call:

  lm(formula = result ~ group, data = fabric)

Residuals:

  Min      1Q  Median      3Q     Max

-0.2600 -0.0700  0.0150  0.0625  0.2600

Coefficients:

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

(Intercept)  2.19000    0.07050  31.062 7.79e-13 ***

  groupa2      0.49000    0.09971   4.914 0.000357 ***

  groupa3      0.23000    0.09971   2.307 0.039710 * 

  groupa4      0.12000    0.09971   1.204 0.251982   

---

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

Residual standard error: 0.141 on 12 degrees of freedom

Multiple R-squared:  0.6871,           Adjusted R-squared:  0.6089

F-statistic: 8.785 on 3 and 12 DF,  p-value: 0.002353

> anova(out)

Analysis of Variance Table

Response: result

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

group     3 0.5240 0.174667  8.7846 0.002353 **

Residuals 12 0.2386 0.019883                   

---

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

분산분석표에서 group의 F(df1=3,df2=12)=8.7846 (p-value=0.002353)으로 세 그룹의 평균들이 모두 같다는 귀무가설을 기각한다.

귀무가설 H0: μa1 = μa2 = μa3 = μa4

대립가설 H1: not H0

검정통계량 F(df1=3,df2=12)=8.7846

p-value = 0.002353 < 0.05

결정 : 모든 그룹의 평균이 같다는 귀무가설을 기각한다. 4개의 그룹 중 어느 쌍에 차이가 있는지 다중비교로 좀 더 진단할 수 있다.

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

> plot(out)

> abline(out,col='blue')

Normal Q-Q Plot이 직선에 가깝고 다른 그림들도 별다른 추세를 보이지 않으므로 문제가 없다는 결론을 내린다.

> shapiro.test(resid(out))

Shapiro-Wilk normality test

data:  resid(out)

W = 0.97318, p-value = 0.8872

Shapiro.test()를 이용하여 잔차의 정규성검정 결과는 p-value=0.8872로 정규분포 가정을 만족시킨다.

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

> library(multcomp)

> dunnett=glht(out,linfct=mcp(group='Dunnett'))

> dunnett

General Linear Hypotheses

Multiple Comparisons of Means: Dunnett Contrasts

Linear Hypotheses:

  Estimate

a2 - a1 == 0     0.49

a3 - a1 == 0     0.23

a4 - a1 == 0     0.12

> summary(dunnett)

Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Dunnett Contrasts

Fit: lm(formula = result ~ group, data = fabric)

Linear Hypotheses:

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

a2 - a1 == 0  0.49000    0.09971   4.914   <0.001 ***

a3 - a1 == 0  0.23000    0.09971   2.307   0.0962 . 

a4 - a1 == 0  0.12000    0.09971   1.204   0.5085   

---

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

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

> plot(dunnett)

a2-a1의 p-value만 0.001로 유의

> tukey=glht(out,linfct=mcp(group='Tukey'))

> tukey

General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts

Linear Hypotheses:

  Estimate

a2 - a1 == 0     0.49

a3 - a1 == 0     0.23

a4 - a1 == 0     0.12

a3 - a2 == 0    -0.26

a4 - a2 == 0    -0.37

a4 - a3 == 0    -0.11

> summary(tukey)

Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts

Fit: lm(formula = result ~ group, data = fabric)

Linear Hypotheses:

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

a2 - a1 == 0  0.49000    0.09971   4.914  0.00187 **

a3 - a1 == 0  0.23000    0.09971   2.307  0.15096  

a4 - a1 == 0  0.12000    0.09971   1.204  0.63632  

a3 - a2 == 0 -0.26000    0.09971  -2.608  0.09215 .

a4 - a2 == 0 -0.37000    0.09971  -3.711  0.01375 *

a4 - a3 == 0 -0.11000    0.09971  -1.103  0.69444  

---

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

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

> plot(tukey)

a2-a1, a4-a2간 95% 신뢰구간이 0을 포함하지 않으므로 유의한 차이가 있음을 알 수 있다.

분산분석은 모든 그룹의 평균이 같은지를 유의수준 0.05로 한 번에 검정한다. 그러나 여러 쌍의 평균이 같은지 다시 검정하면 비교하는 쌍이 하나 이상으로 늘어나 실제로 유의하지 않은데 유의하게 나올 가능성이 높아지므로 유의수준을 0.05보다 낮추거나, p-value를 높게 조종해야 한다. 아래 표처럼 실제 Tukey의 p-value가 Dunnett의 p-value보다 높게 조정되었음을 알 수 있다.

출처: 실험 계획과 응용, R로 하는 통계 분석

2.3 짝지어진 비교 Paired T-Test

Paired T-Test는 쌍을 이룬 두 변수의 차이를 보는 검정이다. 가장 흔한 예는 한 집단을 대상으로 어떤 개입의 효과를 보기 위해 개입 전-후 값을 비교하여 개입의 효과를 측정하는 것이다.

Paired t-test는 다음 순서를 따른다.

 1) 정규분포를 따르는지 검정 with(데이터, shapiro.test(사후-사전)

 2) 정규분포를 따르면 사후-사전 차이의 평균이 0인지 검정하는 one sample t-test를 이용

   with(데이터, t.test(사후-사전))

 3) 정규분포를 따르지 않으면 비모수적 방법을 적용 with(데이터, wilcox.test(사후-사전)

예) 운동화 밑창에 사용되는 B 재질과 A 재질 간 마모도에 차이가 있는지 알아보고자 한다. 이를 위하여 10명의 아이들을 대상으로 임의로 한 쪽 발에는 재질 A의 밑창을 단 운동화를 한 쪽 발에는 재질 B의 밑창을 단 운동화를 신겨 일정기간 사용하도록 한 다음 밑창이 얼마나 마모되었는지를 측정한 결과 다음 자료를 얻었다.

> shoes<-read.csv('shoes.csv',header=T)

> head(shoes)

A    B

1 13.2 14.0

2  8.2  8.8

3 10.9 11.2

4 14.3 14.2

5 10.7 11.8

6  6.6  6.4

#paired t-test 1. 정규분포를 따르는지 검정

> with(shoes,shapiro.test(B-A))

Shapiro-Wilk normality test

data:  B - A

W = 0.96132, p-value = 0.8009

귀무가설 H0: B-A의 차이가 정규분포를 따른다.

대립가설 H1: B-A의 차이가 정규분포를 따르지 않는다.

검정통계량 W = 0.96132

p-value = 0.8009

결정 : 귀무가설 기각할 수 없다. B-A의 차이가 정규 분포를 따른다.

#paired t-test 2. 정규분포를 따르면 one sample t-test를 이용

> with(shoes,t.test(B-A))

One Sample t-test

data:  B - A

t = 3.3489, df = 9, p-value = 0.008539

alternative hypothesis: true mean is not equal to 0

95 percent confidence interval:

  0.1330461 0.6869539

sample estimates:

  mean of x

0.41

귀무가설 H0: μ = 0

대립가설 H1: μ ≠ 0

검정통계량 t(df=9) = 3.3489

p-value = 0.008539

결정 : 귀무가설 기각한다. 두 신발 재질의 마모도에 유의한 차이가 있다.

#paired t-test 3. 정규분포를 따르지 않는 경우 wilcox.test를 사용한다.

> with(shoes,wilcox.test(B-A)

귀무가설 H0: 전후 차이의 median은 0이다

대립가설 H1: 전후 차이의 median은 0이 아니다

검정통계량 V =

p-value =

결정: if p-value > 0.05 then 귀무가설 기각할 수 없다, 전후 차이가 없다.

출처: 실험계획과 응용, R로 하는 통계 분석

2.2 독립표본을 이용한 두 모평균 차이에 대한 추론. Two-sample T-test

독립표본을 바탕으로 두 개의 모집단의 평균을 비교. 가장 흔한 실험 연구는 실험군과 대조군에 서로 다른 개입(intervention)을 적용시킨 후 두 집단의 평균이 같은지를 비교하여 개입 효과의 차이를 평가하는 것이다. 이 경우 two-sample t-test를 사용하는데, 서로 독립적인 두 변수 간에 차이의 평균이 0인지를 검정한다.

Two-sample t-test는 다음 순서를 따른다.

 1) 두 집단의 분산이 같은지 검정한다. var.test(y~그룹변수)

 2) 분산이 다르면 Welch의 t-test를 적용한다. t.test(y~그룹변수)

 3) 분산이 같으면 pooled variance를 이용한 t-test를 적용한다. t.test(y~그룹변수, var.equal=TRUE)

예제) 제약회사에서 어떤 약을 오래 보관해도 약효가 지속되는지를 검사하려고 한다. 표본1과 2를 랜덤추출한 결과가 아래와 같다.

> medical<-read.csv('medical.csv',header=T)

> head(medical,3)

sample result

1 sample1   10.2

2 sample1   10.5

3 sample1   10.3

> tail(medical,3)

sample result

18 sample2    9.6

19 sample2    9.8

20 sample2    9.9

> boxplot(result~sample,data=medical)

-> 해석: Sample1의 분산이 sample2보다 더 큼을 알 수 있다.

#two sample test 1. 등분산 검정

> var.test(result~sample,data=medical)

F test to compare two variances

data:  result by sample

F = 1.7965, num df = 9, denom df = 9, p-value = 0.3959

alternative hypothesis: true ratio of variances is not equal to 1

95 percent confidence interval:

  0.4462364 7.2328801

sample estimates:

  ratio of variances

1.796545  

귀무가설 H0: σ²1 = σ²2

대립가설 H1: σ²1 ≠ σ²2

검정통계량 F(df1=9, df2=9) = 1.7965

p-value = 0.3959

결정 : 귀무가설 기각할 수 없다. 등분산을 가정한다.

> 1/1.7965  #F값 1.7965의 역수는 0.55로 sample2의 분산이 sample1대비 0.55배임을 알 수 있다.

[1] 0.5566379

#two sample test 2. 분산이 같은 경우, pooled variance사용

> t.test(result~sample,var.equal=TRUE,data=medical)

Two Sample t-test

data:  result by sample

t = 3.8511, df = 18, p-value = 0.00117

alternative hypothesis: true difference in means is not equal to 0

95 percent confidence interval:

  0.222688 0.757312

sample estimates:

  mean in group sample1 mean in group sample2

10.32                  9.83

귀무가설 H0: μ1 = μ2

대립가설 H1: μ1 ≠ μ2

검정통계량 t(df=18) = 3.8511

p-value = 0.00117

결정 : 귀무가설 기각한다. 두 모집단의 평균이 다르다.

#two sample test 3. 분산이 다른 경우, Welch의 t-test를 한다.

> t.test(result~sample,data=medical)

출처: 실험계획과 응용, R로 하는 통계 분석

2.4 표준화된 중회귀분석

표준화된 중회귀모형

Y_i^* = α₁Z_i1 + α₂Z_i2 + …+ α_kZ_ik + ε_i’

가 되며, 표준화된 중회귀모형에서 추정된 회귀계수 α_i의 절대값이 크면 클수록 설명변수 X_i가 반응변수 Y_i에 주는 영향이 크게 됨.

> library(lm.beta)

> market2.lm<-lm(Y~X1+X2,data=market2)

> market2.beta<-lm.beta(market2.lm)

> market2.beta

Call:

lm(formula = Y ~ X1 + X2, data = market2)

Standardized Coefficients::

(Intercept)          X1          X2

  0.0000000   0.7015566   0.3376137

> summary(market2.beta)

Call:

lm(formula = Y ~ X1 + X2, data = market2)

Residuals:

     Min       1Q   Median       3Q      Max

-1.30465 -0.69561 -0.01755  0.69003  1.53127

Coefficients:

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

(Intercept)  0.85041      0.00000    0.84624   1.005 0.334770   

X1         1.55811      0.70156    0.14793  10.532 2.04e-07 ***

X2         0.42736      0.33761    0.08431   5.069 0.000276 ***

---

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

Residual standard error: 0.9318 on 12 degrees of freedom

Multiple R-squared:  0.9799,    Adjusted R-squared:  0.9765

F-statistic: 292.5 on 2 and 12 DF,  p-value: 6.597e-11

-> R 결과에서 변수 X1, X2에 대한 표준화 계수는 0.7016, 0.3376이 됨을 알 수 있다. 따라서 적합된 표준화된계수 모형은

Ŷ* = 0.7016Z₁ + 0.3376Z₂

가 된다.여기서 X1의 계수가 X2의 계수보다 크므로 상대적으로 X1의 영향이 더 큼을 알 수 있다.

2.5 추정과 검정

Market2 데이터에 대하여

(1) X1=10, X2=10에 E(Y)=95% , 99%를 신뢰구간으로 추정하고

(2) H₀ : β­₁ = 0, H₀ : β­₂ = 0에 대하여 유의수준 α =0.05로 가설검정을 해 보자

#95%의 신뢰구간

> pred.x<-data.frame(X1=10,X2=10)

> pred.x

  X1 X2

1 10 10

> pc=predict(market2.lm,int='c',newdata=pred.x)

> pc

       fit      lwr      upr

1 20.70503 19.95796 21.45209

> round(pc,3)

     fit    lwr    upr

1 20.705 19.958 21.452

-> X1=10, X2=10 에서의 추정값은 20.705이고 95%의 신뢰구간은 (19.958,21.452)가 된다.

#99%의 신뢰구간

> pc99<-predict(market2.lm,int='c',level=0.99,newdata=pred.x)

> pc99

       fit      lwr      upr

1 20.70503 19.65769 21.75236

> round(pc99,3)

     fit    lwr    upr

1 20.705 19.658 21.752

-> X1=10, X2=10 에서의 추정값은 20.705이고 99%의 신뢰구간은 (19.658,21.752)로 95%의 신뢰구간은 (19.958,21.452)보다 더 넓게 형성된다.

H₀ : β­₁ = 0, H₀ : β­₂ = 0에 대한 가설 검정은 회귀적합 결과를 이용하면 된다.

> summary(market2.lm)

Call:

lm(formula = Y ~ X1 + X2, data = market2)

Residuals:

     Min       1Q   Median       3Q      Max

-1.30465 -0.69561 -0.01755  0.69003  1.53127

Coefficients:

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

(Intercept)  0.85041    0.84624   1.005 0.334770   

X1         1.55811    0.14793  10.532 2.04e-07 ***

X2         0.42736    0.08431   5.069 0.000276 ***

---

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

Residual standard error: 0.9318 on 12 degrees of freedom

Multiple R-squared:  0.9799,    Adjusted R-squared:  0.9765

F-statistic: 292.5 on 2 and 12 DF,  p-value: 6.597e-11

Beta hat 1 값은 1.55811 이고, 표준오차는 0.14793이므로

t-값은 1.55811 / 0.14793 = 10.53275

유의확률 p-값 2.04 X 10⁷ 이 되므로 귀무가설 H₀ : β­₁ = 0은 기각한다. R에서 ***는 p-값이 0.001보다 작은 경우로 표시되었음을 알 수 있다. H₀ : β­₂ = 0 의 가설 역시 p-값이 0.000276이므로 귀각한다.

2.6 변수 추가

중회귀모형을 적합할 때 어떤 특정한 변수를 회귀모형에 포함시키는 것이 바람직한가를 결정하고 싶은 경우가 있다. 이러한 경우 이 변수를 포함시키지 않고 구한 회귀제곱함(SSR)에서 이 변수를 포함시키고 구한 회귀제곱함(SSR)이 추가적으로 어느 정도 커졌는가를 검토하는것이 좋을 것이다. 이와 같은 경우에 추가적으로 증가된 제곱함을 추가제곱합(extra sum of squares)라고 부른다.

Health data 불러오기

X1: 몸무게(파운드),

X2: 분당 맥박수,

X3: 근력(들어올릴 수 있는 무게: 파운드),

X4: 1/4 마일 달리는 시간(초),

Y: 1마일 달리는 시간(초)

> health<-read.table('health.csv',header=T)

> head(health,2)

   번호.X1.X2.X3.X4.Y

1 1,217,67,260,91,481

2 2,141,52,190,66,292

> health<-read.table('health.csv',header=T,sep=',')

> head(health,2)

  번호  X1 X2  X3 X4   Y

1    1 217 67 260 91 481

2    2 141 52 190 66 292

> colnames(health[1])

[1] "번호"

> colnames(health)[1]<-'ID'

> colnames(health)[1]

[1] "ID"

> head(health)

  ID  X1 X2  X3 X4   Y

1  1 217 67 260 91 481

2  2 141 52 190 66 292

회귀분석 실시

> h4.lm=lm(Y~X1+X2+X3+X4,data=health)

> summary(h4.lm)

Call:

lm(formula = Y ~ X1 + X2 + X3 + X4, data = health)

Residuals:

   Min     1Q Median     3Q    Max

-54.49 -21.91  -5.10  18.66  45.69

Coefficients:

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

(Intercept)  -6.8606    59.0070  -0.116   0.9084   

X1            1.3826     0.2933   4.713 7.84e-05 ***

X2           -0.3745     0.8955  -0.418   0.6794   

X3           -0.5302     0.2571  -2.062   0.0497 * 

X4            3.6202     0.7573   4.781 6.58e-05 ***

---

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

Residual standard error: 29.96 on 25 degrees of freedom

Multiple R-squared:  0.8396,    Adjusted R-squared:  0.814

F-statistic: 32.72 on 4 and 25 DF,  p-value: 1.332e-09

> anova(h4.lm)

Analysis of Variance Table

Response: Y

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

X1         1  89117   89117 99.2926 3.444e-10 ***

X2         1   4680    4680  5.2142   0.03117 * 

X3         1   3165    3165  3.5260   0.07213 . 

X4         1  20513   20513 22.8548 6.578e-05 ***

Residuals 25  22438     898                     

---

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

모형 적합 후 anova 함수 사용

> h1.lm=lm(Y~X1,data=health)

> h2.lm=lm(Y~X1+X4,data=health)

> h3.lm=lm(Y~X1+X3+X4,data=health)

> h4.lm=lm(Y~X1+X2+X3+X4,data=health)

> anova(h1.lm,h2.lm)

Analysis of Variance Table

Model 1: Y ~ X1

Model 2: Y ~ X1 + X4

  Res.Df   RSS Df Sum of Sq      F    Pr(>F)    

1     28 50795                                 

2     27 26419  1     24376 24.912 3.119e-05 ***

---

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

-> 모형1의 잔차제곱합 50795, 모형2의 잔차제곱함26419, 추가제곱합 50795-26419 = 24376 이고, 검정통계랑 F₀는 24.912가 된다. 이에 대한 유의확률 p-값은 3.119 X 10⁵이므로 변수4가 유의한 변수임을 알 수 있다.

> anova(h2.lm,h3.lm)

Analysis of Variance Table

Model 1: Y ~ X1 + X4

Model 2: Y ~ X1 + X3 + X4

  Res.Df   RSS Df Sum of Sq      F  Pr(>F) 

1     27 26419                             

2     26 22595  1    3824.4 4.4007 0.04579 *

---

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

> anova(h3.lm,h4.lm)

Analysis of Variance Table

Model 1: Y ~ X1 + X3 + X4

Model 2: Y ~ X1 + X2 + X3 + X4

  Res.Df   RSS Df Sum of Sq      F Pr(>F)

1     26 22595                          

2     25 22438  1    156.93 0.1748 0.6794

-> (X1,X4)인 모형에 X3을 추가하는 경우의 추가제곱합 3842.4 p-값은 0.04579로 X3은 유의한 변수임을 알 수 있고, (X1,X3,X4)인 모형에 X2을 추가하는 경우는 추가제곱합 156.93, p-값은 0.6794로 X2은 유의미한 변수가 아님을 알 수 있다.

추가변수그림

새로운 변수의 효과를 그래프로 표현할 수 있는데, 이러한 그래 중의 하나가 추가변수그림(added variable plot)이다.

> library(car)

> h4.lm=lm(Y~X1+X2+X3+X4,data=health)

> summary(h4.lm)

Call:

lm(formula = Y ~ X1 + X2 + X3 + X4, data = health)

Residuals:

   Min     1Q Median     3Q    Max

-54.49 -21.91  -5.10  18.66  45.69

Coefficients:

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

(Intercept)  -6.8606    59.0070  -0.116   0.9084   

X1            1.3826     0.2933   4.713 7.84e-05 ***

X2           -0.3745     0.8955  -0.418   0.6794   

X3           -0.5302     0.2571  -2.062   0.0497 * 

X4            3.6202     0.7573   4.781 6.58e-05 ***

---

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

Residual standard error: 29.96 on 25 degrees of freedom

Multiple R-squared:  0.8396,    Adjusted R-squared:  0.814

F-statistic: 32.72 on 4 and 25 DF,  p-value: 1.332e-09

> avPlots(h4.lm)

-> 해석: 중회귀모형에서 어떤 특정한 회귀모형에 포함시키조가 할 떄, 변수 선택은 기존의 모형이 설명하지 못하는 부분을 새로운 변수가 들어옴으로써 추가설명력이 얼마나 유의한가에 따라 결정된다. 변수 X1의 추가변수그림에서 회귀계수는 1.3826, X4는 3.6202, 변수X1와 X에 대하 추가변수 그림이 선형성이 강한 것을 볼 수 이씅며, 따라서 이 두 변수가 회귀모형에 매우 유의하다는 것을 알 수있다. 

목표변수가 연속형인 경우 -> 선형 회귀모델, ex) 광고비 투입 대비 매출액

목표변수가 두 개의 범주를 가진 이항형인 경우 -> 로지스틱 회귀모형, ex) 좋다1, 나쁘다0

독일신용평가 데이터 셋

독일신용평가 데이터(German Credit Data)는 머신러닝 저장소에 탑재되어 있는 데이터로 분류의 예제로 많이 활용된다.

1. 데이터 불러오기

2. 로지스틱 회귀 분석 시작

->  <!--[endif]-->해석: fit.step = step(fit.all, direction='both'를 통해 AIC가 가장 작은 모형을 찾는다.

check는 4개의 범주(checkA11 계좌 없음 / A12 잔액 없음 / A13 잔액 200 이하 / A14 잔액 200 이상)를 가지므로 3개의 가변 수 생성. 추정된 회귀계수는 모두 양수이므로, A12~A14 즉 계좌가 있는 경우 계좌 없음(A11)대비 신용이 좋을 확률(Y=1)이 더 높다. 대출기간인 duration은 마이너스의 값을 지니므로 대출 기간이 오래 될 수록 신용도는 낮아진다. 모델의 AIC는 980.91로, AIC가 클 경우 그 모형은 적합하지 않기 때문에, 여러 후보 모형 중에서 AIC가 가장 작은 모형을 선택한다.

3. 예측함수 및 정오분류표 작성

->  해석: 1로 예측할 확률이 임계치(threshold) 0.5보다 클 경우에는 1로, 0.5이하일 경우에는 0으로 예측. 실제로는 0인데 0으로 예측한 경우가 158개, 1인데 1로 분류한 경우가 618개이다.반면에 0인데 1로 오분류한 경우가 142개, 1인데 0으로 오분류한 경우가 82개이다.

4. 예측력 측도

-> 해석: 민감도는 실제 양성(Y=1)일 때 양성으로 예측할 확률, 특이도는 실제 음성(Y=0)일 때 음성으로 예측할 확률이다. 예측정확도(prediction accuracy)는 실제 양서일 때 양성으로, 음성일 때 음성으로 제대로 예측할 확률로 민감도와 특이도의 가중평균이다. 오분류율(misclassification rate)는 양성일 때 음성으로, 음성일 때 양성으로 잘못 예측할 확률이다.

5. ROC 곡선 및 AUC 생성

->  민감도와 특이도는 임계치에 다라 달라지고 임계치는 상황에 따라 다르게 결정할 수 이다. 여러 가능한 임계치에 대해 ‘1-특이도(Specificity)’를 가로축에, 민감도를 세로축에 놓고 그린 그래프를 ROC(Receiver operating characteristic) 곡선이라 한다. 민감도와 특이도가 높을수록 예측력이 좋다고 할 수 있기 때문에 ROC 곡선이 좌상단에 가까울수록 ROC 곡선 아래 면적인 AUC(area under the ROC curve)가 커지고, 예측력이 좋다고 할 수 있다.이 독일신용평가 데이터에 적합한 로지스틱 회귀모형에 대한 예측력의 측도인 AUC는 최대치 1보다 다소 작은 0.831로 상당히 높음을 알 수 있다.

목표변수가 연속형인 경우 -> 선형 회귀모델, ex) 광고비 투입 대비 매출액

목표변수가 두 개의 범주를 가진 이항형인 경우 -> 로지스틱 회귀모형, ex) 좋다1, 나쁘다0

보스턴 하우징 데이터 Housing Values in Suburbs of Boston

(출처: http://127.0.0.1:31866/library/MASS/html/Boston.html)

1. 데이터 불러오기

2. 선형 회귀 모형 만들기