데이터마이너를 꿈꾸며

제목 그대로 R vs 파이썬, 무엇이 더 좋을까?

나도 파이썬을 배우기 전에는 R을 먼저 공부하였다. 2016년으로 기억한다. 당시는 방통대 통계학에 편입하여 공부하고 있었는데, 교재에서 많은 데이터들이 R을 다루기도 하였고, 인터넷으로 검색을 해보니 R 강의도 많았다.

파이썬이 없었던 것은 아니였지만, 처음 교재로 R을 접하다보니 자연스롭게 노출된 언어를 먼저 공부하게 되었다. 당시에도 R이냐 파이썬이냐를 두고 갑론을박을 했던 기억이 난다. 기억으로는 통계를 할거면 R, 통계를 넘어 개발까지 고민한다면 파이썬을 공부하라는 전문가의 말도 있었고, 어떤 프로그램 언어에 채용 공고가 많은지 보라는 의견도 있었다. 실제로 지금도 해당 내용은 인터넷 검색이 쉽게 되는데 당연하지만 파이썬의 채용공고가 더 많다. 그리고 캐글. 캐글은 데이터쟁이들의 놀이터인데 압도적인 비율로 파이썬을 사용한 결과값을 많이 제출하였고, R만 알던 나에게 파이썬도 해야겠구나...라고 생각하게 만들었던 생각이 난다. 

문득 챗 GPT는 뭐라고 말할지 궁금해졌다. 두 프로그램 언어 중 무엇이 좋은지 물어보았다.

R의 장점으로는 통게 및 데이터 분석에 최적화 되있다는 점이다. 통계적 가설 검정을 많이 했던 기억이 난다. 그리고 일본 저자가 쓴 교재를 통해 시각화 공부도 많이 했었다. 파이썬은 역시나 데이터 분석 외 이야기들이 많이 나온다. R만으로 머신러닝이 되지 않는건 아닌데, 파이썬은 R에서 제공하지 못하는 많은 분석 외적 영역을 제공해주곤 하였다. 그리고 이게 지금 내가 파이썬을 공부하는 가장 큰 이유인데 유저가 많다는 것이다. 당연히 구글 검색을 하면 관련된 결과가 더 많이 나와 공부하기 쉬운 환경이다.

물어보다 파이썬으로 했을 때 R 대비 장점을 한번 더 물어보았다.

1. 자료 가져오기 

2. 다양한 형태의 자료 가져오기

3. 통계 패키지 자료 가져오기

 - 카이제곱 = 27.458 , p - value = 1.6 X 10마이너스 7승. 자극과 반응간의 유의한 관계가 있음.

4. RData 저장 및 가져오기

참고문헌: 

[1] 고급 R 활용 (심송용, 이윤동, 이은경, 김성수 공저, KNOU PRESS)

#데이터 불러오기

getwd()[1] "C:/Users/amore/Documents/FOR ME/Data Scientist/Udacity/Data Analysis with R/eda-course-materials/lesson3"
list.files()
[1] "lesson3.Rmd"         "lesson3_student.rmd" "pseudo_facebook.tsv"
pf <- read.csv('pseudo_facebook.tsv', sep='\t')
names(pf)
 [1] "userid"                "age"                   "dob_day"              
 [4] "dob_year"              "dob_month"             "gender"               
 [7] "tenure"                "friend_count"          "friendships_initiated"
[10] "likes"                 "likes_received"        "mobile_likes"         
[13] "mobile_likes_received" "www_likes"             "www_likes_received"   

#Histogram of User's Birthdays

library(ggplot2)qplot(x=dob_day, data = pf) 

Here's some things that I noticed. On the first day of the month I see this huge bin of almost 8,000 people. This seems really unusual since I would expect most people to have the same number of birthday's across every day of the month.  

library(ggplot2)qplot(x=dob_day, data = pf) + 
  scale_x_continuous(breaks=1:31)

#Faceting

library(ggplot2)qplot(x=dob_day, data = pf) + 
  scale_x_continuous(breaks=1:31)+
  facet_wrap(~dob_month, ncol=4) 

Now, you may have noticed some peaks in May or perhaps in October, but I think what's really interesting is this huge spike on January first. There's almost 4,000 users in this bin. Now, this could be because of the default settings that Facebook uses or perhaps users are choosing the first choice in the drop down menus. Another idea is that some users may want to protect their privacy and so they just go with January first by default. Whatever the case may be, I think it's important that we make our considerations in the context of our data. We want to look out for these types of anomalies. 

#Friend Count

ggplot(aes(x = friend_count),data = pf) +
   geom_histogram()

#Limiting the Axes

qplot(x = friend_count, data = pf, xlim = c(0,1000))


## alternative solution
qplot(x = friend_count, data = pf) + 
  scale_x_continuous(limits = c(0,1000))

#adjust the bin width

 qplot(x = friend_count, data = pf) + 
   scale_x_continuous(limits = c(0,1000), breaks = seq(0,1000,50))

## alternative solution
qplot(x = friend_count, data = pf) + 
  scale_x_continuous(limits = c(0,1000))

##splits up the data by gender 

qplot(x = friend_count, data = pf) + 
  scale_x_continuous(limits = c(0,1000), breaks = seq(0,1000,50)) + 
  facet_grid(gender~.)




#Omitting NA Obervations
qplot(x = friend_count, data = subset(pf, !is.na(gender)), 
      binwidth = 10) + 
  scale_x_continuous(limits = c(0,1000), breaks = seq(0,1000,50)) + 
  facet_grid(gender~.)

#equivalent ggplot syntax:

ggplot(aes(x = friend_count), data = pf) +
  geom_histogram() +
  scale_x_continuous(limits = c(0, 1000), breaks = seq(0, 1000, 50)) +
  facet_wrap(~gender)




ggplot(aes(x = friend_count), data = subset(pf, !is.na(gender))) +
  geom_histogram() +
  scale_x_continuous(limits = c(0, 1000), breaks = seq(0, 1000, 50)) +
  facet_wrap(~gender)

##statistics by gender

> table(pf$gender)female   male 
 40254  58574 
> by(pf$friend_count, pf$gender, summary)
pf$gender: female
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
      0      37      96     242     244    4923 
--------------------------------------------------------------------- 
pf$gender: male
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
      0      27      74     165     182    4917 

### Tenure

qplot(x=tenure, data = pf, binwidth = 30,
      color = I('black'), fill = I('#099DD9'))



#Equivalent ggplot syntax:
ggplot(aes(x = tenure), data = pf) +
  geom_histogram(binwidth = 30, color = 'black', fill = '#099DD9')

## How would you create a histogram of tenure by year?

#create a histogram of tenure measured in years rather than in days





#It looks like the bulk of our users had less than two and a half years on Facebook.qplot(x = tenure / 365, data = pf, binwidth = 0.25, 
      color = I('black'), fill = I('#099DD9')) +
  scale_x_continuous(breaks = c(1,7,1), limits = c(0,7))
#Equivalent ggplot syntax:
ggplot(aes(x = tenure/365), data = pf) +
  geom_histogram(binwidth = .25, color = 'black', fill = '#F79420')

## Labeling Plots

qplot(x = tenure / 365, data = pf, binwidth = 0.25, 
      xlab = 'Number of years using Facebook',      ylab = 'Number of users in sample',
      color = I('black'), fill = I('#099DD9')) +
  scale_x_continuous(breaks = c(1,7,1), limits = c(0,7))


#Equivalent ggplot syntax:

ggplot(aes(x = tenure / 365), data = pf) +
  geom_histogram(color = 'black', fill = '#F79420') +
  scale_x_continuous(breaks = seq(1, 7, 1), limits = c(0, 7)) +
  xlab('Number of years using Facebook') +
  ylab('Number of users in sample')

## User Ages

> range(pf$age)
[1]  13 113
qplot(x = age, data = pf,binwidth = 5,
      xlab = 'User ages',
      ylab = 'Number of users in sample',
      color = I('black'), fill = I('blue')) +
  scale_x_continuous(breaks = c(0,113,5),limits= c(13,113))


#Equivalent ggplot syntax:

ggplot(aes(x = age), data = pf) +
  geom_histogram(binwidth = 1, fill = '#5760AB') +
  scale_x_continuous(breaks = seq(0, 113, 5))

## Transforming Data

> summary(pf$friend_count)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
    0.0    31.0    82.0   196.4   206.0  4923.0 
> summary(log10(pf$friend_count))
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   -Inf   1.491   1.914    -Inf   2.314   3.692 
> summary(log10(pf$friend_count + 1))
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  0.000   1.505   1.919   1.868   2.316   3.692 
> summary(sqrt(pf$friend_count))
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  0.000   5.568   9.055  11.090  14.350  70.160 
#Learn about how to use scales and how to create multiple plots on one page.  You'll also need to install and load the package gridExtra.
#Create 1 column thwi the folling thres histograms, friend count / transformed using  log 10 and square root.
> library(gridExtra)
> p1 <- qplot(x = friend_count, data = pf)
> p2 <- qplot(x = log10(friend_count+1), data = pf)
> p3 <- qplot(x = sqrt(friend_count), data = pf)
> grid.arrange(p1,p2,p3, ncol = 1)
#Equivalent ggplot syntax: we get the same outputs as we had before
pp1 <- ggplot(aes(x=friend_count), data = pf) + geom_histogram()
pp2 <- pp1 + scale_x_log10()
pp3 <- pp1 + scale_x_sqrt()

grid.arrange(pp1,pp2,pp3, ncol = 1)

## Add a Scaling Layer

> logScale <- qplot(x = log10(friend_count), data = pf)
> countScale <- ggplot(aes(x = friend_count), data = pf) +
  geom_histogram() +
  scale_x_log10()
> grid.arrange(logScale, countScale, ncol=2)

-> At the two plots, we can see that the difference is really in the x axis labeling. 
   Using scale_x_log10 will label the axis in actual friend_count. Where as using the 
   log10 wrapper will label the x axis in the log units. In general it is easier to 
   think about actual counts, so that's why people prefer using the scale_x_log10 as 
   a layer.

qplot(x=friend_count, data = pf) +
    scale_x_log10()

## Frequency Polygons : who has more friends on average men or women?

#This allows us to see the shape and the peaks of our distribution in more detail.qplot(x = friend_count, data = subset(pf, !is.na(gender)),
      binwidth = 10) +
  scale_x_continuous(lim = c(0,1000), breaks = seq(0,1000,50)) +
  facet_wrap(~gender)


qplot(x = friend_count, data = subset(pf, !is.na(gender)),
      binwidth = 10, geom='freqpoly', color = gender) +
  scale_x_continuous(lim = c(0,1000), breaks = seq(0,1000,50)) 



#But again, this plot doesn't really answer our question who has more friends on 
average men or women. Let's change the y-axis to show proportions instead of raw 
counts. 


qplot(x = friend_count, y = ..count.. / sum(..count..),
      data = subset(pf, !is.na(gender)),
      xlab = 'Friend Count',
      ylab = 'Proportion of Users with that firned count',
      binwidth = 10, geom='freqpoly', color = gender) +
  scale_x_continuous(lim = c(0,1000), breaks = seq(0,1000,50)) 

## Likes on the Web

qplot(x=www_likes, data = subset(pf, !is.na(gender)),
      geom = 'freqpoly', color = gender) +
  scale_x_continuous()+
  scale_x_log10()



# what's the www_like count for males?The first question is asking how many www_likes 
  there are in the entire data set for males.

> by(pf$www_likes, pf$gender, sum)

pf$gender: female
[1] 3507665
--------------------------------------------------------------------- 
pf$gender: male

## Box Plots

qplot(x = gender, y = friend_count,      data = subset(pf, !is.na(gender)),
      geom='boxplot')

## Adjust the code to focus on users who have friend counts between 0 and 1000.

qplot(x = gender, y = friend_count,
      data = subset(pf, !is.na(gender)),
      geom='boxplot',ylim = c(0,1000))




#same way
qplot(x = gender, y = friend_count,
      data = subset(pf, !is.na(gender)),
      geom = 'boxplot') +
scale_y_continuous(limits=c(0,1000))




-> Notice how the top of the box is just below 250, so it might be around 230. But this value might not be accurate for all of our data. use the ylim parameter or the scale_y_continious layer, we actually remove data points from calculations. 

# So a better way to do this is tu use the cord Cartesian layer to set the y limits instead.qplot(x = gender, y = friend_count,
      data = subset(pf, !is.na(gender)),
      geom = 'boxplot') +
  coord_cartesian(ylim= c(0,1000))



-> Here we will set the y limts from 0 to a 1000, notice how the top of the box has moved slightly closer to 250 for females. 

## Box Plots, Quartiles, and Friendships

qplot(x = gender, y = friend_count,
      data = subset(pf, !is.na(gender)),
      geom = 'boxplot') +
  coord_cartesian(ylim = c(0,250))


> by(pf$friend_count, pf$gender, summary)
pf$gender: female
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
      0      37      96     242     244    4923 
---------------------------------------------------------------------------- 
pf$gender: male
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
      0      27      74     165     182    4917 
-> The third quartile of the 75% mark is at 244 and that's all the way up 
here(그래프를 보며). This means that 75% of female users have friend count below 244.
Or another way to say this is that 25% of female user have more than 244 friends.

## On average, who initiated more friendships in our sample: men or women?

> names(pf)
 [1] "userid"                "age"                   "dob_day"              
 [4] "dob_year"              "dob_month"             "gender"               
 [7] "tenure"                "friend_count"          "friendships_initiated"
[10] "likes"                 "likes_received"        "mobile_likes"         
[13] "mobile_likes_received" "www_likes"             "www_likes_received"   
qplot(x = gender, y = friendships_initiated,
      data = subset(pf, !is.na(gender)), geom = 'boxplot') +
  coord_cartesian(ylim = c(0,150))

> by(pf$friendships_initiated, pf$gender, summary)
pf$gender: female
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
    0.0    19.0    49.0   113.9   124.8  3654.0 
---------------------------------------------------------------------------- 
pf$gender: male
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
    0.0    15.0    44.0   103.1   111.0  4144.0 

-> On average, who initiated more friendships in our sample: men or women?
Women.

## Getting Logical : What percent of check in using mobile?

> head(pf$mobile_likes)
[1] 0 0 0 0 0 0
> summary(pf$mobile_likes)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
    0.0     0.0     4.0   106.1    46.0 25110.0 
> summary(pf$mobile_likes > 0)
   Mode   FALSE    TRUE    NA's 
logical   35056   63947       0 
> mobile_check_in <- NA
> pf$mobile_check_in <- ifelse(pf$mobile_likes>0, 1, 0)
> pf$mobile_check_in <- factor(pf$mobile_check_in)
> summary(pf$mobile_check_in)
    0     1 
35056 63947 
> sum(pf$mobile_check_in == 1) / length(pf$mobile_check_in) 
[1] 0.6459097
#what percent of check in using mobile? 65%

#reference

1. https://cn.udacity.com/course/data-analysis-with-r--ud651

2. http://www.cookbook-r.com/Graphs/Facets_(ggplot2)

3. https://en.wikipedia.org/wiki/Web_colors

4. http://ggplot2.tidyverse.org/reference/theme.html

5. http://lightonphiri.org/blog/ggplot2-multiple-plots-in-one-graph-using-gridextra

6. http://ggplot2.tidyverse.org/reference/scale_continuous.html

7. https://en.wikipedia.org/wiki/Linear_regression#Assumptions

8. https://en.wikipedia.org/wiki/Normal_distribution

9. https://www.r-statistics.com/2013/05/log-transformations-for-skewed-and-wide-distributions-from-practical-data-science-with-r/

# 원하는 데이터만 발라서 보기

#reddit data 보기, 범주형 변수 그래프 그리기

#오더 순서 정하기 setting levles of ordered factors solution

#alternative solution

# practice

참고

https://cn.udacity.com/course/data-analysis-with-r--ud651

https://cn.udacity.com/course/data-wrangling-with-mongodb--ud032

http://vita.had.co.nz/papers/tidy-data.pdf

http://courses.had.co.nz.s3-website-us-east-1.amazonaws.com/12-rice-bdsi/slides/07-tidy-data.pdf

http://www.computerworld.com/article/2497143/business-intelligence/business-intelligence-beginner-s-guide-to-r-introduction.html

http://www.statmethods.net/index.html

https://www.r-bloggers.com/

http://www.cookbook-r.com/

http://blog.revolutionanalytics.com/2013/08/foodborne-chicago.html

http://blog.yhat.com/posts/roc-curves.html

https://github.com/corynissen/foodborne_classifier

문제인식: 매스미디어 광고에 의한 신규유저수가 일정치 않다. 이는 매월 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 (사카마키 류지, 사토 요헤이 지음)

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

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

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

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

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

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

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

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

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

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

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

> library(lsmeans)

> setwd('C:/Rwork')

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

> head(drug)

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

1        1      1             37       61 11.3208

2        2      2             37       37 12.9151

3        3      3             45       53 18.8947

4        4      4             41       41 14.6739

5        5      5             57       41  8.6493

6        6      6             49       33  9.5238

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

> colnames(drug)<-colname

> head(drug)

  obs type x1 x2       y

1   1    1 37 61 11.3208

2   2    2 37 37 12.9151

3   3    3 45 53 18.8947

4   4    4 41 41 14.6739

5   5    5 57 41  8.6493

6   6    6 49 33  9.5238

> attach(drug)

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

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

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

> summary(model1)

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

factor(type)     5  250.3    50.1   2.099   0.113   

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

x2               1   54.4    54.4   2.282   0.148   

factor(type):x1  5  160.8    32.2   1.349   0.289   

factor(type):x2  5   42.3     8.5   0.355   0.872   

Residuals       18  429.3    23.9                   

---

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

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

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

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

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

x2 소요시간 – 0.872

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

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

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

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

> summary(model2)

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

factor(type)  5  250.3    50.1   2.216   0.0808 . 

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

x2            1   54.4    54.4   2.409   0.1319   

Residuals    28  632.5    22.6                    

---

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

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

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

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

p-value: 0.0808

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

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

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

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

> summary(model3)

Call:

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

Residuals:

   Min     1Q Median     3Q    Max

-7.222 -2.634 -0.379  1.475  9.646

Coefficients:

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

(Intercept)    36.6370     8.9591   4.089 0.000331 ***

factor(type)2   0.8965     2.7823   0.322 0.749677   

factor(type)3   7.9097     2.7684   2.857 0.007973 **

factor(type)4   3.0722     2.8247   1.088 0.286025   

factor(type)5   9.5434     2.8534   3.345 0.002355 **

factor(type)6   5.8389     2.8391   2.057 0.049149 * 

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

x2              0.1647     0.1061   1.552 0.131868   

---

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

Residual standard error: 4.753 on 28 degrees of freedom

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

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

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

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

대립가설 not H0

p-value : 0.0001641

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

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

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

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

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

> summary(model2)

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

factor(type)  5  250.3    50.1   2.216   0.0808 . 

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

x2            1   54.4    54.4   2.409   0.1319   

Residuals    28  632.5    22.6                    

---

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

OR

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

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

factor(type)  5  250.3    50.1   2.216   0.0808 . 

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

x2            1   54.4    54.4   2.409   0.1319   

Residuals    28  632.5    22.6                    

---

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

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

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

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

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

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

x2            1   62.8    62.8   2.779   0.1067   

factor(type)  5  422.0    84.4   3.736   0.0102 * 

Residuals    28  632.5    22.6                    

---

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

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

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

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

> summary(model3)

Call:

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

Residuals:

   Min     1Q Median     3Q    Max

-7.222 -2.634 -0.379  1.475  9.646

Coefficients:

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

(Intercept)    36.6370     8.9591   4.089 0.000331 ***

factor(type)2   0.8965     2.7823   0.322 0.749677   

factor(type)3   7.9097     2.7684   2.857 0.007973 **

factor(type)4   3.0722     2.8247   1.088 0.286025   

factor(type)5   9.5434     2.8534   3.345 0.002355 **

factor(type)6   5.8389     2.8391   2.057 0.049149 * 

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

x2              0.1647     0.1061   1.552 0.131868   

---

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

Residual standard error: 4.753 on 28 degrees of freedom

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

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

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

ŷ_1j = β^₀ + α^₁ + β^1x_{11j  +} β^2x_21j  = 36.637 + 0 - 0.761 x_11j + 0.165x_21j

ŷ_2j = β^₀ + α^₂ + β^1x_{12j  +} β^2x_21j = 36.637 + 0.897- 0.761 x_12j + 0.165x_22j

ŷ_3j = β^₀ + α^₃ + β^1x_{13j  +} β^2x_21j = 36.637 + 7.908 - 0.761 x_13j + 0.165x_23j

ŷ_4j = β^₀ + α^₄ + β^1x_{14j  +} β^2x_21j = 36.637 + 3.072 - 0.761 x_14j + 0.165x_24j

ŷ_5j = β^₀ + α^₅ + β^1x_{15j  +} β^2x_21j = 36.637 +9.544 - 0.761 x_15j + 0.165x_25j

ŷ6_j = β^₀ + α^₆ + β^1x_{16j  +} β^2x_26j = 36.637 + 5.839 - 0.761 x_16j + 0.165x_26j

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

 6. LSMEANS(Adjusted means) 계산

> lsmeans(model3,~type)

 type    lsmean       SE df  lower.CL  upper.CL

    1  5.674385 1.991087 28  1.595828  9.752942

    2  6.570912 1.946518 28  2.583650 10.558174

    3 13.584122 1.969292 28  9.550210 17.618034

    4  8.746633 1.955440 28  4.741095 12.752171

    5 15.217831 1.987611 28 11.146394 19.289268

    6 11.513284 1.980761 28  7.455879 15.570690

Confidence level used: 0.95

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

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

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

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

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

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

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

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

6. 처리 간 다중비교

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

> print(model3.lsm,omit=1)

$lsmeans

 type    lsmean       SE df  lower.CL  upper.CL

    1  5.674385 1.991087 28  1.595828  9.752942

    2  6.570912 1.946518 28  2.583650 10.558174

    3 13.584122 1.969292 28  9.550210 17.618034

    4  8.746633 1.955440 28  4.741095 12.752171

    5 15.217831 1.987611 28 11.146394 19.289268

    6 11.513284 1.980761 28  7.455879 15.570690

Confidence level used: 0.95

$contrasts

 contrast   estimate       SE df t.ratio p.value

 1 - 2    -0.8965269 2.782316 28  -0.322  0.9995

 1 - 3    -7.9097368 2.768397 28  -2.857  0.0771

 1 - 4    -3.0722481 2.824671 28  -1.088  0.8822

 1 - 5    -9.5434459 2.853395 28  -3.345  0.0257

 1 - 6    -5.8388992 2.839070 28  -2.057  0.3378

 2 - 3    -7.0132099 2.783087 28  -2.520  0.1525

 2 - 4    -2.1757212 2.753630 28  -0.790  0.9669

 2 - 5    -8.6469191 2.802553 28  -3.085  0.0468

 2 - 6    -4.9423723 2.758561 28  -1.792  0.4869

 3 - 4     4.8374887 2.801787 28   1.727  0.5266

 3 - 5    -1.6337091 2.782316 28  -0.587  0.9911

 3 - 6     2.0708376 2.840329 28   0.729  0.9766

 4 - 5    -6.4711978 2.783344 28  -2.325  0.2178

 4 - 6    -2.7666511 2.753890 28  -1.005  0.9125

 5 - 6     3.7045468 2.831753 28   1.308  0.7781

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

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

> names(model3.lsm)

[1] "lsmeans"   "contrasts"

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

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

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

Dunnett vs Tukey

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

다중 분석 시작

> library(multcomp)

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

> summary(tukey)

         Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts

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

Linear Hypotheses:

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

2 - 1 == 0   0.8965     2.7823   0.322   0.9995 

3 - 1 == 0   7.9097     2.7684   2.857   0.0771 .

4 - 1 == 0   3.0722     2.8247   1.088   0.8821 

5 - 1 == 0   9.5434     2.8534   3.345   0.0257 *

6 - 1 == 0   5.8389     2.8391   2.057   0.3374 

3 - 2 == 0   7.0132     2.7831   2.520   0.1524 

4 - 2 == 0   2.1757     2.7536   0.790   0.9669 

5 - 2 == 0   8.6469     2.8026   3.085   0.0469 *

6 - 2 == 0   4.9424     2.7586   1.792   0.4867 

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