데이터마이너를 꿈꾸며

 회귀모형의 적합결과 및 분산분석

> setwd("~/Documents/Rwork/KNOU/2017.2학기/R을 이용한 다변량 분석")

> hald.data = read.table('hald.txt',header=T)

> hald.data

       y x1 x2 x3 x4

1   78.5  7 26  6 60

2   74.3  1 29 15 52

3  104.3 11 56  8 20

4   87.6 11 31  8 47

5   95.9  7 52  6 33

6  109.2 11 55  9 22

7  102.7  3 71 17  6

8   72.5  1 31 22 44

9   93.1  2 54 18 22

10 115.9 21 47  4 26

11  83.8  1 40 23 34

12 113.3 11 66  9 12

13 109.4 10 68  8 12

> hald.ls = lm(y~.,data=hald.data)

> hald.ls

Call:

lm(formula = y ~ ., data = hald.data)

Coefficients:

(Intercept)           x1           x2           x3           x4  

    62.4054       1.5511       0.5102       0.1019      -0.1441  

> summary(hald.ls)

Call:

lm(formula = y ~ ., data = hald.data)

Residuals:

    Min      1Q  Median      3Q     Max 

-3.1750 -1.6709  0.2508  1.3783  3.9254 

Coefficients:

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

(Intercept)  62.4054    70.0710   0.891   0.3991  

x1            1.5511     0.7448   2.083   0.0708 .

x2            0.5102     0.7238   0.705   0.5009  

x3            0.1019     0.7547   0.135   0.8959  

x4           -0.1441     0.7091  -0.203   0.8441  

---

Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.446 on 8 degrees of freedom

Multiple R-squared:  0.9824, Adjusted R-squared:  0.9736 

F-statistic: 111.5 on 4 and 8 DF,  p-value: 4.756e-07

> anova(hald.ls)

Analysis of Variance Table

Response: y

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

x1         1 1450.08 1450.08 242.3679 2.888e-07 ***

x2         1 1207.78 1207.78 201.8705 5.863e-07 ***

x3         1    9.79    9.79   1.6370    0.2366    

x4         1    0.25    0.25   0.0413    0.8441    

Residuals  8   47.86    5.98                       

---

데이터 생성

> a = c(1:10)

> a

 [1]  1  2  3  4  5  6  7  8  9 10

> a = c(20:10)

> a

 [1] 20 19 18 17 16 15 14 13 12 11 10

> a %/% 3

 [1] 6 6 6 5 5 5 4 4 4 3 3

> a %% 3

 [1] 2 1 0 2 1 0 2 1 0 2 1

> a = c(1:3) == c(3:1)

> a

[1] FALSE  TRUE FALSE

SEQ 함수

> x = -5:5

> x

 [1] -5 -4 -3 -2 -1  0  1  2  3  4  5

> x[x<0]

[1] -5 -4 -3 -2 -1

> mean(x)

[1] 0

> x[x<0]=0

> mean(x)

[1] 1.363636

> consec=(1:10)*10

> consec

 [1]  10  20  30  40  50  60  70  80  90 100

> consec[2]

[1] 20

> consec[c(3,6,8,10)]

[1]  30  60  80 100

> consec[c(-2,-4)]

[1]  10  30  50  60  70  80  90 100

> x=c(10:1)

> consec<60

 [1]  TRUE  TRUE  TRUE  TRUE  TRUE FALSE FALSE FALSE FALSE FALSE

> x[consec<60] #첫번째, 두번째, 세번째 값

[1] 10  9  8  7  6

> seq(-pi,pi,1) #파이값을 1씩 증가

[1] -3.1415927 -2.1415927 -1.1415927 -0.1415927  0.8584073  1.8584073  2.8584073

> seq(-pi,pi,length=10) #파이값을 10개로 나눔

 [1] -3.1415927 -2.4434610 -1.7453293 -1.0471976 -0.3490659  0.3490659  1.0471976  1.7453293

 [9]  2.4434610  3.1415927

난수생성

> rnorm(10,-5,2.5) #평균이 -5, 표준편차가 2.5인 정규분포를 따르는 난수 10개

 [1] -5.66696754 -2.73338407 -3.51628322 -5.24746144 -7.56506891 -4.47006154 -6.19141094 -9.51129455

 [9] -0.08720922 -8.05256413

> rnorm(20) #평균이 0, 표준편차가 1인 정규분포를 따르는 난수 20개

 [1] -0.43086982 -0.40359489  1.93710961 -0.84071208 -0.10802247 -0.01138187  0.92168668 -2.05161187

 [9] -0.12887195 -0.27202527  1.73687933 -1.09407623  0.46089730  1.74135668  0.69786046 -0.38488478

[17]  1.05465527 -1.22405759  0.18905127  1.26340571

행렬연산

> x=c(1:12)

> x=matrix(x,ncol=4,byrow=T)

> x

     [,1] [,2] [,3] [,4]

[1,]    1    2    3    4

[2,]    5    6    7    8

[3,]    9   10   11   12

> y=matrix(x,ncol=4,byrow=F)

> y

     [,1] [,2] [,3] [,4]

[1,]    1    2    3    4

[2,]    5    6    7    8

[3,]    9   10   11   12

> y=c(1:12)

> y=matrix(y,ncol=4,byrow=F)

> y

     [,1] [,2] [,3] [,4]

[1,]    1    4    7   10

[2,]    2    5    8   11

[3,]    3    6    9   12

> x

     [,1] [,2] [,3] [,4]

[1,]    1    2    3    4

[2,]    5    6    7    8

[3,]    9   10   11   12

> x[,c(1:3)]

     [,1] [,2] [,3]

[1,]    1    2    3

[2,]    5    6    7

[3,]    9   10   11

> matrix(1,nrow=4,byrow=T)

     [,1]

[1,]    1

[2,]    1

[3,]    1

[4,]    1

> matrix(1,nrow=4,ncol=4)

     [,1] [,2] [,3] [,4]

[1,]    1    1    1    1

[2,]    1    1    1    1

[3,]    1    1    1    1

[4,]    1    1    1    1

> x[c(1:3),]

     [,1] [,2] [,3] [,4]

[1,]    1    2    3    4

[2,]    5    6    7    8

[3,]    9   10   11   12

> x[c(1:3),-2]

     [,1] [,2] [,3]

[1,]    1    3    4

[2,]    5    7    8

[3,]    9   11   12

> a=c(1:3,5:7)

> a

[1] 1 2 3 5 6 7

> x=matrix(a,ncol=3,byrow=T)

> x

     [,1] [,2] [,3]

[1,]    1    2    3

[2,]    5    6    7

> y=matrix(c(1:4),ncol=2,byrow=F)

> y

     [,1] [,2]

[1,]    1    3

[2,]    2    4

> t(x)%*%y

     [,1] [,2]

[1,]   11   23

[2,]   14   30

[3,]   17   37

apply함수

> apply(x,1,mean) #1은 행

[1] 2 6

> apply(x,2,mean) #2는 열

[1] 3 4 5

Function 문

> squareF=function(x){x*x}

> View(squreF)

> squareF(5)

[1] 25

기술통계량 구하기

> survey.data=read.table('survey.txt',header=T)

> survey.data

   seq sex marriage age job edu salary

1    1   1        1  21   1   4     60

2    2   1        1  22   5   5    100

3    3   1        1  33   1   4    200

4    4   2        2  33   7   4    200

5    5   1        2  28   1   4     70

6    6   1        1  21   5   5     80

7    7   2        2  39   7   4    190

8    8   1        1  32   1   4    100

9    9   1        2  44   3   1    120

10  10   1        2  55   4   4    110

11  11   2        2  46   7   5    150

12  12   1        1  20   1   4     50

13  13   1        2  31   6   4    210

14  14   1        1  27   1   4     60

15  15   2        1  21   5   5     80

> str(survey.data)

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

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

 $ sex     : int  1 1 1 2 1 1 2 1 1 1 ...

 $ marriage: int  1 1 1 2 2 1 2 1 2 2 ...

 $ age     : int  21 22 33 33 28 21 39 32 44 55 ...

 $ job     : int  1 5 1 7 1 5 7 1 3 4 ...

 $ edu     : int  4 5 4 4 4 5 4 4 1 4 ...

 $ salary  : int  60 100 200 200 70 80 190 100 120 110 ...

> survey.data$sex = as.factor(survey.data$sex)

> survey.data$marriage = as.factor(survey.data$marriage)

> survey.data$job = as.factor(survey.data$job)

> survey.data$edu = as.factor(survey.data$edu)

> str(survey.data)

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

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

 $ sex     : Factor w/ 2 levels "1","2": 1 1 1 2 1 1 2 1 1 1 ...

 $ marriage: Factor w/ 2 levels "1","2": 1 1 1 2 2 1 2 1 2 2 ...

 $ age     : int  21 22 33 33 28 21 39 32 44 55 ...

 $ job     : Factor w/ 6 levels "1","3","4","5",..: 1 4 1 6 1 4 6 1 2 3 ...

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

 $ salary  : int  60 100 200 200 70 80 190 100 120 110 ...

> attach(survey.data)

> mean(age)

[1] 31.53333

> sd(age)

[1] 10.57535

> summary(age)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 

  20.00   21.50   31.00   31.53   36.00   55.00 

> tapply(age,sex,mean)

       1        2 

30.36364 34.75000 

> tapply(age,sex,sd)

       1        2 

10.82841 10.59481 

> tapply(age,marriage,mean)

       1        2 

24.62500 39.42857 

> tapply(age,marriage,sd)

       1        2 

5.316752 9.571784 

> table(sex,marriage)

   marriage

sex 1 2

  1 7 4

  2 1 3

> sex.marriage=list(sex,marriage)

> table(sex.marriage)

              sex.marriage.2

sex.marriage.1 1 2

             1 7 4

             2 1 3

> sex.edu=list(survey.data$sex,survey.data$edu)

> sex.edu.tb=table(sex.edu)

> sex.edu.tb

         sex.edu.2

sex.edu.1  1  2  3  4  5

        1  1  1  1 13 11

        2  0  0  2  6  5

> rownames(sex.edu.tb)=c('남성','여성')

> colnames(sex.edu.tb)=c('무학','초졸','중졸','고졸','대졸')

> barplot(sex.edu.tb)

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

> pie(sex.edu.tb[1,])

> title('Education of Male')

> pie(sex.edu.tb[2,])

> title('Education of Female')

#산점도 사례

> math <-c(68,64,48,46,78,60,90,50,66,70)

> physics<-c(70,68,46,48,84,64,92,52,68,72)

> cor(math,physics) #두 변수의 상관계수 값

[1] 0.9921232

> plot(math,physics,main='수학과 물리 산점도')

##설문조사 사례##

문항 1. 귀하의 성별은?

1) 남자 2) 여자 문항

2. 결혼하셨습니까?

1) 미혼 2) 기혼 3) 이혼 문항

3. 귀하의 나이는?  (단위: 세)

4. 귀하의 직업은?

1) 회사원 2) 공무원 3) 노무자 4) 정치가 5) 학생 6) 기업가 7) 주부 8) 기타 문항

5. 귀하의 학력은?

1) 중졸이하 2)고졸 3) 대졸 4) 대학원 5) 그 이상

6. 가족의 월수입은? (단위: 만원)

> survey.data = read.table('survey.csv',sep=',',header=T)

> head(survey.data)

  seq sex marriage age job edu salary

1   1   1        1  21   1   4     60

2   2   1        1  22   5   5    100

3   3   1        1  33   1   4    200

4   4   2        2  33   7   4    120

5   5   1        2  28   1   4     70

6   6   1        1  21   5   5     80

> attach(survey.data)

> mean(age)

[1] 34.275

> sd(age)

[1] 11.60236

> summary(age)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.

  20.00   24.75   32.00   34.28   42.50   59.00

#그룹화

> tapply(survey.data$age,survey.data$sex,mean) #성별 나의 평균

       1        2

33.96296 34.92308

> tapply(survey.data$age,survey.data$sex,sd)

       1        2

11.96945 11.24323

> tapply(survey.data$age,survey.data$marriage,mean) #결혼 별 나이 평균

       1        2        3

24.66667 39.13043 50.50000

> tapply(survey.data$age,survey.data$marriage,sd)

        1         2         3

 4.151879 10.467718 12.020815

> sex.marriage=list(survey.data$sex,survey.data$marriage) #그룹

> table(sex.marriage) #분할표

              sex.marriage.2

sex.marriage.1  1  2  3

             1 10 15  2

             2  5  8  0

> tapply(survey.data$age,sex.marriage,mean)

     1        2    3

1 24.8 37.86667 50.5

2 24.4 41.50000   NA

> tapply(survey.data$age,sex.marriage,sd)

         1         2        3

1 4.709329 11.230486 12.02082

2 3.209361  9.071147       NA

> table(survey.data$sex) #빈도표

 1  2

27 13

> table(edu)

edu

 1  2  3  4  5

 1  1  3 19 16

> SexEdu=table(survey.data$sex,survey.data$edu)

> SexEdu

     1  2  3  4  5

  1  1  1  1 13 11

  2  0  0  2  6  5

어느 집단에서 표본을 10명 추출하여 다음과 같은 4개 문항에 대하여 설문조사를 실시하였다.

문항 1. 귀하의 성별은? 1) 남자 2) 여자 문항

2. 귀하의 나이는? (단위 : 세) 문항

3. 교육정도는? 1) 중졸이하 2) 고졸 3)대졸 및 그 이상 문항

4. 월수입(단위: 만원)

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

> colnames(ex8)<-c('sex','age','edu','salary')

> sex.tb=table(ex8$sex)

> sex.tb 

1 2

6 4

> edu.tb=table(ex8$edu)

> edu.tbe 

1 2 3

2 4 4

> rownames(edu.tb)=c('중졸이하','고졸이하','대졸이상')

> edu.tb

중졸이하 고졸이하 대졸이상

2        4        4

> barplot(edu.tb)

> plot(ex8$age,ex8$salary)

> title('Scatter plot of (Age,Salary)')

> plot(ex8$age,ex8$salary,type="n")

> points(ex8$age[sex=='1'],ex8$salary[sex=='1'],pch='M',col='blue')

> points(ex8$age[sex=='2'],ex8$salary[sex=='2'],pch='F',col='red')

> title("나이, 월수입 산점도:남녀별 구분")

 위의 그래프가 내가 그린 것이고, 아래는 교재. 왜 안되지?.. 좀더 공부해보자.

> edt.tb = table(survey.data$edu)

> edt.tb

 1  2  3  4  5

 1  1  3 19 16

> rownames(edt.tb) = c('무학','초졸','중졸','고졸','대졸')

> barplot(edu.tb)

> title('교육정도 막대그림')

> pie(edu.tb)

> title('교육정도 원그림')

Sex를 그룹변수로 하여 교육(edu)의 막대그림을 그려보도록 하자.

> sex_edu = list(survey.data$sex, survey.data$edu)

> sex_edu.tb = table(sex_edu)

> sex_edu.tb

         sex_edu.2

sex_edu.1  1  2  3  4  5

        1  1  1  1 13 11

        2  0  0  2  6  5

> colnames(sex_edu.tb) = c('무학','초졸','중졸','고졸','대졸')

> rownames(sex_edu.tb) = c('남자','여자')

> barplot(sex_edu.tb)

> title('성별 교육 정도 막대그림')

> boxplot(salary~sex,data=survey.data,xlab='성별',ylab='소득')

R에 내장된 데이터 co2를 이용, 선으로 그려보자(lines(smooth(co2)). col에서 4는 blue를 의미.

> plot(co2)

> lines(smooth(co2),col=4)

> x<-seq(0,20,0.1)

> y<-exp(-x/10)*cos(2*x)

> plot(x,y,type='line') #or type=’l’

> library(MVA)

> library(HSAUR2)

> data(USairpollution)

> head(USairpollution)

            SO2 temp manu popul wind precip predays

Albany       46 47.6   44   116  8.8  33.36     135

Albuquerque  11 56.8   46   244  8.9   7.77      58

Atlanta      24 61.5  368   497  9.1  48.34     115

Baltimore    47 55.0  625   905  9.6  41.31     111

Buffalo      11 47.1  391   463 12.4  36.11     166

Charleston   31 55.2   35    71  6.5  40.75     148

> plot(USairpollution$manu,USairpollution$popul)

> x=USairpollution[,c(3,4)]

> bvbox(x,xlab='manu',ylab='popul')

> title('bivariate boxplot') # bivariate 이변량 

> identify(x)

> rownames(x)[c(7,9,14,30)]

[1] "Chicago"      "Cleveland"    "Detroit"      "Philadelphia"

> plot(wind~temp,data=USairpollution,pch=9)

> with(USairpollution,symbols(temp,wind,circles=SO2,inches=0.5,add=T))

> title('Bubble Plot')

등고선 그림(contour)과 이미지 그림(image plot)

기상자료, 지리 위치 등을 나타내는 자료는 좌표상에 z변수를 등고선이나, 이미지 등으로 표시하게 되면 매우 효율적임

- 등고선 자료:

자료를 나타내기 위하 두 좌표의 격자(grad)위에 등고선으로 면을 그려줌

- 이미지 그림:

색깔을 이용하여 z변수의 이미지를 표현한 그래프

> x<--6:16

> x

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

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

> contour(outer(x,x),method='edge',vfont=c('sans serif','plain'))

> z<-outer(x,sqrt(abs(x)),FUN='/')

> image(x,x,z)

> contour(x,x,z,col='pink',add=TRUE,method='edge',vfont=c('sans serif','plain'))

> contour(x,x,z,ylim=c(1,6),method='simple',laboex=1)

Warning messages:

1: In plot.window(xlim, ylim, ...) : "laboex" is not a graphical parameter

2: In title(...) : "laboex" is not a graphical parameter

3: In axis(side = side, at = at, labels = labels, ...) :

  "laboex" is not a graphical parameter

4: In axis(side = side, at = at, labels = labels, ...) :

  "laboex" is not a graphical parameter

5: In box(...) : "laboex" is not a graphical parameter

> contour(x,x,z,ylim=c(-6,6),nlev=20,lty=2,method='simple')

겨냥도 그림(Perspective Plot)

등간격의 격자상에 높이 값을 갖는 행렬 자료에 대한 3차원 표현 방법 높이는 선으로 연결

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

> x<-seq(-10,10,length=30)

> y<-x

> f<-function(x,y){r<-sqrt(x^2+y^2);10*sin(r)/r}

> z<-outer(x,y,f)

> z[is.na(z)]<-l

Error: object 'l' not found

> z[is.na(z)]<-1

> persp(x,y,z,theta=30,phi=30,expand=0.5,col='lightblue')

> persp(x,y,z,theta=30,phi=30,expand=0.5,col='lightblue',ltheta=120,shade=0.75,thicktype='detailed',xlab='X',ylab='Y',zlab='Sinc(r)')

Warning message:

In persp.default(x, y, z, theta = 30, phi = 30, expand = 0.5, col = "lightblue",  :

  "thicktype" is not a graphical parameter

> edu<-read.csv('education.csv')

> edu

   grade y1970 y1990

1 초등졸  23.8  21.1

2   중졸  24.7  20.2

3   고졸  36.0  43.6

4   대졸  11.0  15.1

> name=edu[,1]

> grade70=edu[,2]

> grade90=edu[,3]

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

> pie(grade70,label=name,clockwise=T,col=c(2:5),main='1970년도 여성학력')

> mtext(outer=F,'(a)',side=1,line=-6)

> pie(grade90,label=name,clockwise=T,col=c(2:5),main='1990년도 여성학력')

> mtext(outer=F,'(b)',side=1,line=-6)

연습문제

1. 13개 시중은행에 대한 편리성, 신속성, 친절, 능률, 쾌적, 자동화등의 점수. 출처: 유종렬 외 3인, S-Plus를 이용한 통계계산, 박영사, 1997

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

> head(bank)

         X convenience accuracy kindness efficiency pleasant automatic

1  kookmin        70.5     59.4     63.7       54.3     66.9      62.6

2  enterpr        64.8     70.3     68.6       55.2     68.0      64.1

3    boram        67.1     79.6     78.5       62.4     79.8      62.4

4 commerce        61.1     65.0     65.6       54.4     64.5      63.9

5    seoul        63.4     66.5     67.9       65.0     59.7      62.0

6  shinhan        72.3     69.1     74.2       60.0     70.1      68.2

> bank<-bank[,-1]

> name<-c('kookmin','enterpr','boram','commerce','seoul','shinhan','city','exchange','first','chohung','hana','hanil','house')

> head(bank)

  convenience accuracy kindness efficiency pleasant automatic

1        70.5     59.4     63.7       54.3     66.9      62.6

2        64.8     70.3     68.6       55.2     68.0      64.1

3        67.1     79.6     78.5       62.4     79.8      62.4

4        61.1     65.0     65.6       54.4     64.5      63.9

5        63.4     66.5     67.9       65.0     59.7      62.0

6        72.3     69.1     74.2       60.0     70.1      68.2

(1) 각 변수들의 히스토그램을 그리고 설명

> library(dplyr)

> nrow(bank)

[1] 13

> bank$convenience %>% summary

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.

  61.10   63.50   64.80   65.98   68.40   72.30

> bank$convenience %>% mean

[1] 65.97692

> bank$convenience %>% sd

[1] 3.279521

>  bank$convenience %>% hist

> con.h<-bank %>% filter(convenience > mean(bank$convenience))

> con.l<-bank %>% filter(convenience < mean(bank$convenience))

> con.h

  convenience accuracy kindness efficiency pleasant automatic

1        70.5     59.4     63.7       54.3     66.9      62.6

2        67.1     79.6     78.5       62.4     79.8      62.4

3        72.3     69.1     74.2       60.0     70.1      68.2

4        68.4     67.5     67.3       51.3     71.3      65.8

5        66.1     66.5     67.3       50.7     63.4      63.3

6        69.0     74.3     80.5       63.6     75.7      55.9

> con.h$convenience %>% summary

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.

  66.10   67.42   68.70   68.90   70.12   72.30

> con.l$convenience %>% summary

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.

  61.10   63.30   63.50   63.47   64.15   64.80

> con.h %>% summarise('평균 친절도' = mean(kindness))

  평균 친절도

1    71.91667

> con.l %>% summarise('평균 친절도'=mean(kindness))

  평균 친절도

1        67.7

2) 산점도 행렬

> plot(bank)

3) 별그림과 얼굴그림

> rownames(bank) = name

> head(bank)

         convenience accuracy kindness efficiency pleasant automatic

kookmin         70.5     59.4     63.7       54.3     66.9      62.6

enterpr         64.8     70.3     68.6       55.2     68.0      64.1

boram           67.1     79.6     78.5       62.4     79.8      62.4

commerce        61.1     65.0     65.6       54.4     64.5      63.9

seoul           63.4     66.5     67.9       65.0     59.7      62.0

shinhan         72.3     69.1     74.2       60.0     70.1      68.2

> stars(bank)

> library(aplpack)

Loading required package: tcltk

> faces(bank,face.type=1)

effect of variables:

 modified item       Var         

 "height of face   " "convenience"

 "width of face    " "accuracy"  

 "structure of face" "kindness"  

 "height of mouth  " "efficiency"

 "width of mouth   " "pleasant"  

 "smiling          " "automatic" 

 "height of eyes   " "convenience"

 "width of eyes    " "accuracy"  

 "height of hair   " "kindness"  

 "width of hair   "  "efficiency"

 "style of hair   "  "pleasant"  

 "height of nose  "  "automatic" 

 "width of nose   "  "convenience"

 "width of ear    "  "accuracy"  

 "height of ear   "  "kindness"

출처: R을 이용한 데이터 마이닝(김성수, 김현중, 정성석, 이용구 공저)

연속형 목표 변수.

목표변수가 연속형인 경우 모형의 예측력 측도로는 주로 예측평균제곱오차(PMSE: Prediction Mean Squared Error)를 사용한다.

PMSE = sigma(yi-yhat i)²/n

관측치(yi)와 예측치(yhat i) 사이에 차이가 적을수록 예측을 잘한 것으로므로 예측평균제곱오차가 작을수록 그 모형의 예측력이 높다고 할 수 있다. 또한 관측치와 예측치를 가로축 및 세로축으로 놓고 그린 산점도(scatter plot)가 45도 대각선을 중심으로 모여 있으면 예측력이 좋다고 할 수 있다.

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

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

<데이터 준비>

> library(MASS)

> str(Boston)

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

> dim(Boston)

[1] 506  15

> Boston<-Boston[,-15]

> str(Boston)

'data.frame':   506 obs. of  14 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 ...

1. 선형회귀모형

> i=sample(1:nrow(Boston),round(nrow(Boston)*0.7))

> Boston.train = Boston[i,] #70% for training data

> Boston.test = Boston[-i,] #30% for test data

> fit.reg = lm(medv~.,data = Boston.train) #선형회귀모형

> summary(fit.reg) #선형회귀모형 요약

Call:

lm(formula = medv ~ ., data = Boston.train)

Residuals:

    Min      1Q  Median      3Q     Max

-9.9375 -2.7260 -0.4588  1.7937 24.3220

Coefficients:

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

(Intercept)  41.027525   6.644569   6.175 1.93e-09 ***

crim         -0.098730   0.052682  -1.874 0.061794 . 

zn            0.059161   0.016208   3.650 0.000304 ***

indus         0.073377   0.075763   0.969 0.333492   

chas1         2.314962   1.023899   2.261 0.024409 * 

nox         -18.307923   4.638939  -3.947 9.67e-05 ***

rm            3.014437   0.498806   6.043 4.04e-09 ***

age          -0.002066   0.016414  -0.126 0.899922   

dis          -1.564686   0.235506  -6.644 1.25e-10 ***

rad2          1.485691   1.707949   0.870 0.384999   

rad3          5.371463   1.513283   3.550 0.000441 ***

rad4          2.763349   1.365438   2.024 0.043791 * 

rad5          3.407697   1.380288   2.469 0.014057 * 

rad6          1.840557   1.687536   1.091 0.276204   

rad7          3.790169   1.832307   2.069 0.039362 * 

rad8          4.691472   1.677890   2.796 0.005473 **

rad24         8.224882   2.114575   3.890 0.000121 ***

tax          -0.009025   0.004555  -1.981 0.048382 * 

ptratio      -1.021842   0.175757  -5.814 1.43e-08 ***

black         0.008734   0.003633   2.404 0.016749 * 

lstat        -0.596699   0.059397 -10.046  < 2e-16 ***

---

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

Residual standard error: 4.658 on 333 degrees of freedom

Multiple R-squared:  0.7468,    Adjusted R-squared:  0.7316

F-statistic: 49.11 on 20 and 333 DF,  p-value: < 2.2e-16

> fit.step.reg = step(fit.reg,direction = 'both', trace = FALSE) #단계적 선택법, trace=T일 경우 AIC를 볼 수 있다.

> summary(fit.step.reg) #조정된 모델 요약

Call:

lm(formula = medv ~ crim + zn + chas + nox + rm + dis + rad +

    tax + ptratio + black + lstat, data = Boston.train)

Residuals:

    Min      1Q  Median      3Q     Max

-9.8479 -2.6319 -0.4838  1.7225 24.3723

Coefficients:

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

(Intercept)  40.432105   6.597006   6.129 2.49e-09 ***

crim         -0.102050   0.052455  -1.945 0.052555 . 

zn            0.058012   0.015903   3.648 0.000307 ***

chas1         2.428119   1.014994   2.392 0.017296 * 

nox         -17.025847   4.237976  -4.017 7.27e-05 ***

rm            2.969435   0.489146   6.071 3.45e-09 ***

dis          -1.604642   0.221607  -7.241 3.08e-12 ***

rad2          1.763045   1.681098   1.049 0.295051   

rad3          5.319042   1.507369   3.529 0.000476 ***

rad4          2.871626   1.358651   2.114 0.035290 * 

rad5          3.377174   1.377514   2.452 0.014731 * 

rad6          1.669858   1.672013   0.999 0.318656   

rad7          3.831541   1.823028   2.102 0.036322 * 

rad8          4.551858   1.667065   2.730 0.006659 **

rad24         7.906985   2.076704   3.807 0.000167 ***

tax          -0.007384   0.004213  -1.753 0.080567 . 

ptratio      -0.996843   0.173266  -5.753 1.98e-08 ***

black         0.008481   0.003612   2.348 0.019458 *  

lstat        -0.593839   0.055467 -10.706  < 2e-16 ***

---

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

Residual standard error: 4.65 on 335 degrees of freedom

Multiple R-squared:  0.7461,    Adjusted R-squared:  0.7324

F-statistic: 54.69 on 18 and 335 DF,  p-value: < 2.2e-16

> fit.step.reg$anova #제거된 변수 보기. age와 indus가 제거

     Step Df   Deviance Resid. Df Resid. Dev      AIC

1         NA         NA       333   7223.934 1109.614

2   - age  1  0.3436196       334   7224.278 1107.631

3 - indus  1 20.2559598       335   7244.534 1106.622

> yhat.reg = predict(fit.step.reg,newdata=Boston.test,type='response') #데이터 예측, 목표값 예측 시 type=’response’

> mean((Boston.test$medv-yhat.reg)^2) #예측 모형의 평균 제곱 오차 PMSE

[1] 23.68968

2. 회귀 나무모형

> library(rpart)

> set.seed(1234)

> i=sample(1:nrow(Boston),round(nrow(Boston)*0.7))

> Boston.train = Boston[i,]

> Boston.test = Boston[-i,]

> my.control = rpart.control(xval=10, cp=0, minsplit=nrow(Boston.train)*0.05)

> fit.tree = rpart(medv~.,data=Boston.train, method='anova',control=my.control)#method=anova 선택 시 회귀나무, class 선택 시 분류나무

> which.min(fit.tree$cp[,4])

17

17

> fit.tree$cp[17]

[1] 0.00193237

> fit.tree$cp[17,]

         CP      nsplit   rel error      xerror        xstd

 0.00193237 17.00000000  0.12723520  0.23200911  0.04362543

> fit.tree$cp[17,1]

[1] 0.00193237

> ii = which.min(fit.tree$cp[,4])

> fit.prune.tree = prune(fit.tree,cp=fit.tree$cp[ii,1])

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

> text(fit.prune.tree,col='blue',cex=0.7)

> yhat.tree = predict(fit.prune.tree,newdata=Boston.test,type='vector')

> mean((Boston.test$medv - yhat.tree)^2)

[1] 24.17478

 3. 신경망모형

> library(neuralnet)

> Boston2 = Boston

> for(i in 1:ncol(Boston2))if(!is.numeric(Boston2[,i]))Boston2[,i]=as.numeric(Boston[,i])

> max1 = apply(Boston2,2,max)

> min1 = apply(Boston2,2,min)

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

> sdat=as.data.frame(sdat)

> View(sdat)

> set.seed(1234)

> i=sample(1:nrow(Boston2),round(nrow(Boston2)*0.7))

> Boston2.train = sdat[i,]

> Boston2.test = sdat[-i,]

> vname = names(Boston2.train)

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

> form

medv ~ crim + zn + indus + chas + nox + rm + age + dis + rad +

    tax + ptratio + black + lstat

> fit.nn = neuralnet(form,data=Boston2.train,hidden=c(3,2),linear.output=T)#linear.output=T 회귀

> plot(fit.nn) 

> str(Boston.test)

'data.frame':   152 obs. of  14 variables:

 $ crim   : num  0.0273 0.0299 0.1446 0.6274 1.0539 ...

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

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

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

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

 $ rm     : num  7.18 6.43 6.17 5.83 5.93 ...

 $ age    : num  61.1 58.7 96.1 56.5 29.3 98.1 89.2 91.7 85.7 90.3 ...

 $ dis    : num  4.97 6.06 5.95 4.5 4.5 ...

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

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

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

 $ black  : num  393 394 397 396 387 ...

 $ lstat  : num  4.03 5.21 19.15 8.47 6.58 ...

 $ medv   : num  34.7 28.7 27.1 19.9 23.1 13.6 19.6 15.2 13.9 16.6 ...

> pred=compute(fit.nn,Boston2.test[,1:13])#14번째 변수가 medv 목표변수.compute는 작성된 신경망모형을 이용하여 새로운 예에 적용하여 결과를 도출.

> summary(pred)

           Length Class  Mode  

neurons      3    -none- list  

net.result 152    -none- numeric

> yhat.nn = pred$net.result*(max(Boston2$medv)-min(Boston2$medv))+min(Boston2$medv)#back transformation, 목표변수를 조정전 변수값으로 환원하는 과정. 원래 가지고 있는 값으로 환원. 역변환 과정

> head(cbind(Boston2[i,14],yhat.nn),50)

    [,1]        [,2]

2   31.6 21.17115946

3   23.8 30.30760801

10  28.2 20.06686203

13  21.6 21.29808480

14  16.1 20.28145030

18  23.8 18.30584489

24  36.2 16.10587573

26  21.2 16.16710531

27  17.1 17.71479369

28  20.9 16.71171622

33  31.2 12.87970542

35  20.7 14.84095463

36  17.8 22.05195863

42  14.1 33.16705465

43  15.6 26.86733582

46  17.2 21.39641163

47  14.0 20.62187769

49  19.2 13.55437741

51  22.0 21.63602446

57  18.7 23.04343250

59  19.4 23.33368068

72  15.6 21.68452217

75  20.0 22.95006783

76  18.2 21.37727858

81  19.5 27.88596126

90  11.5 31.97882170

94  29.6 24.56223312

98   8.4 41.94592034

102  8.5 23.20282812

103 19.6 27.89450422

104 28.7 19.47439349

105 21.4 19.90631048

111 11.8 20.49063003

120 23.3 20.15009581

123 26.6 19.19281571

128 21.7 16.99190292

135 20.6 14.54799220

148 20.3 13.09572781

150 21.4 15.35938684

151 13.3 18.92860416

161 50.0 31.83885454

165 24.8 20.14597079

167 29.8 44.60929619

169 23.2 20.70938882

170 15.3 20.60174521

171 31.7 19.12045877

172 22.1 19.63836241

173 27.5 20.36004622

176 22.8 32.60724125

177 26.6 24.24612266

Warning message:

In cbind(Boston2[i, 14], yhat.nn) :

  number of rows of result is not a multiple of vector length (arg 1)

> Boston2.test$medv=Boston2.test$medv*(max(Boston2$medv)-min(Boston2$medv))+min(Boston2$medv)#back transformation

> mean((Boston2.test$medv-yhat.nn)^2)

[1] 19.34365205

4. 랜덤 포레스트

> library(randomForest)

> library(MASS)

> str(Boston)

'data.frame':   506 obs. of  14 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 ...

> set.seed(1234)

> i=sample(1:nrow(Boston),round(nrow(Boston)*0.7))

> Boston.train = Boston[i,]

> Boston.test = Boston[-i,]

> fit.rf = randomForest(medv~., data=Boston.train, ntree=100, mtry=5, importance = T, na.action=na.omit)

> fit.rf

Call:

 randomForest(formula = medv ~ ., data = Boston.train, ntree = 100,      mtry = 5, importance = T, na.action = na.omit)

               Type of random forest: regression

                     Number of trees: 100

No. of variables tried at each split: 5

          Mean of squared residuals: 11.60448807

                    % Var explained: 85.6

> summary(fit.rf)

                Length Class  Mode    

call              7    -none- call    

type              1    -none- character

predicted       354    -none- numeric 

mse             100    -none- numeric 

rsq             100    -none- numeric 

oob.times       354    -none- numeric 

importance       26    -none- numeric 

importanceSD     13    -none- numeric 

localImportance   0    -none- NULL    

proximity         0    -none- NULL    

ntree             1    -none- numeric 

mtry              1    -none- numeric 

forest           11    -none- list    

coefs             0    -none- NULL    

y               354    -none- numeric 

test              0    -none- NULL    

inbag             0    -none- NULL    

terms             3    terms  call    

> head(fit.rf$predicted,30)

         58         315         308         314         433         321

31.38358586 22.70753968 28.77652778 21.11594595 20.61043860 23.92682540

          5         117         332         256         345         270

32.56341667 20.64287500 20.64851351 20.71283333 30.29218391 21.12653153

        140         456         144         412         141         131

15.02797980 15.33454167 14.36800926 14.37462963 15.01537037 19.89743590

         92         114         154         147          77          20

23.40651042 19.04477778 15.95588235 16.05761261 21.11859649 19.79196970

        106         390         253         439         398          22

19.35319820 10.55704563 31.93869369  9.22796875 12.18531532 18.03844203

> importance(fit.rf,type=1)

              %IncMSE

crim     6.3818661727

zn       1.5083317623

indus    4.9357684617

chas     0.6538995939

nox      6.8146688870

rm      15.9443905699

age      3.6616801321

dis      7.0604314940

rad      6.1703029034

tax      3.9762473150

ptratio  7.1788230837

black    3.4813158641

lstat   13.5182724091

> yhat.rf=predict(fit.rf,newdata=Boston.test,type='response')

> mean((Boston.test$medv-yhat.rf)^2)

[1] 12.41384092

각 모형 별 에측평균제곱오차

관측치와 예측치의 산점도 생성

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

>plot(Boston.test$medv,yhat.reg,xlab='Observed',ylab='Predicted',main='Regression',xlim=c(0,55),ylim=c(0,55))

> abline(0,1,lty=2)

> plot(Boston.test$medv,yhat.tree,xlab='Observed',ylab='Predicted',main='Decision tree',xlim=c(0,55),ylim=c(0,55))

> abline(0,1,lty=2)

> plot(Boston.test$medv,yhat.nn,xlab='Observed',ylab='Predicted',main='Neural Network',xlim=c(0,55),ylim=c(0,55))

> abline(0,1,lty=2)

> plot(Boston.test$medv,yhat.rf,xlab='Observed',ylab='Predicted',main='Random Forest',xlim=c(0,55),ylim=c(0,55))

> abline(0,1,lty=2) 

 -> 랜덤포레스트의 그래프가 가장 조밀하게 형성되어, 가장 좋은 예측 결과를 가진다고 할 수 있다.

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

문제인식: 매스미디어 광고에 의한 신규유저수가 일정치 않다. 이는 매월 TV광고와 잡지 광고의배분이 일정하지 않기 때문이다. 이에 TV, 잡지 광고비와 신규 유저수의 관계를 파악한다.

해결 방법: TV, 잡지 광고비와 신규 유저수 데이터를 기반으로 중회귀분석을 실시한다.

R

1. 데이터 읽어 들이기

> ad.data <- read.csv('ad_result.csv',header=T,stringsAsFactors = F)

> ad.data

month  tvcm magazine install

1  2013-01  6358     5955   53948

2  2013-02  8176     6069   57300

3  2013-03  6853     5862   52057

4  2013-04  5271     5247   44044

5  2013-05  6473     6365   54063

6  2013-06  7682     6555   58097

7  2013-07  5666     5546   47407

8  2013-08  6659     6066   53333

9  2013-09  6066     5646   49918

10 2013-10 10090     6545   59963

2. TV 광고의 광고비용과 신규 유저수의 산점도 그리기

> library(ggplot2)
> library(scales)

> ggplot(ad.data,aes(x=tvcm,y=install))

> ggplot(ad.data,aes(x=tvcm,y=install))+geom_point()

> ggplot(ad.data,aes(x=tvcm,y=install))+geom_point()+xlab('TV 광고비')+ylab('신규유저수')

> ggplot(ad.data,aes(x=tvcm,y=install))+geom_point()+xlab('TV 광고비')+ylab('신규유저수')+scale_x_continuous(label=comma)+scale_y_continuous(label=comma)

3. 잡지 광고의 광고비용과 신규 유저수의 산점도 그리기

> ggplot(ad.data,aes(x=magazine,y=install))+geom_point()+xlab('잡지 광고비')+ylab('신규유저수')+scale_x_continuous(label=comma)+scale_y_continuous(label=comma)

4. 회귀분석 실행

> fit <-lm(install~.,data=ad.data[,c('install','tvcm','magazine')])

> fit

Call:

  lm(formula = install ~ ., data = ad.data[, c("install", "tvcm",

                                               "magazine")])

Coefficients:

  (Intercept)         tvcm     magazine 

188.174        1.361        7.250 

-> 이상 내용으로부터 아래와 같은 모델을 만들 수 있다.

신규 유저수 = 1.361 X TV광고비 + 7.25 X 잡지광고비 + 188.174

이라는 관계가 있으며, 신규 유저는 광고를 실시하지 않을 때 월 188명 정도이다. (아래 summary(fit)을 통해 유의하지 않음을 알 수 있다) 그리고 TV 광고에 1만원을 투입하면 양 1.3609명의 신규 고객을, 잡지 광고에 1만원을 투자하면 약 7.2498명의 신규 유저를 확보할 수 있다.

5. 회귀분석의 결과 해석

> summary(fit)

Call:

  lm(formula = install ~ ., data = ad.data[, c("install", "tvcm", "magazine")])

Residuals:

  Min       1Q   Median       3Q      Max

-1406.87  -984.49   -12.11   432.82  1985.84

Coefficients:

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

(Intercept)  188.1743  7719.1308   0.024  0.98123  

tvcm         1.3609     0.5174   2.630  0.03390 *

magazine     7.2498     1.6926   4.283  0.00364 **

---

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

Residual standard error: 1387 on 7 degrees of freedom

Multiple R-squared:  0.9379,           Adjusted R-squared:  0.9202

F-statistic: 52.86 on 2 and 7 DF,  p-value: 5.967e-05

잔차(Residuals) : 잔차(예측값과 측정값의 차이)분포를 사분위수로 표현한 것으로, 데이터의 치우침이 있는지 확인할 수 있다. 1Q의 절대값이 3Q의 절대값보다 커서 약간 치우침이 있어 보인다.

Coefficients : 절편과 기울기에 관한 개요.

Adjusted R-Squared : 0.9202로 이 모델로 전체 데이터의 약 92.02%를 설명할 수 있다.

또는 아래와 같이 회귀분석 모델을 만들 수 있다.

> fit2<-lm(install~tvcm+magazine,data=ad.data)

> summary(fit2)

Call:

  lm(formula = install ~ tvcm + magazine, data = ad.data)

Residuals:

  Min       1Q   Median       3Q      Max

-1406.87  -984.49   -12.11   432.82  1985.84

Coefficients:

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

(Intercept)  188.1743  7719.1308   0.024  0.98123  

tvcm           1.3609     0.5174   2.630  0.03390 *

magazine       7.2498     1.6926   4.283  0.00364 **

---

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

Residual standard error: 1387 on 7 degrees of freedom

Multiple R-squared:  0.9379,              Adjusted R-squared:  0.9202

F-statistic: 52.86 on 2 and 7 DF,  p-value: 5.967e-05

출처: 비지니스 활용 사레로 배우는 데이터 분석: R (사카마키 류지, 사토 요헤이 지음)

문제인식: 자사 제품의 구매율이 다른 앱과 비교하여 낮음을 확인, 해당 배너 광고의 클릭율이 낮음을 밝힘.

해결 방법: 분석 위한 데이터가 없기 때문에  A/B 테스트를 통해 분석용 로그를 출력하는 데이터 수집. 이후 A/B 테스트 한 후 배너광고 A와 B의 클릭률 산출하여 x² 검정 실행

R 분석 내용

1. 데이터 읽어 들이기

> setwd("C:/파일/temp/r_temp")

> ab.test.imp <-read.csv('section5-ab_test_imp.csv',header=T,stringsAsFactors = F)

> head(ab.test.imp,2)

log_date app_name  test_name test_case user_id transaction_id

1 2013-10-01  game-01 sales_test         B   36703          25622

2 2013-10-01  game-01 sales_test         A   44339          25623

> ab.test.goal <- read.csv('section5-ab_test_goal.csv',header=T,stringsAsFactors = F)

> head(ab.test.goal,2)

log_date app_name  test_name test_case user_id transaction_id

1 2013-10-01  game-01 sales_test         B   15021          25638

2 2013-10-01  game-01 sales_test         B     351          25704

2. ab.test.imp에 ab.test.goal 결합시키기

> ab.test.imp <- merge(ab.test.imp, ab.test.goal, by='transaction_id', all.x=T, suffixes=c('','.g'))

> head(ab.test.imp,2)

-> all.x - logical; if TRUE, then extra rows will be added to the output, one for each row in x that has no matching row in y. These rows will have NAs in those columns that are usually filled with values from y. The default is FALSE, so that only rows with data from both x and y are included in the output.

3. 클릭했는지 하지 않았는지 나타내는 플래그 작성하기

> ab.test.imp$is.goal <-ifelse(is.na(ab.test.imp$user_id.g),0,1)

-> 항목 user_id.g의 값이 없음(NA)인 경우 0, 그렇지 않을 경우 1 입력

4. 클릭률 집계하기

> library(plyr)

> ddply(ab.test.imp,

           .(test_case),

           summarize,

           cvr=sum(is.goal)/length(user_id))

test_case        cvr

1         A 0.08025559

2         B 0.11546015

-> cvr = 클릭한 사람의 합계 / 배너광고가 표시된 유저수

    배너광고 A의 클릭률 보다 B가 전반적으로 높다는 것을 통계적으로 알 수 있다.

5. 카이제곱 검정

> chisq.test(ab.test.imp$test_case,ab.test.imp$is.goal)

Pearson's Chi-squared test with Yates' continuity correction

data:  ab.test.imp$test_case and ab.test.imp$is.goal

X-squared = 308.38, df = 1, p-value < 2.2e-16

-> p-value가 매주 작으므로 (0.05보다 작음) A와 B의 차이는 우연히 발생한 것이 아님

6. 날짜별,테스트 케이스별로 클릭률 산출하기

> ab.test.imp.summary<- ddply(ab.test.imp,

                                 .(log_date,test_case),

                                 summarize,

                                 imp=length(user_id),

                                 cv=sum(is.goal),

                                 cvr=sum(is.goal)/length(user_id))

> head(ab.test.imp.summary)

log_date test_case  imp  cv        cvr

1 2013-10-01         A 1358  98 0.07216495

2 2013-10-01         B 1391 176 0.12652768

3 2013-10-02         A 1370  88 0.06423358

4 2013-10-02         B 1333 212 0.15903976

5 2013-10-03         A 1213 170 0.14014839

6 2013-10-03         B 1233 185 0.15004055

-> imp = 배너 광고가 얼마나 표시되었는지 확인

cv = 클릭한 사람의 합계

cvr = 클릭률

7. 테스트 케이스별로 클릭률 산출하기

> ab.test.imp.summary<- ddply(ab.test.imp.summary,

                                 .(test_case),

                                 transform,

                                 cvr.avg=sum(cv/sum(imp)))

> head(ab.test.imp.summary)

log_date test_case  imp  cv        cvr    cvr.avg

1 2013-10-01         A 1358  98 0.07216495 0.08025559

2 2013-10-02         A 1370  88 0.06423358 0.08025559

3 2013-10-03         A 1213 170 0.14014839 0.08025559

4 2013-10-04         A 1521  89 0.05851414 0.08025559

5 2013-10-05         A 1587  56 0.03528670 0.08025559

6 2013-10-06         A 1219 120 0.09844135 0.08025559

-> transform: 기존 집계툴에 데이터 추가

   cvr.avg: test_case 별 추가

8. 테스트 케이스별 클릭률의 시계열 추이 그래프

> str(ab.test.imp.summary)

'data.frame':          62 obs. of  6 variables:

 $ log_date : chr  "2013-10-01" "2013-10-02" "2013-10-03" "2013-10-04" ...

$ test_case: chr  "A" "A" "A" "A" ...

$ imp      : int  1358 1370 1213 1521 1587 1219 1595 1401 1648 1364 ...

$ cv       : num  98 88 170 89 56 120 194 99 199 114 ...

$ cvr      : num  0.0722 0.0642 0.1401 0.0585 0.0353 ...

$ cvr.avg  : num  0.0803 0.0803 0.0803 0.0803 0.0803 ...

> ab.test.imp.summary$log_date <- as.Date(ab.test.imp.summary$log_date)

> library(ggplot2)

> library(scales)

> limits<-c(0,max(ab.test.imp.summary$cvr))

> ggplot(ab.test.imp.summary,aes(x=log_date,y=cvr,

                                    col=test_case,lty=test_case,shape=test_case))+

                                    geom_line(lwd=1)+geom_point(size=4)+geom_line(aes(y=cvr.avg,col=test_case))+

                                    scale_y_continuous(label=percent,limits=limits)

-> 배너광고 A의 클릭률 보다 B가 전반적으로 높다는 것을 통계적으로 알 수 있다

##ggplot 그래프 이해하기

> library(ggplot2)

> library(scales)

> ggplot(ab.test.imp.summary,aes(x=log_date,y=cvr,col=test_case,lty=test_case,shape=test_case))

> ggplot(ab.test.imp.summary,aes(x=log_date,y=cvr,col=test_case,lty=test_case,shape=test_case))+geom_line(lwd=1)

> ggplot(ab.test.imp.summary,aes(x=log_date,y=cvr,col=test_case,lty=test_case,shape=test_case))+geom_line(lwd=1)+geom_point(size=4)

> ggplot(ab.test.imp.summary,aes(x=log_date,y=cvr,col=test_case,lty=test_case,shape=test_case))+geom_line(lwd=1)+geom_point(size=4)+geom_line(aes(y=cvr.avg,col=test_case))

> limits<-c(0,max(ab.test.imp.summary$cvr))

> max(ab.test.imp.summary$cvr)

[1] 0.1712159

> ggplot(ab.test.imp.summary,aes(x=log_date,y=cvr,col=test_case,lty=test_case,shape=test_case))+geom_line(lwd=1)+geom_point(size=4)+geom_line(aes(y=cvr.avg,col=test_case))+ scale_y_continuous(label=percent,limits=limits)

출처: 비지니스 활용 사레로 배우는 데이터 분석: R (사카마키 류지, 사토 요헤이 지음)

1. csv 파일 읽어 들이기

> setwd("C:/파일/temp/r_temp")

> dau<-read.csv('section4-dau.csv',header=T,stringsAsFactors = F)

> head(dau)

     log_date app_name user_id

1 2013-08-01  game-01   33754

2 2013-08-01  game-01   28598

3 2013-08-01  game-01   30306

4 2013-08-01  game-01     117

5 2013-08-01  game-01    6605

6 2013-08-01  game-01     346

> user.info<-read.csv('section4-user_info.csv',header=T,stringsAsFactors = F)

> head(user.info)

i     nstall_date app_name user_id gender generation device_type

1   2013-04-15  game-01       1      M         40         iOS

2   2013-04-15  game-01       2      M         10     Android

3   2013-04-15  game-01       3      F         40         iOS

4   2013-04-15  game-01       4      M         10     Android

5   2013-04-15  game-01       5      M         40         iOS

6   2013-04-15  game-01       6      M         40         iOS

2. DAU에 user.info 결합하기

> dau.user.info<-merge(dau,user.info,by=c('user_id','app_name'))

> head(dau.user.info)

     user_id app_name   log_date   install_date  gender  generation  device_type

1       1  game-01 2013-09-06   2013-04-15      M         40         iOS

2       1  game-01 2013-09-05   2013-04-15      M         40         iOS

3       1  game-01 2013-09-28   2013-04-15      M         40         iOS

4       1  game-01 2013-09-12   2013-04-15      M         40         iOS

5       1  game-01 2013-09-11   2013-04-15      M         40         iOS

6       1  game-01 2013-09-08   2013-04-15      M         40         iOS

3. 월 항목 추가

> dau.user.info$month<-substr(dau.user.info$log_date,1,7)

> head(dau.user.info)

     user_id app_name   log_date   install_date  gender generation  device_type  month

1       1  game-01 2013-09-06   2013-04-15      M         40         iOS 2013-09

2       1  game-01 2013-09-05   2013-04-15      M         40         iOS 2013-09

3       1  game-01 2013-09-28   2013-04-15      M         40         iOS 2013-09

4       1  game-01 2013-09-12   2013-04-15      M         40         iOS 2013-09

5       1  game-01 2013-09-11   2013-04-15      M         40         iOS 2013-09

6       1  game-01 2013-09-08   2013-04-15      M         40         iOS 2013-09

> colnames(dau.user.info)[8]<-'log_month'

> head(dau.user.info)

user_id app_name   log_date  install_date   gender generation  device_type log_month

1       1  game-01 2013-09-06   2013-04-15      M         40         iOS   2013-09

2       1  game-01 2013-09-05   2013-04-15      M         40         iOS   2013-09

3       1  game-01 2013-09-28   2013-04-15      M         40         iOS   2013-09

4       1  game-01 2013-09-12   2013-04-15      M         40         iOS   2013-09

5       1  game-01 2013-09-11   2013-04-15      M         40         iOS   2013-09

6       1  game-01 2013-09-08   2013-04-15      M         40         iOS   2013-09

4. 성별로 집계

> table(dau.user.info[,c('log_month','gender')])

gender

log_month     F     M

2013-08 47343 46842

2013-09 38027 38148

> table(dau.user.info[,c('gender','log_month')])

log_month

gender 2013-08 2013-09

F   47343   38027

M   46842   38148

-> 전체적인 수치는 떨어졌지만 남녀간의 구성비율은 변화가 없으므로, 성별과 게임이용율 간 연관관계가 크지 않음.

5. 연령별로 집계

> table(dau.user.info[,c('log_month','generation')])

generation

log_month    10    20    30    40    50

2013-08 18785 33671 28072  8828  4829

2013-09 15391 27229 22226  7494  3835

> table(dau.user.info[,c('generation','log_month')])

log_month

generation 2013-08 2013-09

10   18785   15391

20   33671   27229

30   28072   22226

40    8828    7494

50    4829    3835

-> 전체적인 수치는 떨어졌지만 연령간의 구성비율은 변화가 없으므로, 연령과 게임이용율 간 연관관계가 크지 않음.

6. 세그먼트 분석(성별과 연령대를 조합해서 집계)

> library(reshape2)

> dcast(dau.user.info,log_month~gender+generation,value.var='user_id',length)

log_month F_10  F_20  F_30 F_40 F_50 M_10  M_20  M_30 M_40 M_50

1   2013-08 9091 17181 14217 4597 2257 9694 16490 13855 4231 2572

2   2013-09 7316 13616 11458 3856 1781 8075 13613 10768 3638 2054

-> 전체 이용율 하락에 영햐을 끼친 세그먼트는 찾아내지 못함

7. 세그먼트 분석(단말기별로 집계)

> head(dau.user.info)

user_id app_name   log_date   install_date  gender  generation device_type log_month

1       1  game-01 2013-09-06   2013-04-15      M         40         iOS   2013-09

2       1  game-01 2013-09-05   2013-04-15      M         40         iOS   2013-09

3       1  game-01 2013-09-28   2013-04-15      M         40         iOS   2013-09

4       1  game-01 2013-09-12   2013-04-15      M         40         iOS   2013-09

5       1  game-01 2013-09-11   2013-04-15      M         40         iOS   2013-09

6       1  game-01 2013-09-08   2013-04-15      M         40         iOS   2013-09

> table(dau.user.info[,c('log_month','device_type')])

device_type

log_month Android   iOS

2013-08   46974  47211

2013-09   29647  46528

-> ios의 변화율은 없는 반면, android의 변화율이 큼을 알 수 있다.

8. 세그먼트 분석 결과 시각화하기

#날짜별로 단말기별 유저수 산출하기

> library(plyr)

> dau.user.info.device.summary<-ddply(dau.user.info,

.(log_date,device_type),

                                                   summarize,

                                                    dau=length(user_id))

> head(dau.user.info.device.summary)

log_date  device_type  dau

1 2013-08-01     Android 1784

2 2013-08-01         iOS   1805

3 2013-08-02     Android 1386

4 2013-08-02         iOS   1451

5 2013-08-03     Android 1295

6 2013-08-03         iOS   1351

#날짜별 데이터 형식으로 변환하기

> str(dau.user.info.device.summary)

'data.frame':          122 obs. of  3 variables:

$ log_date   : chr  "2013-08-01" "2013-08-01" "2013-08-02" "2013-08-02" ...

$ device_type: chr  "Android" "iOS" "Android" "iOS" ...

$ dau        : int  1784 1805 1386 1451 1295 1351 1283 1314 2002 2038 ...

> dau.user.info.device.summary$log_date<-as.Date(dau.user.info.device.summary$log_date)

> str(dau.user.info.device.summary)

'data.frame':  122 obs. of  3 variables:

$ log_date   : Date, format: "2013-08-01" "2013-08-01" "2013-08-02" ...

$ device_type: chr  "Android" "iOS" "Android" "iOS" ...

$ dau        : int  1784 1805 1386 1451 1295 1351 1283 1314 2002 2038 ... 

9. 시계열 트렌드 그래프 그리기

> library(ggplot2)

> library(scales)

> limits<-c(0,max(dau.user.info.device.summary$dau))

> ggplot(dau.user.info.device.summary,aes(x=log_date,y=dau,col=device_type,lty=device_type,shape=device_type))+ geom_line(lwd=1)+geom_point(size=4)+scale_y_continuous(label=comma,limits = limits)

-> 9월부터 안드로이드 유저가 빠르게 감소함을 알 수 있다. 따라서 안드로이드 시스템을 점검해야 한다는 결론을 얻을 수 있다.

#ggplot 그래프 구조 이해하기

> library(ggplot2)

> library(scales)

> head(dau.user.info.device.summary)

log_date device_type  dau

1 2013-08-01     Android 1784

2 2013-08-01         iOS 1805

3 2013-08-02     Android 1386

4 2013-08-02         iOS 1451

5 2013-08-03     Android 1295

6 2013-08-03         iOS 1351

> ggplot(dau.user.info.device.summary)

> ggplot(dau.user.info.device.summary,aes(x=log_date,y=dau,col=device_type,lty=device_type,shape=device_type))

> ggplot(dau.user.info.device.summary,aes(x=log_date,y=dau,col=device_type,lty=device_type,shape=device_type))+geom_line(lwd=1)

> ggplot(dau.user.info.device.summary,aes(x=log_date,y=dau,col=device_type,lty=device_type,shape=device_type))+geom_line(lwd=1)+geom_point(size=4)

> ggplot(dau.user.info.device.summary,aes(x=log_date,y=dau,col=device_type,lty=device_type,shape=device_type))+geom_line(lwd=1)+geom_point(size=4)+scale_y_continuous(label=comma)

> limits<-c(0,max(dau.user.info.device.summary$dau))

> max(dau.user.info.device.summary$dau)

[1] 2133

> ggplot(dau.user.info.device.summary,aes(x=log_date,y=dau,col=device_type,lty=device_type,shape=device_type))+geom_line(lwd=1)+geom_point(size=4)+scale_y_continuous(label=comma,limits = limits)

출처: 비지니스 활용 사레로 배우는 데이터 분석: R (사카마키 류지, 사토 요헤이 지음)

6.5 준모수적 방법

Cox의 비례 모형

모수적 모형은 가정이 타당할 때는 상당히 효율적이지만 어떤 모형이 적당한가에 대한 지식이 없다면 함부로 사용하기가 곤란하다. 이에 반해 Cox(1972)의 비례 위험 모형(proportional hazards model)은 준모수적(semi-parametric)방법으로서 생존 시간의 분포에 대한 가정을 필요로 하지 않는다. 또한 시간에 따라 바뀌는 공변량(time-dependent variable)의 경우에도 분석할 수 있다는 장점이 있어, 생존자료의 분석에 매우 자주 사용된다.

[예] Prenctice(1973)에 소개된 것으로 40명의 폐암 환자의 생존시간을 조사한 것이다. 40명의 환자 중 21명은 기존 치료방법인 처리1에, 나머지 19명은 새로운 치료방법은 처리2에 할당되었으며, 생존시간에 영향을 미칠 것으로 생각되는 공변량은 다음과 같다.

X1: 진단시의 환자상태(Perfermance Status : 0~100 점)

X2: 환자의 나이(단위 : 년)

X3: 진단 후 연구 참여시까지의 시간(단위 : 월)

Trt: 치료 방법 (1 – 기존 치료, 2 – 새로운 치료)

TYPE: 종양의 유형 (squamous, small, adeno, large)

1. 데이터 입력

> setwd('c:/Rwork/')

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

> colnames(lung)<-c('time','status','x1','x2','x3','trt','type')

> head(lung)

  time status x1 x2 x3 trt type

1  411      1 70 64  5   1    1

2  126      1 60 63  9   1    1

3  118      1 70 65 11   1    1

4   82      1 40 69 10   1    1

5    8      1 40 63 58   1    1

6   25      0 70 48  9   1    1

2. Cox비례위험모형에 근거한 생존시간 분석

> library(survival)

> coxfit1 = coxph(Surv(time,status)~x1+x2+x3+factor(trt)+factor(type),data=lung)

> summary(coxfit1)

Call:

coxph(formula = Surv(time, status) ~ x1 + x2 + x3 + factor(trt) +

    factor(type), data = lung)

  n= 40, number of events= 37

                   coef exp(coef)  se(coef)      z Pr(>|z|)   

x1            -0.060281  0.941500  0.013777 -4.375 1.21e-05 ***

x2            -0.015086  0.985027  0.022340 -0.675   0.4995   

x3             0.001201  1.001201  0.011886  0.101   0.9195   

factor(trt)2  -0.448171  0.638795  0.431302 -1.039   0.2988   

factor(type)2  0.279682  1.322709  0.547259  0.511   0.6093   

factor(type)3  1.418190  4.129638  0.625283  2.268   0.0233 * 

factor(type)4  0.361145  1.434971  0.479210  0.754   0.4511   

---

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

              exp(coef) exp(-coef) lower .95 upper .95

x1               0.9415     1.0621    0.9164    0.9673

x2               0.9850     1.0152    0.9428    1.0291

x3               1.0012     0.9988    0.9781    1.0248

factor(trt)2     0.6388     1.5654    0.2743    1.4876

factor(type)2    1.3227     0.7560    0.4525    3.8663

factor(type)3    4.1296     0.2422    1.2125   14.0655

factor(type)4    1.4350     0.6969    0.5610    3.6707

Concordance= 0.764  (se = 0.058 )

Rsquare= 0.524   (max possible= 0.994 )

Likelihood ratio test= 29.66  on 7 df,   p=0.0001097

Wald test            = 26.29  on 7 df,   p=0.0004479

Score (logrank) test = 30.67  on 7 df,   p=7.138e-05

- coxph: cox의 비례위험 모형을 ㅈ거합함.

-> 결과 해석:

- 진단시의 환자상태(x1)에 대한 회귀계수가 -0.060281로 유의하다(p-value < 0.001). 위험률(exp(coef))를 보면 진단시의 환자상태에 대한 점수가 1점 더 높은 환자는 다른 조건이 동일한 환자와 비교할 때 0.9415(exp(-0.060281))배 위험하다는 것을 알 수 있다.

exp(-coef) = 1.0621은 exp(coef)의 역수인 1/0.9415이다.

- 나이(x2), 진단 후 연구 참여시까지의 시간(x3), 처리(trt)는 모두 p-value가 0.05보다 커서 유의한 영향을 미치지 않는다.

- TYPE: 종양의 유형 (squamous, small, adeno, large)을 보면 factor(type)3인 adeno의 경우 squamous와 차이를 보인다.(p-value: 0.0233)

회귀계수가 0이라는 가설에 대한 검정통계량. 차례로 우도비 검정(-2 LOG L), Wald 검정통계량, Score 검정이 있다. 세 통계랑 모두 가설을 기각하게 되므로 (p-value < 0.05), 적어도 하나의 공변량은 의미가 있다는 것(= 적어도 모든 변수가 0이라는 귀무가설을 기각)을 알 수 있다.

Likelihood ratio test= 29.66  on 7 df,   p=0.0001097

Wald test            = 26.29  on 7 df,   p=0.0004479

Score (logrank) test = 30.67  on 7 df,   p=7.138e-05

3. 생존함수추정치

> fit4 = survfit(coxfit1)

> summary(fit4)

Call: survfit(formula = coxfit1)

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

    1     40       1 9.90e-01 1.10e-02     9.68e-01        1.000

    2     39       1 9.79e-01 1.67e-02     9.47e-01        1.000

    8     38       2 9.54e-01 2.71e-02     9.03e-01        1.000

   10     36       1 9.37e-01 3.32e-02     8.74e-01        1.000

   11     35       1 9.20e-01 3.87e-02     8.47e-01        0.999

   12     34       2 8.83e-01 4.84e-02     7.93e-01        0.984

   15     32       1 8.64e-01 5.33e-02     7.65e-01        0.975

   16     31       1 8.44e-01 5.79e-02     7.38e-01        0.965

   18     30       1 8.22e-01 6.27e-02     7.07e-01        0.954

   19     29       1 7.96e-01 6.77e-02     6.74e-01        0.941

   20     28       1 7.67e-01 7.27e-02     6.37e-01        0.924

   21     27       1 7.34e-01 7.73e-02     5.97e-01        0.903

   43     25       1 6.97e-01 8.22e-02     5.54e-01        0.879

   44     24       1 6.61e-01 8.62e-02     5.12e-01        0.854

   51     23       1 6.26e-01 8.94e-02     4.73e-01        0.828

   54     22       1 5.84e-01 9.23e-02     4.28e-01        0.796

   56     21       1 5.44e-01 9.39e-02     3.87e-01        0.763

   82     20       1 5.06e-01 9.46e-02     3.50e-01        0.730

   84     19       1 4.64e-01 9.54e-02     3.11e-01        0.694

   90     18       1 4.25e-01 9.51e-02     2.74e-01        0.659

  100     17       1 3.81e-01 9.46e-02     2.34e-01        0.620

  118     15       1 3.35e-01 9.36e-02     1.94e-01        0.580

  126     14       1 2.93e-01 9.12e-02     1.59e-01        0.540

  153     13       1 2.53e-01 8.79e-02     1.28e-01        0.500

  164     12       1 2.14e-01 8.37e-02     9.97e-02        0.461

  177     11       1 1.80e-01 7.81e-02     7.66e-02        0.421

  200     10       1 1.40e-01 7.08e-02     5.20e-02        0.377

  201      9       1 1.07e-01 6.19e-02     3.41e-02        0.333

  231      8       1 8.00e-02 5.25e-02     2.21e-02        0.289

  250      6       1 5.18e-02 4.19e-02     1.06e-02        0.253

  287      5       1 2.88e-02 2.98e-02     3.81e-03        0.218

  340      4       1 9.18e-03 1.50e-02     3.74e-04        0.225

  411      3       1 2.19e-03 5.10e-03     2.30e-05        0.210

  991      2       1 1.28e-04 5.20e-04     4.54e-08        0.362

  999      1       1 3.36e-10 5.39e-09     7.43e-24        1.000

 > names(fit4)

 [1] "n"         "time"      "n.risk"    "n.event"   "n.censor"  "surv"     

 [7] "type"      "cumhaz"    "std.err"   "upper"     "lower"     "conf.type"

[13] "conf.int"  "call"     

> fit4$surv

 [1] 9.895862e-01 9.787704e-01 9.543895e-01 9.372414e-01 9.195316e-01

 [6] 8.834541e-01 8.637107e-01 8.439345e-01 8.215162e-01 7.962148e-01

[11] 7.670584e-01 7.344065e-01 7.344065e-01 6.974181e-01 6.609451e-01

[16] 6.256915e-01 5.835925e-01 5.436422e-01 5.055202e-01 4.643798e-01

[21] 4.249339e-01 3.809762e-01 3.809762e-01 3.352450e-01 2.932721e-01

[26] 2.533265e-01 2.142532e-01 1.795911e-01 1.401508e-01 1.065182e-01

[31] 8.001270e-02 5.176743e-02 2.881620e-02 9.176925e-03 2.194374e-03

[36] 1.281990e-04 3.359491e-10

4. 생존함수그래프

> plot(survfit(coxfit1),xlab='time',ylab='Survival function',xlim=c(0,998.9))

> legend(500,1.0,c('누적한계추정치','95%신뢰구간'),lty=c(1,2))

5. 누적함수그래프

누적위험함수의 추정치 그래프를 통해 위험함수의 형태를 짐작할 수 있다. 예를 들어 단조 증가하는 직선형태는 위험함수가 시간에 대해 일정하다는 것을 의미하며, 위쪽으로 휘는 모양이면 시간이 지남에 따라 일정하다는 것을 의미하며, 위쪽으로 휘는 모양이면 시간이 지남에 따라 위험함수가 증가하고, 아래 방향으로 휘면 감소한다는 것을 의미한다.

> H.hat = -log(fit4$surv)

> H.hat = c(H.hat,tail(H.hat,1)) 

> plot(c(fit4$time,1100),H.hat,xlab='time',ylab='comulative hazard function',type='s')

 - tail: 벡터, 매트릭스 데이터에서 마지막 n개의 행들을 선택함. 예를 들어 tail(H.hat,1)라고 입력하면 H.hat 벡터의 마지막 1개의 성분을 선택. 해당 처리를 하는 이유는 위험함수 곡선은 시간에 대한 상승곡선이므로, 마지막 함수값을 추가함을써 곡선의 모양을 자연스럽게 하기 위한 처리임.

6. 비례성 검토를 위한 로그-로그 그림

비례성의 가정이 타당한 것인지 검토하는 방법에 대해 알아보자. 폐암자료에서 두 처리그룹별로 시간에 따른 log(-log S^ (t)) 의 그래프를 그렸을 때 평행하게 되는 가를 볼 수 있다.

> coxfit2 = coxph(Surv(time,status)~x1+x2+x3+strata(trt)+factor(type),data=lung) 

> plot(survfit(coxfit2),fun='cloglog',lty=1:2,col=c('red','blue'))#fun='cloglog', 두 그래프가 대체적으로 평행하므로 비례성이 타당하다고 할 수 있다.

> legend('topleft',c('처리1','처리2'),lty=1:2,col=c('red','blue'))

- strata(trt): 처리를 층으로 입력해주면, 처리를 층으로 입력하여, 처리1과 처리2로 나누어 Cox 비례 위험 모형을 추정함.

- fun='cloglog'를 입력하면 로그-로그 그림을 그릴 수 있다.

- 비례성검토를 위한 로그-로그 그림(log-log plot)은 그래프가 평행하므로 비례성의 가정이 타당하다 할 수 있다.

7. 로그 – 랭크 테스트 log – rank test

1) 종양 유형 별 생존함수의 차이

> survdiff(Surv(time,status)~factor(type),data=lung)

Call:

survdiff(formula = Surv(time, status) ~ factor(type), data = lung)

                N Observed Expected (O-E)^2/E (O-E)^2/V

factor(type)=1 14       12    16.72     1.330     2.825

factor(type)=2 11       10     6.93     1.356     1.771

factor(type)=3  5        5     2.18     3.651     4.099

factor(type)=4 10       10    11.17     0.123     0.187

 Chisq= 7.4  on 3 degrees of freedom, p= 0.0614

X²(df=3)=7.4이고, p-value가 0.0614로 0.05보다는 약간 크지만 상당히 의미 있음을 알 수 있다.

2) 치료법 생존함수의 차이

> survdiff(Surv(time,status)~x1,data=lung)

Call:

survdiff(formula = Surv(time, status) ~ x1, data = lung)

      N Observed Expected (O-E)^2/E (O-E)^2/V

x1=20 2        2    0.128   27.3119   28.6351

x1=30 4        4    1.565    3.7893    4.2568

x1=40 7        7    1.932   13.2948   15.1966

x1=50 4        3    3.298    0.0270    0.0305

x1=60 7        7    6.642    0.0193    0.0247

x1=70 9        7   12.420    2.3650    3.7914

x1=80 6        6    6.982    0.1382    0.1761

x1=90 1        1    4.033    2.2811    3.7772

 Chisq= 58.8  on 7 degrees of freedom, p= 2.65e-10

X²(df=7)=58.8이고, p-value가 0.05보다 매우 작기 때문에 의미 있음을 알 수 있다.

8. anova를 이용하여 모형 비교

> coxfit1 = coxph(Surv(time,status)~x1+x2+x3+factor(trt)+factor(type),data=lung)

> coxfit2 = coxph(Surv(time,status)~x1+factor(type),data=lung)

> anova(coxfit2,coxfit1)

Analysis of Deviance Table

 Cox model: response is  Surv(time, status)

 Model 1: ~ x1 + factor(type)

 Model 2: ~ x1 + x2 + x3 + factor(trt) + factor(type)

   loglik Chisq Df P(>|Chi|)

1 -88.143                  

2 -87.516 1.254  3    0.7401

anova()는 Cox Regression의 모형을 비교할 때 LRT(Likelihood ratio test)를 사용한다. 환자의 나이(x2), 진단 후 연구 참여시까지의 시간( x3)의 p-value는 0.7401으로 유의하지 않다.

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

6.4 모수적 방법을 이용한 생존함수의 추정과 비교

공학(시멘트의 양, 유리의 버티는 힘), 경영(고객 수), 교통(소방차 수) 모두 모수적 방법을 이용.

분포를 이루기 때문에 많은 분야에서 사용된다.

1. 데이터 입력

> setwd('c:/Rwork')

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

> head(lung)

  癤퓍ime status x1 x2 x3 trt type

1     411      1 70 64  5   1    1

2     126      1 60 63  9   1    1

3     118      1 70 65 11   1    1

4      82      1 40 69 10   1    1

5       8      1 40 63 58   1    1

6      25      0 70 48  9   1    1

> colnames(lung)<-c('time','status','x1','x2','x3','trt','type')

> head(lung)

  time status x1 x2 x3 trt type

1  411      1 70 64  5   1    1

2  126      1 60 63  9   1    1

3  118      1 70 65 11   1    1

4   82      1 40 69 10   1    1

5    8      1 40 63 58   1    1

6   25      0 70 48  9   1    1

> attach(lung)

2. 모수적 모형에 근거한 생존시간 분석

> library(survival)

> weibull = survreg(Surv(time,status)~x1+x2+x3+factor(trt)+factor(type),data=lung,dist='weibull') #공변량factor(trt), 처리효과factor(type), 와이블 분포weibull, gaussian(정규분포), logistic(로지스틱)eh rksmd.

> summary(weibull)

Call:

survreg(formula = Surv(time, status) ~ x1 + x2 + x3 + factor(trt) +

    factor(type), data = lung, dist = "weibull")

                 Value Std. Error      z        p

(Intercept)    1.06044    1.35959  0.780 4.35e-01

x1             0.05420    0.00954  5.680 1.35e-08

x2             0.01168    0.01918  0.609 5.42e-01

x3             0.00379    0.01051  0.361 7.18e-01

factor(trt)2   0.28871    0.36899  0.782 4.34e-01

factor(type)2 -0.49964    0.45323 -1.102 2.70e-01

factor(type)3 -1.25968    0.49732 -2.533 1.13e-02

factor(type)4 -0.40243    0.38726 -1.039 2.99e-01

Log(scale)    -0.13615    0.13146 -1.036 3.00e-01

Scale= 0.873

Weibull distribution

Loglik(model)= -203.4   Loglik(intercept only)= -219.7

        Chisq= 32.59 on 7 degrees of freedom, p= 3.2e-05

Number of Newton-Raphson Iterations: 6

n= 40

> weibull

Call:

survreg(formula = Surv(time, status) ~ x1 + x2 + x3 + factor(trt) +

    factor(type), data = lung, dist = "weibull")

Coefficients:

  (Intercept)            x1            x2            x3  factor(trt)2

  1.060436846   0.054195931   0.011681287   0.003792838   0.288708242

factor(type)2 factor(type)3 factor(type)4

 -0.499640589  -1.259681146  -0.402431957

Scale= 0.8727099

Loglik(model)= -203.4   Loglik(intercept only)= -219.7

        Chisq= 32.59 on 7 degrees of freedom, p= 3.2e-05

n= 40

logT = 1.060 + 0.054x1 + 0.012x2 + 0.004x3 + 0.289x4 - 0.500x5 -1.260x6 - 0.402x7 +0.873e

factor(trt)2   0.28871    0.36899  0.782 4.34e-01

-> x4는 처리그룹을 나타내는 가변수(dummy variable)로서 처리가 standard일 때 1, test 일 때 2의 값을 갖는다.

p-value가 0.434로 유의하지 않으므로 두 처리(standard, test)는 생존시간의 차이를 보이지 않는다.

type은  squamous,small,adeno,large, x5,6,7은 종양의 유형을 가르키는 가변수로서

x5 = 1, if 'small', = 0 o.w.

x6 = 1, if 'adeno', = 0 o.w.

x7 = 1, if 'large', = 0 o.w.

따라서 종양이 squamous인 경우에는 x5=x6=x7 = 0 된다.

factor(type)3 -1.25968    0.49732 -2.533 1.13e-02

-> adeno의 p-value만 0.0113으로 0.05보다 적으므로 유의, squamous와 차이를 보인다고 할 수 있다. adeno의 회귀계수가 -1.25968로 이 유형을 가진 사람들은 상대적으로 생존시간이 짧다는 것을 알 수 있다. type2 small, type4 large의 경우  p-value가 각각 0.270, 0,299로 유의하지 않다, 즉 squamous와 차이를 보이지 않는다.

Loglik(model)= -203.4   Loglik(intercept only)= -219.7

        Chisq= 32.59 on 7 degrees of freedom, p= 3.2e-05

2(logL-logL0) = 32.59 > 14.067 모든 공변량의 회귀계수가 0이라는 귀무가설을 아주 강하게 기각한다.

통상적인 선형모형에서 모형에 대한 F-검정을 하는 것과 같이 여기서도 모든 공변량의 회귀계수가 0이라는 가설에 대하 우도비 검정(likelihood-ratio ttest)를 할 수 있다. 우리가 고려한 모형에서의 로그-우도(log-likelihood)를 logL, 공변량이 전혀 없는 귀무모형에서의 로그-우도를 logL0라고 했을 때 2(logL-logL0)이 귀무가설하에서 근사적으로 x2-분포를 따르게 되므로 이를 이용하여 검정할 수 있다.

이때 자유도는 귀무모형에서 제외되는 공변량의 수와 같게 되며, 위에서는 7이 된다.

LogL: 고려한 모형에서의 로그-우도 Log Likelihood for WEIBULL -203.4는 이 모형에 대한 로그-우도이다.

LogL0: 공변량이 전혀 없는 귀무모형에 대한 로그-우도를 구해보면 Log Likelihood for WEIULL = =219.7을

따라서 Chisq= 32.59는 아래와 같은 계산을 통해 얻을 수 있으며, 이 값이 x2-분포의 임계치인 x2 0.95(7) = 14.067보다 크므로

5%유의수준하에서 모든 공변량의 회귀계수가 0이라는 귀무사설을 기각하게 된다.

> 2*(-203.4-(-219.7))

[1] 32.6

3. 모수적 모형의 적합도 검토

고려한 모형이 타당한가를 검토하는 방법 중 하나는 로그-우도 비교. R을 이용하여 각 모형에 대하 로그-우도를 출력한 다음 이들을 비교하여 절대값이 가장 큰 모형을 택할 수 있다. 또는 AIC(Akaike information criterion) 값을 비교하여 이 값이 더 작은 것을 선택할 수 있다.

> library(flexsurv) #flexsurv library: Flexible parametric survival models

Warning message:

package ‘flexsurv’ was built under R version 3.2.5

> gengamma=flexsurvreg(formula=Surv(time,status)~x1+x2+x3+factor(trt)+factor(type),data=lung,dist='gengamma')

> gengamma

Call:

flexsurvreg(formula = Surv(time, status) ~ x1 + x2 + x3 + factor(trt) +

    factor(type), data = lung, dist = "gengamma")

Estimates:

               data mean  est       L95%      U95%      se        exp(est)

mu                   NA    1.16226  -1.69652   4.02104   1.45859        NA

sigma                NA    0.82367   0.47903   1.41628   0.22778        NA

Q                    NA    1.19119  -0.35449   2.73687   0.78863        NA

x1             56.50000    0.05426   0.03598   0.07254   0.00933   1.05576

x2             56.57500    0.01210  -0.02543   0.04963   0.01915   1.01217

x3             15.65000    0.00494  -0.01741   0.02729   0.01141   1.00495

factor(trt)2    0.50000    0.27185  -0.47815   1.02186   0.38266   1.31239

factor(type)2   0.27500   -0.57430  -1.63425   0.48566   0.54080   0.56310

factor(type)3   0.12500   -1.36101  -2.57769  -0.14434   0.62076   0.25640

factor(type)4   0.25000   -0.48503  -1.45929   0.48923   0.49708   0.61568

               L95%      U95%   

mu                   NA        NA

sigma                NA        NA

Q                    NA        NA

x1              1.03664   1.07524

x2              0.97489   1.05088

x3              0.98274   1.02767

factor(trt)2    0.61993   2.77835

factor(type)2   0.19510   1.62524

factor(type)3   0.07595   0.86560

factor(type)4   0.23240   1.63107

N = 40,  Events: 37,  Censored: 3

Total time at risk: 5784

Log-likelihood = -203.4059, df = 10

AIC = 426.8117

> weibull=flexsurvreg(formula=Surv(time,status)~x1+x2+x3+factor(trt)+factor(type),data=lung,dist='weibull')

> weibull

Call:

flexsurvreg(formula = Surv(time, status) ~ x1 + x2 + x3 + factor(trt) +

    factor(type), data = lung, dist = "weibull")

Estimates:

               data mean  est       L95%      U95%      se        exp(est)

shape                NA    1.14586   0.88555   1.48269   0.15066        NA

scale                NA    2.88763   0.20141  41.39933   3.92317        NA

x1             56.50000    0.05420   0.03551   0.07288   0.00953   1.05569

x2             56.57500    0.01168  -0.02587   0.04924   0.01916   1.01175

x3             15.65000    0.00379  -0.01679   0.02438   0.01050   1.00380

factor(trt)2    0.50000    0.28871  -0.43443   1.01185   0.36896   1.33470

factor(type)2   0.27500   -0.49964  -1.38783   0.38855   0.45317   0.60675

factor(type)3   0.12500   -1.25968  -2.23430  -0.28506   0.49726   0.28374

factor(type)4   0.25000   -0.40243  -1.16137   0.35651   0.38722   0.66869

               L95%      U95%   

shape                NA        NA

scale                NA        NA

x1              1.03615   1.07560

x2              0.97446   1.05047

x3              0.98335   1.02468

factor(trt)2    0.64763   2.75068

factor(type)2   0.24962   1.47484

factor(type)3   0.10707   0.75197

factor(type)4   0.31306   1.42833

N = 40,  Events: 37,  Censored: 3

Total time at risk: 5784

Log-likelihood = -203.4363, df = 9

AIC = 424.8727

일반화감마분포 gengamma분포의 경우 로그-우도 값이 -203.4059, AIC는 426.8117

와이블분포 wibull 분포의 경우 로그-우도 값이 -203.4363, AIC는 424.8727로 큰 차이는 없다. 그래프로 그려보면 아래와 같이 큰 차이가 없음을 알 수 있다.

데이터 시각화

> plot(weibull,xlab='time',ylab='survival function',ci=F)

> lines(gengamma,col='blue',lty=2,ci=F)

> legend('topright',c('weibull','generalized gamma'),lty=1:2,col=c('red','blue'))

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

6.3 비모수적 방법을 이용한 생존함수의 비교

[예] 흑색종(melanoma) 환자들에 대한 BCG와 CP(coryne-bacterium parvum)의 생존지속 효과를 비교하기 위한 연구에서 30명의 흑색종 환자 중 11명은 BCG 처리를 받고 나머지 19명은 CP처리를 받았다고 한다. 중도절단이 포함되어 있는 경우에 이 두 그룹의 생존분포를 비교하기 위한 방법을 알아보자.