데이터마이너를 꿈꾸며

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

> 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 (사카마키 류지, 사토 요헤이 지음)