이번 포스팅에서는 R Shiny를 사용하여

 

- (1) 두 개의 독립 표본 집단 간 평균 차이 t-검정 (independent two-sample t-test) 

- (2) 신뢰수준 별 신뢰구간 (confidence interval by confidence level)

을 구하는 웹 애플리케이션을 만들어보겠습니다. 

 

예제로 사용한 데이터는 MASS 패키지의 Cars93 데이터프레임이며, 'Origin'  변수의 두 개 집단(USA, Non-USA) 별로 Price, RPM, Length, Width 변수 간의 t-test 를 해보겠습니다. 

 

왼쪽 사이드바 패널에는 t-검정을 할 변수를 선택할 수 있는 콤보 박스와, 신뢰수준(confidence level)을 지정할 수 있는 슬라이드 입력(default = 0.95) 화면을 만들었습니다. 

그리고 메인 패널에는 변수 이름, 신뢰수준, t-검정 결과를 텍스트로 보여줍니다. 

 

R Shiny - independent two-sample t-test

 

library(shiny)

# Define UI for application that does t-test

ui <- fluidPage(

   

   # Application title

   titlePanel("Independent two sample t-test"),

   

   # Select a variable to do t-test by Origin (USA vs. non-USA groups)

   sidebarPanel(

     selectInput("var_x", 

                 label = "Select a Variable to do t-test", 

                 choices = list("Price"="Price", 

                                "RPM"="RPM", 

                                "Length"="Length", 

                                "Width"="Width"), 

                 selected = "Price"), 

     

     sliderInput("conf_level", "Confidence Level:", 

                  min = 0.05, max = 0.99, value = 0.95)

     ),

   

   # Show a t-test result

   mainPanel(

     h5("Variable is:"),

     verbatimTextOutput("var_x_name"), 

     h5("Confidence Level is:"),

     verbatimTextOutput("conf_level"), 

     h4("t-test Result by Origin (USA vs. non-USA independent two samples)"), 

     verbatimTextOutput("t_test")

   )

)

# Define server logic required to do t-test

server <- function(input, output) {

  

  library(MASS)

  cars93_sub <- Cars93[, c('Origin', 'Price', 'RPM', 'Length', 'Width')]

  

  output$var_x_name <- renderText({

    as.character(input$var_x)

  })

  

  output$conf_level <- renderText({

    as.character(input$conf_level)

  })

  

  output$t_test <- renderPrint({

    # independent two-sample t-test

    t.test(cars93_sub[,input$var_x] ~ Origin, 

           data = cars93_sub, 

           alternative = c("two.sided"), 

           var.equal = FALSE, 

           conf.level = input$conf_level)

  })

}

 

# Run the application 

shinyApp(ui = ui, server = server)

 

많은 도움이 되었기를 바랍니다. 

 

Posted by R Friend R_Friend

이번 포스팅에서는 R Shiny를 사용하여 interactive하게 편하고 쉽게 두 연속형 변수 간 관계를 알아볼 수 있도록 

(1) 상관계수 (Correlation Coefficient)를 구하고, 

(2) 산점도 (Scatter Plot)을 그린 후, 

(3) 선형 회귀선 (Linear Regression Line) 을 겹쳐서 그려주는 

웹 애플리케이션을 만들어보겠습니다. 

 

 

R Shiny - Correlation Coefficient, Scatter Plot, Linear Regression Line

 

예제로 사용할 데이터는 MASS 패키지에 내장되어 있는 Cars93 데이터프레임으로서, Price, MPG.city, MPG.highway, EngineSize, Horsepower, RPM 의 6개 연속형 변수를 선별하였습니다. 

 

왼쪽 사이드바 패널에는 산점도의 X축, Y축에 해당하는 연속형 변수를 콤보박스로 선택할 수 있도록 하였으며, 상관계수와 선형회귀모형을 적합할 지 여부를 선택할 수 있는 체크박스도 추가하였습니다. 

 

오른쪽 본문의 메인 패널에는 두 연속형 변수에 대한 산점도 그래프에 선형 회귀선을 겹쳐서 그려서 보여주고, 상관계수와 선형회귀 적합 결과를 텍스트로 볼 수 있도록 화면을 구성하였습니다. 

 

library(shiny)

# Define UI for application that analyze correlation coefficients and regression

ui <- fluidPage(

   

   # Application title

   titlePanel("Correlation Coefficients and Regression"),

   

   # Select 2 Variables, X and Y 

   sidebarPanel(

     selectInput("var_x", 

                 label = "Select a Variable of X", 

                 choices = list("Price"="Price", 

                                "MPG.city"="MPG.city", 

                                "MPG.highway"="MPG.highway", 

                                "EngineSize"="EngineSize", 

                                "Horsepower"="Horsepower", 

                                "RPM"="RPM"), 

                 selected = "RPM"),

     

     selectInput("var_y", 

                 label = "Select a Variable of Y", 

                 choices = list("Price"="Price", 

                                "MPG.city"="MPG.city", 

                                "MPG.highway"="MPG.highway", 

                                "EngineSize"="EngineSize", 

                                "Horsepower"="Horsepower", 

                                "RPM"="RPM"), 

                 selected = "MPG.highway"), 

     

     checkboxInput(inputId = "corr_checked", 

                   label = strong("Correlation Coefficients"), 

                   value = TRUE), 

     

     checkboxInput(inputId = "reg_checked", 

                   label = strong("Simple Regression"), 

                   value = TRUE)

      ),

      # Show a plot of the generated distribution

      mainPanel(

        h4("Scatter Plot"), 

        plotOutput("scatterPlot"), 

        

        h4("Correlation Coefficent"), 

        verbatimTextOutput("corr_coef"), 

        

        h4("Simple Regression"), 

        verbatimTextOutput("reg_fit")

      )

   )

# Define server logic required to analyze correlation coefficients and regression

server <- function(input, output) {

  library(MASS)

  scatter <- Cars93[,c("Price", "MPG.city", "MPG.highway", "EngineSize", "Horsepower", "RPM")]

  

  # scatter plot

  output$scatterPlot <- renderPlot({

    var_name_x <- as.character(input$var_x)

    var_name_y <- as.character(input$var_y)

    

    plot(scatter[, input$var_x],

         scatter[, input$var_y],

         xlab = var_name_x,

         ylab = var_name_y,

         main = "Scatter Plot")

    

    # add linear regression line

    fit <- lm(scatter[, input$var_y] ~ scatter[, input$var_x])

    abline(fit)

    })

  

  # correlation coefficient

  output$corr_coef <- renderText({

    if(input$corr_checked){

      cor(scatter[, input$var_x], 

          scatter[, input$var_y])

      }

    })

  

  # simple regression

  output$reg_fit <- renderPrint({

    if(input$reg_checked){

      fit <- lm(scatter[, input$var_y] ~ scatter[, input$var_x])

      names(fit$coefficients) <- c("Intercept", input$var_x)

      summary(fit)$coefficients

    }

  })

}

# Run the application 

shinyApp(ui = ui, server = server)

 

많은 도움이 되었기를 바랍니다. 

 

Posted by R Friend R_Friend

R Shiny를 이용하면 탐색적 데이터 분석(Exploratory Data Analysis)을 할 때 요약 통계량, 그래프를 변수, 기간, 통계량, 그래프 종류 등을 자유롭게 조건 선택하고 조합하여 편리하게 사용할 수 있습니다. 

 

이번 포스팅에서는 R Shiny로 DataFrame 의

(1) 특정 연속형 변수를 선택하고 

(2) 분석하고 싶은 요약 통계량을 선택했을 때 

(3) 요약 통계량 결과를 interactive하게 볼 수 있는 웹 애플리케이션을 만들어보겠습니다. 

 

Output Image는 아래와 같습니다. 

 

[ R Shiny 를 활용한 연속형 변수 요약 통계량 조회 웹 애플리케이션 화면 ]

R Shiny - Summary Statistics

 

먼저 R Studio 에서 'File > New File > Shiny Web App' 을 선택해주세요. 

 

 

다음으로, (a) 애플리케이션의 이름, (b) 애플리케이션 유형 (간단한 것이면 하나의 파일에 UI와 Server를 함께 만드는게 편하구요, 복잡한 것이면 ui.R, server.R 의 두 개의 파일로 분리하는 것이 편리함), (c) 애플리케이션 저장 경로 (위치)를 지정해주고, Create 해줍니다. 

 

 

이제 R Shiny 코드를 짜볼 텐데요,

먼저 UI 부분에서는 (a) 왼쪽 사이드바에서 분석 대상이 되는 변수를 콤보 박스로 선택하는 화면, (b) 왼쪽 사이드바에서 요약 통계량을 선택하는 체크 박스 화면, (c) 중앙 메인 패널에서 보여줄 제목과 요약 통계량 테이블을 지정해줍니다. 

 

그리고 Server 부분에서는 (a) MASS 패키지의 Cars93 데이터 프레임으로 부터 일부 6개의 연속형 변수 데이터 읽어오기, (b) 요약 통계량 계산하기 (단, 결측값은 제외하기), (c) 요약 통계량들을 하나의 데이터 프레임으로 만들기, (d) 데이터 프레임 테이블을 렌더링 해줍니다. 

 

#------------------------------------------------------------------ 
# Summary Statistics of a Continuous Variable using RShiny 
#------------------------------------------------------------------

library(shiny)

# Define UI for application that calculate summary statistics
ui <- fluidPage(
   
   # Application title
   titlePanel("Descriptive Statistics of Automobiles"),
   
   # select Input from lists 
   sidebarPanel(
     selectInput("var", "Select a Variable to Analyze", 
                 choice = c("Price"=1, 
                            "MPG.city"=2, 
                            "MPG.highway"=3, 
                            "EngineSize"=4, 
                            "Horsepower"=5, 
                            "RPM"=6), 
                 selectize = F), 
     
     checkboxGroupInput("checked", "Check summary statistics", 
                        list(mean = "mean", 
                             std_dev = "std_dev", 
                             median = "median", 
                             min = "min", 
                             max = "max", 
                             obs_num = "obs_num", 
                             quartile_1 = "quartile_1", 
                             quartile_3 = "quartile_3"), 
                        selected = "mean")
   ), 
   mainPanel(
     h4("Summary Statistics of a variable"),
     tableOutput("desc_statistics")
   )
)

# Define server logic required to calculate summary statistics
server <- function(input, output) {
  # read dataset
  library(MASS)
  descriptive <- Cars93[,c("Price", "MPG.city", "MPG.highway", "EngineSize", "Horsepower", "RPM")]

  output$desc_statistics <- renderTable({
    # summary statistics
    variable <- colnames(descriptive)[as.numeric(input$var)]
    mean <- mean(descriptive[, as.numeric(input$var)], na.rm=TRUE)
    std_dev <- sd(descriptive[, as.numeric(input$var)], na.rm=TRUE)
    median <- median(descriptive[, as.numeric(input$var)], na.rm=TRUE)
    min <- min(descriptive[, as.numeric(input$var)], na.rm=TRUE)
    max <- max(descriptive[, as.numeric(input$var)], na.rm=TRUE)
    obs_num <- nrow(descriptive)
    quartile_1 <- quantile(descriptive[, as.numeric(input$var)], 0.25, na.rm=TRUE)
    quartile_3 <- quantile(descriptive[, as.numeric(input$var)], 0.75, na.rm=TRUE)
    
    # make a DataFrame
    desc_summary <- data.frame(variable, 
                               mean, std_dev, 
                               median, min, max, 
                               obs_num, 
                               quartile_1, quartile_3)
    
    # render Table
    desc_summary[, c("variable", input$checked), drop = FALSE]
  })
}

 


# Run the application 
shinyApp(ui = ui, server = server)

 

 

마지막으로, shinyApp(ui = ui, server = server) 로 UI와 Server를 지정해주고, 저장한 후에, RStudio 우측 상단의 'Run App' 단추를 눌러서 실행시켜주면 됩니다. 

 

많은 도움이 되었기를 바랍니다. 

 

Posted by R Friend R_Friend

Greenplum DB에 R이나 Python, Perl, Java 등의 Procedural Language Extention을 설치해서 대용량 데이터를 In-Database 분산 병렬 처리, 분석할 수 있습니다. 

 

이번 포스팅에서는 인터넷이 되는 환경에서 다운로드한 R 패키지들을 회사/ 기관 정책상 폐쇄망으로 운영하는 환경에서 Greenplum DB에 설치하는 방법을 소개하겠습니다. 

 

1. Greenplum PL/R Extention (Procedural Language R) 설치 방법

2. Greenplum DB에 R 패키지 설치 방법

 

PL/R on Greenplum Database

 

1. Greenplum PL/R Extention 설치 방법

PL/R 은 procedural language 로서, PL/R Extension을 설치하면 Greenplum DB에서 R 프로그래밍 언어, R 패키지의 함수와 데이터셋을 사용할 수 있습니다. 

 

Greenplum DB에 PL/R 확장 언어 설치 방법은 https://gpdb.docs.pivotal.io/5180/ref_guide/extensions/pl_r.html 를 참고하였습니다. 웹 페이지의 상단에서 사용 중인 Greenplum DB version을 선택해주세요. (아래 예는 GPDB v5.18 선택 시 화면)

 

PL/R 은 패키지 형태로 되어 있으며, Pivotal Network(https://network.pivotal.io/products/pivotal-gpdb)에서 다운로드 할 수 있고, Greenplum Package Manager (gppkg) 를 사용해서 쉽게 설치할 수 있습니다. 

 

Greenplum Package Manager (gppkg) 유틸리티는 Host와 Cluster에 PL/R 과 의존성있는 패키지들을 한꺼번에 설치를 해줍니다. 또한 gppkg는 시스템 확장이나 세그먼트 복구 시에 자동으로 PL/R extension을 설치해줍니다. 

 

Greenplum PL/R Extention 설치 순서는 아래와 같습니다. 

 

(0) 먼저, Greenplum DB 작동 중이고, source greenplum_path.sh 실행,  $MASTER_DATA_DIRECTORY, $GPHOME variables 설정 완료 필요합니다. 

psql에서 Greenplum DB 버전을 확인합니다. 

psql # sql -c “select version;”

 

master host에서 gpadmin 계정으로 작업 디렉토리를 만듭니다.

(예: /home/gpadmin/packages)

 

(1) Pivotal Network에서 사용 중인 Greenplum DB version에 맞는  PL/R Extension을 다운로드 합니다. 

(예: plr-2.3.3-gp5-rhel7-x86_64.gppkg)

 

(2) 다운로드 한 PL/R Extension Package를  scp 나 sftp 를 이용해서 Greenplum DB master host로 복사합니다. (아마 회사 정책 상 DBA만 root 권한에 접근 가능한 경우가 대부분일 것이므로, 그런 경우에는 DBA에게 복사/설치 요청을 하셔야 합니다). 

$ scp plr-2.3.3-gp5-rhel7-x86_64.gppkg root@mdw:~/packages

 

(3) PL/R Extension Package를 gppkg 커맨드를 실행하여 설치합니다. (아래 예는 Linux에서 실행한 예)

$ gppkg -i plr-2.3.3-gp5-rhel7-x86_64.gppkg

 

(4) Greenplum DB를 재실행 합니다.

(GPDB를 껐다가 켜는 것이므로 DBA에게 반드시 사전 통보, 허락 받고 실행 필요합니다!)

$ gpstop -r

 

(5) Source the file $GPHOME/greenplum_path.sh

# source /usr/local/greenplum-db/greenplum_path.sh

 

R extension과 R 환경은 아래 경로에 설치되어 있습니다. 

$ GPHOME/ext/R-2.3.3/

 

(6) 각 데이터베이스가 PL/R 언어를 사용하기 위해서는 SQL 문으로 CREATE LANGUAGE  또는 createlang 유틸리티로 PL/R을 등록해주어야 합니다. (아래는 testdb 데이터베이스에 등록하는 예)

$ createlang plr -d testdb

이렇게 하면 PL/R이 untrusted language 로 등록이 되었습니다. 

 

 

참고로, Database 확인은 psql 로 \l 해주면 됩니다. 

psql # \l

 

 

 

2. Greenplum DB에 R 패키지 설치 방법 (Installing external R packages)

 

(0) 필요한 R 패키지, 그리고 이에 의존성이 있는 R 패키지를 한꺼번에 다운로드 합니다. (=> https://rfriend.tistory.com/441 참조)

 

(1) 다운로드한 R 패키지들을 압축하여 Greenplum DB 서버로 복사합니다. 

 

다운로드한 R 패키지들 조회해보겠습니다. 

[root@mdw /]# find . | grep sp_1.3-1.tar.gz
./home/gpadmin/r-pkg/sp_1.3-1.tar.gz
[root@mdw /]# exit
logout
[gpadmin@mdw tmp]$ cd ~
[gpadmin@mdw ~]$ cd r-pkg
[gpadmin@mdw r-pkg]$ ls -la
total 47032
drwxrwxr-x 2 gpadmin gpadmin    4096 Apr 23 13:17 .
drwx------ 1 gpadmin gpadmin    4096 Apr 23 13:14 ..
-rw-rw-r-- 1 gpadmin gpadmin  931812 Apr 23 12:55 DBI_1.0.0.tar.gz
-rw-rw-r-- 1 gpadmin gpadmin  794622 Apr 23 12:55 LearnBayes_2.15.1.tar.gz
-rw-rw-r-- 1 gpadmin gpadmin  487225 Apr 23 12:55 MASS_7.3-51.3.tar.gz
-rw-rw-r-- 1 gpadmin gpadmin 1860456 Apr 23 12:55 Matrix_1.2-17.tar.gz
-rw-rw-r-- 1 gpadmin gpadmin   31545 Apr 23 12:55 R6_2.4.0.tar.gz
-rw-rw-r-- 1 gpadmin gpadmin 3661123 Apr 23 12:55 Rcpp_1.0.1.tar.gz
-rw-rw-r-- 1 gpadmin gpadmin   21810 Apr 23 12:55 abind_1.4-5.tar.gz
-rw-rw-r-- 1 gpadmin gpadmin  231855 Apr 23 12:55 boot_1.3-20.tar.gz
-rw-rw-r-- 1 gpadmin gpadmin   17320 Apr 23 12:55 classInt_0.3-1.tar.gz
-rw-rw-r-- 1 gpadmin gpadmin   19757 Apr 23 12:55 class_7.3-15.tar.gz
-rw-rw-r-- 1 gpadmin gpadmin   73530 Apr 23 12:55 coda_0.19-2.tar.gz
-rw-rw-r-- 1 gpadmin gpadmin  658694 Apr 23 12:55 crayon_1.3.4.tar.gz
-rw-rw-r-- 1 gpadmin gpadmin   80772 Apr 23 12:55 deldir_0.1-16.tar.gz
-rw-rw-r-- 1 gpadmin gpadmin  128553 Apr 23 12:55 digest_0.6.18.tar.gz
-rw-rw-r-- 1 gpadmin gpadmin  582415 Apr 23 12:55 e1071_1.7-1.tar.gz
-rw-rw-r-- 1 gpadmin gpadmin  137075 Apr 23 12:55 expm_0.999-4.tar.gz
-rw-rw-r-- 1 gpadmin gpadmin  347295 Apr 23 12:55 foreign_0.8-71.tar.gz
-rw-rw-r-- 1 gpadmin gpadmin 1058430 Apr 23 12:55 gdata_2.18.0.tar.gz
-rw-rw-r-- 1 gpadmin gpadmin  758133 Apr 23 12:55 geosphere_1.5-7.tar.gz
-rw-rw-r-- 1 gpadmin gpadmin   33783 Apr 23 12:55 gmodels_2.18.1.tar.gz
-rw-rw-r-- 1 gpadmin gpadmin   12577 Apr 23 12:55 goftest_1.1-1.tar.gz
-rw-rw-r-- 1 gpadmin gpadmin  187516 Apr 23 12:55 gtools_3.8.1.tar.gz
-rw-rw-r-- 1 gpadmin gpadmin   45408 Apr 23 12:55 htmltools_0.3.6.tar.gz
-rw-rw-r-- 1 gpadmin gpadmin 1758514 Apr 23 12:55 httpuv_1.5.1.tar.gz
-rw-rw-r-- 1 gpadmin gpadmin 1052728 Apr 23 12:55 jsonlite_1.6.tar.gz
-rw-rw-r-- 1 gpadmin gpadmin   40293 Apr 23 12:55 later_0.8.0.tar.gz
-rw-rw-r-- 1 gpadmin gpadmin  359031 Apr 23 12:55 lattice_0.20-38.tar.gz
-rw-rw-r-- 1 gpadmin gpadmin  200504 Apr 23 12:55 magrittr_1.5.tar.gz
-rw-rw-r-- 1 gpadmin gpadmin 1581592 Apr 23 12:55 maptools_0.9-5.tar.gz
-rw-rw-r-- 1 gpadmin gpadmin  915991 Apr 23 12:55 mgcv_1.8-28.tar.gz
-rw-rw-r-- 1 gpadmin gpadmin   12960 Apr 23 12:55 mime_0.6.tar.gz
-rw-rw-r-- 1 gpadmin gpadmin   79619 Apr 23 12:55 polyclip_1.10-0.tar.gz
-rw-rw-r-- 1 gpadmin gpadmin  106866 Apr 23 12:55 promises_1.0.1.tar.gz
-rw-rw-r-- 1 gpadmin gpadmin  255244 Apr 23 12:55 rgeos_0.4-2.tar.gz
-rw-rw-r-- 1 gpadmin gpadmin  858992 Apr 23 12:55 rlang_0.3.4.tar.gz
-rw-rw-r-- 1 gpadmin gpadmin  639286 Apr 23 12:55 rpart_4.1-15.tar.gz
-rw-rw-r-- 1 gpadmin gpadmin 8166770 Apr 23 12:55 sf_0.7-3.tar.gz
-rw-rw-r-- 1 gpadmin gpadmin 2991469 Apr 23 12:55 shiny_1.3.2.tar.gz
-rw-rw-r-- 1 gpadmin gpadmin   24155 Apr 23 12:55 sourcetools_0.1.7.tar.gz
-rw-rw-r-- 1 gpadmin gpadmin 3485268 Apr 23 12:55 spData_0.3.0.tar.gz
-rw-rw-r-- 1 gpadmin gpadmin 1133621 Apr 23 12:55 sp_1.3-1.tar.gz
-rw-rw-r-- 1 gpadmin gpadmin 2861828 Apr 23 12:55 spatstat.data_1.4-0.tar.gz
-rw-rw-r-- 1 gpadmin gpadmin   65106 Apr 23 12:55 spatstat.utils_1.13-0.tar.gz
-rw-rw-r-- 1 gpadmin gpadmin 6598638 Apr 23 12:55 spatstat_1.59-0.tar.gz
-rw-rw-r-- 1 gpadmin gpadmin 1227625 Apr 23 12:55 spdep_1.1-2.tar.gz
-rw-rw-r-- 1 gpadmin gpadmin    2518 Apr 23 12:55 tensor_1.5.tar.gz
-rwxr-xr-x 1 gpadmin gpadmin    2326 Apr 23 13:17 test.sh
-rw-rw-r-- 1 gpadmin gpadmin  917316 Apr 23 12:55 units_0.6-2.tar.gz
-rw-rw-r-- 1 gpadmin gpadmin  564589 Apr 23 12:55 xtable_1.8-4.tar.gz

 

R 패키지들이 들어있는 폴더를 r-pkg.tar 이름으로 압축해보겠습니다. 

[gpadmin@mdw r-pkg]$ pwd
/home/gpadmin/r-pkg
[gpadmin@mdw r-pkg]$ cd ..
[gpadmin@mdw ~]$ tar cf r-pkg.tar r-pkg
[gpadmin@mdw ~]$ ls -lrt
total 47000
drwxr-xr-x 2 gpadmin gpadmin     4096 Aug 13  2018 gpconfigs
drwxr-xr-x 2 root    root        4096 Mar 22 07:02 gppkgs
drwxrwxr-x 1 gpadmin gpadmin     4096 Apr 23 12:48 gpAdminLogs
-rw-rw-r-- 1 gpadmin gpadmin      983 Apr 23 13:14 pkg.r
drwxrwxr-x 2 gpadmin gpadmin     4096 Apr 23 13:17 r-pkg
-rw-rw-r-- 1 gpadmin gpadmin 48107520 Apr 25 01:52 r-pkg.tar

 

명령 프롬프트 창에서 GPDB Docker 에서 압축한 파일을 로커로 복사 후에 ==> 다른 GPDB 서버로 복사하고 압축을 풀어줍니다. (저는 Docker 환경에서 하다보니 좀 복잡해졌는데요, 만약 로컬에서 R 패키지 다운받았으면 로컬에서 바로 GPDB 서버로 복사하면 됩니다. 압축한 R패키지 파일을 scp로 복사하거나 sftp로 업로드할 수 있으며, 권한이 없는 경우 DBA에게 요청하시구요.) 아래는 mdw에서 root 계정으로 시작해서 다운로드해서 압축한 R 패키지 파일을 scp로  /root/packages 경로에 복사하는 스크립트입니다. 

-- GPDB Docker에서 압축한 파일을 로컬로 복사하기
-- 다른 명령 프롬프트 창에서 복사해오고 확인하기

ihongdon-ui-MacBook-Pro:Downloads ihongdon$ docker cp gpdb-ds:/home/gpadmin/r-pkg.tar /Users/ihongdon/Downloads/r-pkg.tar
ihongdon-ui-MacBook-Pro:Downloads ihongdon$
ihongdon-ui-MacBook-Pro:Downloads ihongdon$
ihongdon-ui-MacBook-Pro:Downloads ihongdon$ ls -lrt
-rw-rw-r--   1 ihongdon  staff  48107520  4 25 10:52 r-pkg.tar

-- 다른 GPDB 서버로 복사하기
ihongdon-ui-MacBook-Pro:Downloads ihongdon$ scp r-pkg.tar root@mdw:~/package

-- 압축 해제
$ tar -xvf r-pkg.tar

 

Greenplum DB에 R 패키지를 설치하려면 모든 Greenplum 서버에 R이 이미 설치되어 있어야 합니다. 

여러개의 Segments 에 동시에 R 패키지들을 설치해주기 위해서 배포하고자 하는 host list를 작성해줍니다. 

# source /usr/local/greenplum-db/greenplum_path.sh
# vi hostfile_packages

 

vi editor 창이 열리면 아래처럼 R을 설치하고자 하는 host 이름을 등록해줍니다. (1개 master, 3개 segments 예시)

-- vi 편집창에서 --
smdw
sdw1
sdw2
sdw3
~
~
~
esc 누르고 :wq!

 

명령 프롬프트 창에서 mdw로 부터 root 계정으로 각 노드에 package directory 를 복사해줍니다. 

# gpscp -f hostfile_packages -r packages =:/root

 

hostfile_packages를 복사해서 hostfile_all 을 만들고, mdw를 추가해줍니다. 

-- copy
$ cp hostfile_packages  hostfile_all

-- insert mdw
$ vi hostfile_all
-- vi 편집창에서 --
mdw
smdw
sdw1
sdw2
sdw3
~
~
~
esc 누르고 :wq!

 

mdw를 포함한 모든 서버에 R packages 를 설치하는 'R CMD INSTALL r_package_name' 명령문을 mdw에서 실행합니다. (hostfile_all 에 mdw, smdw, sdw1, sdw2, sdw3 등록해놓았으므로 R이 모든 host에 설치됨)

$ pssh -f hostfile_all -v -e 'R CMD INSTALL ./DBI_1.0.0.tar.gz 
LearnBayes_2.15.1.tar.gz MASS_7.3-51.3.tar.gz Matrix_1.2-17.tar.gz 
R6_2.4.0.tar.gz Rcpp_1.0.1.tar.gz 
abind_1.4-5.tar.gz boot_1.3-20.tar.gz classInt_0.3-1.tar.gz
class_7.3-15.tar.gz coda_0.19-2.tar.gz crayon_1.3.4.tar.gz
deldir_0.1-16.tar.gz digest_0.6.18.tar.gz e1071_1.7-1.tar.gz
expm_0.999-4.tar.gz foreign_0.8-71.tar.gz gdata_2.18.0.tar.gz
geosphere_1.5-7.tar.gz gmodels_2.18.1.tar.gz goftest_1.1-1.tar.gz
gtools_3.8.1.tar.gz htmltools_0.3.6.tar.gz httpuv_1.5.1.tar.gz
jsonlite_1.6.tar.gz later_0.8.0.tar.gz lattice_0.20-38.tar.gz
magrittr_1.5.tar.gz maptools_0.9-5.tar.gz mgcv_1.8-28.tar.gz
mime_0.6.tar.gz polyclip_1.10-0.tar.gz promises_1.0.1.tar.gz
rgeos_0.4-2.tar.gz rlang_0.3.4.tar.gz rpart_4.1-15.tar.gz
sf_0.7-3.tar.gz shiny_1.3.2.tar.gz sourcetools_0.1.7.tar.gz
spData_0.3.0.tar.gz sp_1.3-1.tar.gz spatstat.data_1.4-0.tar.gz
spatstat.utils_1.13-0.tar.gz spatstat_1.59-0.tar.gz spdep_1.1-2.tar.gz
tensor_1.5.tar.gz units_0.6-2.tar.gz xtable_1.8-4.tar.gz'

 

특정 R 패키지를 설치하려고 할 때, 만약 의존성 있는 패키지 (dependencies packages) 가 이미 설치되어 있지 않다면 특정 R 패키지는 설치가 되지 않습니다. 따라서 위의 'R CMD INSTALL r-package-names' 명령문을 실행하면 설치가 되는게 있고, 안되는 것(<- 의존성 있는 패키지가 먼저 설치된 이후에나 설치 가능)도 있게 됩니다. 따라서 이 설치 작업을 수작업으로 반복해서 여러번 돌려줘야 합니다. loop 돌리다보면 의존성 있는 패키지가 설치가 먼저 설치가 될거고, 그 다음에 이전에는 설치가 안되었던게 의존성 있는 패키지가 바로 전에 설치가 되었으므로 이제는 설치가 되고, ...., ....., 다 설치 될때까지 몇 번 더 실행해 줍니다. 

 

많은 도움이 되었기를 바랍니다. 

Posted by R Friend R_Friend

회사의 보안 정책 상 사내 폐쇄망으로 운영함에 따라 R 패키지를 인터넷으로부터 바로 다운로드하지 못하는 경우가 있습니다. 이처럼 폐쇄망일 경우 R 패키지를 설치하는 것이 '지옥 그 자체!' 일 수 있습니다. 왜냐하면 특정 R 패키지를 설치하려면 그와 의존성이 있는 다른 패키지를 설치해야만 하고, 그 의존성이 있는 패키지는 또 다른 의존성이 있는 패키지가 설치되어야만 하고.... 꼬리에 꼬리를 무는 의존성 있는 패키지들을 설치하다보면 입에서 투덜거림이 궁시렁 궁시렁 나오게 됩니다. -_-;;; 

 

이럴 때 편리하게 사용할 수 있는 'R 패키지 다운로드 & 의존성 있는 R 패키지도 같이 한꺼번에 다운로드 하기' 하는 방법을 소개하겠습니다. 

 

 

먼저 R 패키지를 다운로드 받아 놓을 'r-pkt' 폴더를 만들어보겠습니다. 

> mainDir <- "/Users/ihongdon/Downloads"

> subDir <- "r-pkg"

dir.create(file.path(mainDir, subDir), showWarnings = FALSE)

 

다음으로 tools 패키지의 package_dependencies() 함수를 이용하여 다운도르하려는 특정 R 패키지와 이 패키지의 의존성 있는 패키지(dependencies of package)들 정보를 가져오는 사용자 정의 함수 getPackages() 를 만들어보겠습니다. 

# UDF for get packages with dependencies

getPackages <- function(packs){

  packages <- unlist(

    # Find (recursively) dependencies or reverse dependencies of packages.

    tools::package_dependencies(packs, available.packages(),

                         which=c("Depends", "Imports"), recursive=TRUE)

  )

  packages <- union(packs, packages)

  

  return(packages)

}

 

 

자, 이제 준비가 되었으니 예제로 "dplyr"와 "ggplot2" 의 두 개 패키지와 이들과 의존성을 가지는 패키지들을 getPackages() 사용자 정의 함수로 정보를 가져온 후에, download.packages() 함수로 '/Users/ihongdon/Downloads/r-pkg' 폴더로 다운로드 해보도록 하겠습니다. 

> packages <- getPackages(c("dplyr", "ggplot2"))

download.packages(packages, destdir=file.path(mainDir, subDir))

 

trying URL 'https://cran.rstudio.com/src/contrib/dplyr_0.8.0.1.tar.gz'

Content type 'application/x-gzip' length 1075146 bytes (1.0 MB)

==================================================

downloaded 1.0 MB

trying URL 'https://cran.rstudio.com/src/contrib/ggplot2_3.1.1.tar.gz'

Content type 'application/x-gzip' length 2862022 bytes (2.7 MB)

==================================================

downloaded 2.7 MB

trying URL 'https://cran.rstudio.com/src/contrib/assertthat_0.2.1.tar.gz'

Content type 'application/x-gzip' length 12742 bytes (12 KB)

==================================================

downloaded 12 KB

trying URL 'https://cran.rstudio.com/src/contrib/glue_1.3.1.tar.gz'

Content type 'application/x-gzip' length 122950 bytes (120 KB)

==================================================

downloaded 120 KB

trying URL 'https://cran.rstudio.com/src/contrib/magrittr_1.5.tar.gz'

Content type 'application/x-gzip' length 200504 bytes (195 KB)

==================================================

downloaded 195 KB

Warning in download.packages(packages, destdir = file.path(mainDir, subDir)) :

  no package 'methods' at the repositories

trying URL 'https://cran.rstudio.com/src/contrib/pkgconfig_2.0.2.tar.gz'

Content type 'application/x-gzip' length 6024 bytes

==================================================

downloaded 6024 bytes

trying URL 'https://cran.rstudio.com/src/contrib/R6_2.4.0.tar.gz'

Content type 'application/x-gzip' length 31545 bytes (30 KB)

==================================================

downloaded 30 KB

trying URL 'https://cran.rstudio.com/src/contrib/Rcpp_1.0.1.tar.gz'

Content type 'application/x-gzip' length 3661123 bytes (3.5 MB)

==================================================

downloaded 3.5 MB

trying URL 'https://cran.rstudio.com/src/contrib/rlang_0.3.4.tar.gz'

Content type 'application/x-gzip' length 858992 bytes (838 KB)

==================================================

downloaded 838 KB

trying URL 'https://cran.rstudio.com/src/contrib/tibble_2.1.1.tar.gz'

Content type 'application/x-gzip' length 311836 bytes (304 KB)

==================================================

downloaded 304 KB

trying URL 'https://cran.rstudio.com/src/contrib/tidyselect_0.2.5.tar.gz'

Content type 'application/x-gzip' length 21883 bytes (21 KB)

==================================================

downloaded 21 KB

Warning in download.packages(packages, destdir = file.path(mainDir, subDir)) :

  no package 'utils' at the repositories

Warning in download.packages(packages, destdir = file.path(mainDir, subDir)) :

  no package 'tools' at the repositories

trying URL 'https://cran.rstudio.com/src/contrib/cli_1.1.0.tar.gz'

Content type 'application/x-gzip' length 40232 bytes (39 KB)

==================================================

downloaded 39 KB

trying URL 'https://cran.rstudio.com/src/contrib/crayon_1.3.4.tar.gz'

Content type 'application/x-gzip' length 658694 bytes (643 KB)

==================================================

downloaded 643 KB

trying URL 'https://cran.rstudio.com/src/contrib/fansi_0.4.0.tar.gz'

Content type 'application/x-gzip' length 266123 bytes (259 KB)

==================================================

downloaded 259 KB

trying URL 'https://cran.rstudio.com/src/contrib/pillar_1.3.1.tar.gz'

Content type 'application/x-gzip' length 103972 bytes (101 KB)

==================================================

downloaded 101 KB

trying URL 'https://cran.rstudio.com/src/contrib/purrr_0.3.2.tar.gz'

Content type 'application/x-gzip' length 373701 bytes (364 KB)

==================================================

downloaded 364 KB

Warning in download.packages(packages, destdir = file.path(mainDir, subDir)) :

  no package 'grDevices' at the repositories

trying URL 'https://cran.rstudio.com/src/contrib/utf8_1.1.4.tar.gz'

Content type 'application/x-gzip' length 218882 bytes (213 KB)

==================================================

downloaded 213 KB

trying URL 'https://cran.rstudio.com/src/contrib/digest_0.6.18.tar.gz'

Content type 'application/x-gzip' length 128553 bytes (125 KB)

==================================================

downloaded 125 KB

Warning in download.packages(packages, destdir = file.path(mainDir, subDir)) :

  no package 'grid' at the repositories

trying URL 'https://cran.rstudio.com/src/contrib/gtable_0.3.0.tar.gz'

Content type 'application/x-gzip' length 368081 bytes (359 KB)

==================================================

downloaded 359 KB

trying URL 'https://cran.rstudio.com/src/contrib/lazyeval_0.2.2.tar.gz'

Content type 'application/x-gzip' length 83482 bytes (81 KB)

==================================================

downloaded 81 KB

trying URL 'https://cran.rstudio.com/src/contrib/MASS_7.3-51.4.tar.gz'

Content type 'application/x-gzip' length 487233 bytes (475 KB)

==================================================

downloaded 475 KB

trying URL 'https://cran.rstudio.com/src/contrib/mgcv_1.8-28.tar.gz'

Content type 'application/x-gzip' length 915991 bytes (894 KB)

==================================================

downloaded 894 KB

trying URL 'https://cran.rstudio.com/src/contrib/plyr_1.8.4.tar.gz'

Content type 'application/x-gzip' length 392451 bytes (383 KB)

==================================================

downloaded 383 KB

trying URL 'https://cran.rstudio.com/src/contrib/reshape2_1.4.3.tar.gz'

Content type 'application/x-gzip' length 36405 bytes (35 KB)

==================================================

downloaded 35 KB

trying URL 'https://cran.rstudio.com/src/contrib/scales_1.0.0.tar.gz'

Content type 'application/x-gzip' length 299262 bytes (292 KB)

==================================================

downloaded 292 KB

Warning in download.packages(packages, destdir = file.path(mainDir, subDir)) :

  no package 'stats' at the repositories

trying URL 'https://cran.rstudio.com/src/contrib/viridisLite_0.3.0.tar.gz'

Content type 'application/x-gzip' length 44019 bytes (42 KB)

==================================================

downloaded 42 KB

trying URL 'https://cran.rstudio.com/src/contrib/withr_2.1.2.tar.gz'

Content type 'application/x-gzip' length 53578 bytes (52 KB)

==================================================

downloaded 52 KB

Warning in download.packages(packages, destdir = file.path(mainDir, subDir)) :

  no package 'graphics' at the repositories

trying URL 'https://cran.rstudio.com/src/contrib/nlme_3.1-139.tar.gz'

Content type 'application/x-gzip' length 793473 bytes (774 KB)

==================================================

downloaded 774 KB

trying URL 'https://cran.rstudio.com/src/contrib/Matrix_1.2-17.tar.gz'

Content type 'application/x-gzip' length 1860456 bytes (1.8 MB)

==================================================

downloaded 1.8 MB

Warning in download.packages(packages, destdir = file.path(mainDir, subDir)) :

  no package 'splines' at the repositories

trying URL 'https://cran.rstudio.com/src/contrib/stringr_1.4.0.tar.gz'

Content type 'application/x-gzip' length 135777 bytes (132 KB)

==================================================

downloaded 132 KB

trying URL 'https://cran.rstudio.com/src/contrib/labeling_0.3.tar.gz'

Content type 'application/x-gzip' length 10722 bytes (10 KB)

==================================================

downloaded 10 KB

trying URL 'https://cran.rstudio.com/src/contrib/munsell_0.5.0.tar.gz'

Content type 'application/x-gzip' length 182653 bytes (178 KB)

==================================================

downloaded 178 KB

trying URL 'https://cran.rstudio.com/src/contrib/RColorBrewer_1.1-2.tar.gz'

Content type 'application/x-gzip' length 11532 bytes (11 KB)

==================================================

downloaded 11 KB

trying URL 'https://cran.rstudio.com/src/contrib/lattice_0.20-38.tar.gz'

Content type 'application/x-gzip' length 359031 bytes (350 KB)

==================================================

downloaded 350 KB

trying URL 'https://cran.rstudio.com/src/contrib/colorspace_1.4-1.tar.gz'

Content type 'application/x-gzip' length 2152594 bytes (2.1 MB)

==================================================

downloaded 2.1 MB

trying URL 'https://cran.rstudio.com/src/contrib/stringi_1.4.3.tar.gz'

Content type 'application/x-gzip' length 7290890 bytes (7.0 MB)

==================================================

downloaded 7.0 MB

      [,1]           [,2]                                                       

 [1,] "dplyr"        "/Users/ihongdon/Downloads/r-pkg/dplyr_0.8.0.1.tar.gz"     

 [2,] "ggplot2"      "/Users/ihongdon/Downloads/r-pkg/ggplot2_3.1.1.tar.gz"     

 [3,] "assertthat"   "/Users/ihongdon/Downloads/r-pkg/assertthat_0.2.1.tar.gz"  

 [4,] "glue"         "/Users/ihongdon/Downloads/r-pkg/glue_1.3.1.tar.gz"        

 [5,] "magrittr"     "/Users/ihongdon/Downloads/r-pkg/magrittr_1.5.tar.gz"      

 [6,] "pkgconfig"    "/Users/ihongdon/Downloads/r-pkg/pkgconfig_2.0.2.tar.gz"   

 [7,] "R6"           "/Users/ihongdon/Downloads/r-pkg/R6_2.4.0.tar.gz"          

 [8,] "Rcpp"         "/Users/ihongdon/Downloads/r-pkg/Rcpp_1.0.1.tar.gz"        

 [9,] "rlang"        "/Users/ihongdon/Downloads/r-pkg/rlang_0.3.4.tar.gz"       

[10,] "tibble"       "/Users/ihongdon/Downloads/r-pkg/tibble_2.1.1.tar.gz"      

[11,] "tidyselect"   "/Users/ihongdon/Downloads/r-pkg/tidyselect_0.2.5.tar.gz"  

[12,] "cli"          "/Users/ihongdon/Downloads/r-pkg/cli_1.1.0.tar.gz"         

[13,] "crayon"       "/Users/ihongdon/Downloads/r-pkg/crayon_1.3.4.tar.gz"      

[14,] "fansi"        "/Users/ihongdon/Downloads/r-pkg/fansi_0.4.0.tar.gz"       

[15,] "pillar"       "/Users/ihongdon/Downloads/r-pkg/pillar_1.3.1.tar.gz"      

[16,] "purrr"        "/Users/ihongdon/Downloads/r-pkg/purrr_0.3.2.tar.gz"       

[17,] "utf8"         "/Users/ihongdon/Downloads/r-pkg/utf8_1.1.4.tar.gz"        

[18,] "digest"       "/Users/ihongdon/Downloads/r-pkg/digest_0.6.18.tar.gz"     

[19,] "gtable"       "/Users/ihongdon/Downloads/r-pkg/gtable_0.3.0.tar.gz"      

[20,] "lazyeval"     "/Users/ihongdon/Downloads/r-pkg/lazyeval_0.2.2.tar.gz"    

[21,] "MASS"         "/Users/ihongdon/Downloads/r-pkg/MASS_7.3-51.4.tar.gz"     

[22,] "mgcv"         "/Users/ihongdon/Downloads/r-pkg/mgcv_1.8-28.tar.gz"       

[23,] "plyr"         "/Users/ihongdon/Downloads/r-pkg/plyr_1.8.4.tar.gz"        

[24,] "reshape2"     "/Users/ihongdon/Downloads/r-pkg/reshape2_1.4.3.tar.gz"    

[25,] "scales"       "/Users/ihongdon/Downloads/r-pkg/scales_1.0.0.tar.gz"      

[26,] "viridisLite"  "/Users/ihongdon/Downloads/r-pkg/viridisLite_0.3.0.tar.gz

[27,] "withr"        "/Users/ihongdon/Downloads/r-pkg/withr_2.1.2.tar.gz"       

[28,] "nlme"         "/Users/ihongdon/Downloads/r-pkg/nlme_3.1-139.tar.gz"      

[29,] "Matrix"       "/Users/ihongdon/Downloads/r-pkg/Matrix_1.2-17.tar.gz"     

[30,] "stringr"      "/Users/ihongdon/Downloads/r-pkg/stringr_1.4.0.tar.gz"     

[31,] "labeling"     "/Users/ihongdon/Downloads/r-pkg/labeling_0.3.tar.gz"      

[32,] "munsell"      "/Users/ihongdon/Downloads/r-pkg/munsell_0.5.0.tar.gz"     

[33,] "RColorBrewer" "/Users/ihongdon/Downloads/r-pkg/RColorBrewer_1.1-2.tar.gz"

[34,] "lattice"      "/Users/ihongdon/Downloads/r-pkg/lattice_0.20-38.tar.gz"   

[35,] "colorspace"   "/Users/ihongdon/Downloads/r-pkg/colorspace_1.4-1.tar.gz"  

[36,] "stringi"      "/Users/ihongdon/Downloads/r-pkg/stringi_1.4.3.tar.gz"

 

위에 다운로드된 패키지들의 리스트를 보면 알 수 있는 것처럼, 'dplyr'과 'ggplot2'의 두 개 패키지를 다운로드 하려고 했더니 이들이 의존성을 가지고 있는 [3] 'assertthat' 패키지부터 ~ [36] 'stringi' 패키지까지 34개의 패키지가 추가로 다운로드 되었습니다. 

 

외부의 인터넷 연결이 되는 환경에서 다운받아 놓은 이들 R packages 파일들을 저장매체에 저장하여 가져가서 폐쇄망으로 되어 있는 서버에 복사를 한 후에, 수작업으로 R package 설치하면 되겠습니다. (이때 의존성 있는 패키지들까지 알아서 다 설치가 되도록 여러번 설치 작업을 반복해주면 됩니다. 

 

##==============================================##

폐쇄망 환경에서 Greenplum DB에 R 패키지 설치하는 방법은 

==> https://rfriend.tistory.com/442 

포스팅을 참고하세요. 

##==============================================##

 

많은 도움이 되었기를 바랍니다. 

Posted by R Friend R_Friend

이번 포스팅에서는 행 기준으로 숫자형 변수들의 합을 구한 다음에, 그룹별로 행 기준 합이 최대인 전체 행을 선별하는 방법을 소개하겠습니다. 


말로 설명한 내용만으로는 얼른 이해가 안 올 수도 있겠는데요, 이번에 해볼 내용은 아래의 이미지를 참고하시면 이해가 쉬울 듯 합니다. 




예제로 사용할 하나의 그룹 변수(V1)와 나머지 9개의 숫자형 변수(V2~V10)로 구성된 간단한 DataFrame을 만들어보겠습니다. 



> ##------------------------------------------------------

> ## selecting distinct object using dplyr chain operator

> ##------------------------------------------------------

> rm(list=ls())

> set.seed(123) # for reproducibility

> V1 <- c(rep("apple", 5), rep("banana", 5), rep("tomato", 5)) # group

> V2 <- sample(x=1:10, size=15, replace=T)

> V3 <- sample(x=1:10, size=15, replace=T)

> V4 <- sample(x=1:10, size=15, replace=T)

> V5 <- sample(x=1:10, size=15, replace=T)

> V6 <- sample(x=1:10, size=15, replace=T)

> V7 <- sample(x=1:10, size=15, replace=T)

> V8 <- sample(x=1:10, size=15, replace=T)

> V9 <- sample(x=1:10, size=15, replace=T)

> V10 <- sample(x=1:10, size=15, replace=T)

> df <- data.frame(V1, V2, V3, V4, V5, V6, V7, V8, V9, V10)

> df

       V1 V2 V3 V4 V5 V6 V7 V8 V9 V10

1   apple  3  9 10  2  7  3  2  9   7

2   apple  8  3 10  3  1  4  7 10   4

3   apple  5  1  7  5  4  7  4  7   4

4   apple  9  4  8  3  3  4  7  5   3

5   apple 10 10  1  9  9  2  4  2   4

6  banana  1  9  5  1  5  3  2 10  10

7  banana  6  7  8  5  9  7  8  4   2

8  banana  9  7  3  8  9  5  1  1   1

9  banana  6 10  4  2  8  8  5 10   2

10 banana  5  7  3  6  5  2  6  8   7

11 tomato 10  8  2  3  8  5  6  2   7

12 tomato  5  6  5  2  7 10  4  6   9

13 tomato  7  6  5  8  8  9  5 10   7

14 tomato  6  3  4  9  1  9 10  6   8

15 tomato  2  2  2  4  5  2  5  5   6

 

> rm(V1, V2, V3, V4, V5, V6, V7, V8, V9, V10)





rowSums() 함수를 사용하여서 행(row) 기준의 숫자형 변수들 모두에 대해 합계를 구하여 'V_sum' 이라는 새로운 변수를 추가해보겠습니다. 



> # summation in a row direction

> df$V_sum <- rowSums(df[,2:10])

> df

       V1 V2 V3 V4 V5 V6 V7 V8 V9 V10 V_sum


1   apple  3  9 10  2  7  3  2  9   7    52

2   apple  8  3 10  3  1  4  7 10   4    50

3   apple  5  1  7  5  4  7  4  7   4    44

4   apple  9  4  8  3  3  4  7  5   3    46

5   apple 10 10  1  9  9  2  4  2   4    51


6  banana  1  9  5  1  5  3  2 10  10    46

7  banana  6  7  8  5  9  7  8  4   2    56

8  banana  9  7  3  8  9  5  1  1   1    44

9  banana  6 10  4  2  8  8  5 10   2    55

10 banana  5  7  3  6  5  2  6  8   7    49


11 tomato 10  8  2  3  8  5  6  2   7    51

12 tomato  5  6  5  2  7 10  4  6   9    54

13 tomato  7  6  5  8  8  9  5 10   7    65

14 tomato  6  3  4  9  1  9 10  6   8    56

15 tomato  2  2  2  4  5  2  5  5   6    33

 




위의 행 기준 합계로 보면 'apple' 그룹에서는 1번째 행의 합이 52로 가장 크며, 'banana' 그룹에서는 7번째 행의 합이 56으로서 가장 크고, 'tomato' 그룹에서는 13번째 행의 합이 65로서 가장 큽니다. 이를 1번째, 7번째, 13번째 전체 행을 선별해보겠습니다. (빨간색으로 표시한 부분)



> library(dplyr)

> df_group_distinct_max <- df %>% 

+   arrange(V1, desc(V_sum)) %>% 

+   group_by(V1) %>% slice(1:1)

> df_group_distinct_max <- data.frame(df_group_distinct_max[1:10])

> df_group_distinct_max

      V1 V2 V3 V4 V5 V6 V7 V8 V9 V10

1  apple  3  9 10  2  7  3  2  9   7

2 banana  6  7  8  5  9  7  8  4   2

3 tomato  7  6  5  8  8  9  5 10   7

 


* dplyr 패키지 사용법은 

https://rfriend.tistory.com/234 , 

https://rfriend.tistory.com/236

참고하시기 바랍니다


많은 도움이 되었기를 바랍니다. 


이번 포스팅이 도움이 되었다면 아래의 '공감~'를 꾹 눌러주세요. ^^



Posted by R Friend R_Friend

일반적으로 R ggplot2로 히스토그램을 그릴 때 동일한 폭의 bin 넓이(bin width)를 설정해줍니다. 


이번 포스팅에서는 R ggplot2로 bin의 넓이를 다르게 설정한 히스토그램 그리는 방법을 소개하겠습니다. (histogram with different bin width using R ggplot2)


간단한 예제 데이터프레임을 만들어서 예를 들어보았습니다. 요점은 사용자가 정의한 binwidth에 해당하는 group 변수를 하나 만들구요, geom_histogram() 함수로 히스토그램을 그릴 때 group 별로 subset을 하고 breaks argument로 사용자가 정의한 binwidth 구간을 설정해주는 것입니다. 



#----------------

# histogram with variable size of bin width and different colors per bins using ggplot2

#----------------


# sample data frame

mydf <- data.frame(var = c(1100, 10000, 100000, 190000, 110000, 220000, 550000, 701000, 790000))


# numeric notation for large numbers

options(scipen = 30)


library("ggplot2")


# fill color with different colors per bins

mydf$group <- ifelse(mydf$var < 10000, 1, 

                          ifelse(mydf$var < 100000, 2, 

                                 ifelse(mydf$var < 200000, 3, 

                                        ifelse(mydf$var < 500000, 4, 5))))


# breaks of bin

bins <- c(1000, 10000, 100000, 200000, 500000, 800000)


# draw histogram with variable size of bin width and different colors per bins

ggplot(mydf, aes(x= var)) +

  geom_histogram(data=subset(mydf, group==1), breaks = c(1000, 10000), fill="black") +

  geom_histogram(data=subset(mydf, group==2), breaks = c(10000, 100000), fill="yellow") +

  geom_histogram(data=subset(mydf, group==3), breaks = c(100000, 200000), fill="green") +

  geom_histogram(data=subset(mydf, group==4), breaks = c(200000, 500000), fill="blue") +

  geom_histogram(data=subset(mydf, group==5), breaks = c(500000, 800000), fill="red") +

  scale_x_continuous(breaks = bins, limits = c(1000, 800000)) +

  xlab("variable 1") + 

  ylab("count") +

  ggtitle("Histogram with different size of bin width and colors") + 

  theme(plot.title = element_text(hjust = 0.5, size = 14))


 



많은 도움이 되었기를 바랍니다. 


이번 포스팅이 도움이 되었다면 아래의 '공감~'를 꾹 눌러주세요. :-)



Posted by R Friend R_Friend

이번 포스팅에서는 R ggplot2 패키지의 함수를 Python에서 그대로 사용하여 그래프를 그릴 수 있게 해주는 Python의 PlotNine library를 소개하겠습니다. 


R ggplot2 패키지의 그래프 기능이 굉장히 강력하고 체계적이며 이쁘기 때문에 많은 분석가들이 시각화 도구로 R ggplot2를 사용하고 있습니다. 저도 R ggplot2와 RShiny 의 굉장한 팬이기도 합니다. R ggplot2에 이미 익숙한 사용자라면 R ggplot2 함수를 Python에서 그대로 사용하여 시각화를 할 수 있다면 매우 편리하고 빠르게 시각화를 할 수 있을 것입니다. 



연도별 월별 항공기 탑승객 수를 저장해놓은 flights 데이터프레임을 사용하여 파이썬의 PlotNine 라이브러리를 사용하여 파이썬에서 R ggplot2 함수로 선 그래프와 히트맵을 그려보겠습니다. 



import matplotlib.pyplot as plt

import seaborn as sns

plt.rcParams['figure.figsize'] = [10, 8]

 



먼저 flights 데이터프레임을 seaborn 라이브러리에서 가져와 로딩하였습니다. 



# data loading

flights = sns.load_dataset('flights')

flights.head()

yearmonthpassengers
01949January112
11949February118
21949March132
31949April129
41949May121

 




PlotNine 라이브러리를 처음 사용한다면 아래처럼 명령 프롬프트에서 pip install plotnine 으로 설치해주세요. 


  명령 프롬프트에서 'PlotNine' 라이브러리 설치하기



Microsoft Windows [Version 10.0.17134.523]

(c) 2018 Microsoft Corporation. All rights reserved.


C:\Users\admin>conda env list

# conda environments:

#

base                  *  C:\Users\admin\Anaconda3

py_v35                   C:\Users\admin\Anaconda3\envs\py_v35

py_v36                   C:\Users\admin\Anaconda3\envs\py_v36



C:\Users\admin>activate py_v36


(py_v36) C:\Users\admin>pip install plotnine

 





  (1) PlotNine 라이브러리를 사용하여 파이썬에서 R ggplot2 함수로 선 그래프 그리기


일단 plotnine에서 import * 로 모든 클래스를 불러오고 나면 그 다음부터는 R ggplot2 함수를 정확히 똑같이 사용해서 jupyter notebook에서 파이썬 데이터프레임으로 선 그래프를 그릴 수 있습니다. 신기하고 놀랍지요?!  PlotNine 라이브러리 만들어주신 분들에게 깊은 감사를 드립니다. ^^



# line graph using R ggplot2 function in Python environment by plotnine library

from plotnine import *


fig = plt.figure()


ggplot(flights, aes(x='year', y='passengers', group='month', color='month')) \

    + geom_line() \

    + ggtitle('Time Series Graph of Flight')

 






  (2) PlotNine 라이브러리를 사용하여 파이썬에서 R ggplot2 함수로 히트맵 그리기



# heatmap using R ggplot2 function in Python environment by plotnine library

from plotnine import *


fig = plt.figure()


ggplot(flights, aes(x='year', y='month')) \

    + geom_tile(aes(fill='passengers')) \

    + scale_fill_gradientn(colors=['#9ebcda','#8c6bb1','#88419d','#6e016b']) \

    + ggtitle('Heatmap of Flights by plotnine')



많은 도움이 되었기를 바랍니다. 


이번 포스팅이 도움이 되었다면 아래의 '공감~'를 꾹 눌러주세요.  :-)



Posted by R Friend R_Friend

이번 포스팅에서는 마지막 네번째 포스팅으로서, 12개의 설명변수를 사용하여 유방암 악성, 양성 여부 진단하는 로지스틱 회귀모형을 훈련 데이터셋으로 적합시켜보고, 테스트셋에 대해 모델 성능을 평가하며, 모델의 회귀계수를 가지고 해석을 해보겠습니다. 


[로지스틱 회귀분석을 통한 유방암 예측(분류) 포스팅 순서]

  1. WDBC(Wisconsin Diagnostic Breast Cancer) dataset 소개 및 분석 목적과 방향 설정
  2. 탐색적 데이터 분석 및 전처리
  3. 1차 변수 선택 및 목표변수와 설명변수 간 관계 분석
  4. 로지스틱 회귀모형 적합 및 모델 평가, 해석



이번 포스팅은 아래의 순서대로 분석을 진행해보겠습니다. 

(탐색적 데이터 분석 및 전처리는 이전 포스팅 참고 => http://rfriend.tistory.com/400)


(1) training set (80%) vs. test set (20%) 분리

(2) 로지스틱 회귀모형 적합 (w/ training set)

(3) 후진 소거법 (backward elimination)

(4) 모델 성능 평가 (w/ test set) : confusion matrix, ROC curve

(5) 모델 해석



  (1) training set (80%) vs. test set (20%) 분리


569개 관측치를 가진 데이터셋을 (a) 모형 적합에 사용할 훈련용 데이터셋(training set), (b) 모형 성능 평가에 사용할 테스트 데이터셋(test set) 으로 나누어 보겠습니다. 이번 포스팅에서는 training : test 를 8:2 비율로 나누어서 분석하였습니다. 


만약 training set으로 훈련한 모델의 정확도는 높은데 test set을 대상으로 한 정확도는 낮다면 과적합(overfitting)을 의심해볼 수 있습니다. 




로지스틱 회귀모형을 적합하는데는 나중에 회귀계수 해석하기에 용이하도록 표준화하기 이전의 원래의 데이터를 사용하였으며, training : test = 8:2 로 나누기 위해 base패키지의 sample() 함수를 사용하여 indexing을 위한 무작위 index를 생성해두었습니다. 


참고로, 아래 코드의 train, test, Y.test에는 관측치 측정치가 들어있는 것이 아니라 나중에 indexing을 위한 index 정보 (랜덤 샘플링 indexing 할 위치 정보)가 들어 있습니다. 다음번 glm() 함수에서 subset argument에 train, test index를 사용할 예정입니다. 



> ##=========================

> ## fitting logistic regression

> #----------

> # data set with 12 input variable and Y variable

> wdbc_12 <- data.frame(Y, wdbc[,x_names_sorted])

> str(wdbc_12)

'data.frame': 569 obs. of  13 variables:

 $ Y                      : num  1 1 1 1 1 1 1 1 1 1 ...

 $ concavity_se           : num  0.0537 0.0186 0.0383 0.0566 0.0569 ...

 $ compactness_se         : num  0.049 0.0131 0.0401 0.0746 0.0246 ...

 $ fractal_dimension_worst: num  0.1189 0.089 0.0876 0.173 0.0768 ...

 $ symmetry_mean          : num  0.242 0.181 0.207 0.26 0.181 ...

 $ smoothness_mean        : num  0.1184 0.0847 0.1096 0.1425 0.1003 ...

 $ concave.points_se      : num  0.0159 0.0134 0.0206 0.0187 0.0188 ...

 $ texture_mean           : num  10.4 17.8 21.2 20.4 14.3 ...

 $ symmetry_worst         : num  0.46 0.275 0.361 0.664 0.236 ...

 $ smoothness_worst       : num  0.162 0.124 0.144 0.21 0.137 ...

 $ perimeter_se           : num  8.59 3.4 4.58 3.44 5.44 ...

 $ area_worst             : num  2019 1956 1709 568 1575 ...

 $ concave.points_mean    : num  0.1471 0.0702 0.1279 0.1052 0.1043 ...

> # split train(0.8) and test set(0.2)

> set.seed(123) # for reproducibility

> train <- sample(1:nrow(wdbc_12), size=0.8*nrow(wdbc_12), replace=F)

> test <- (-train)

> Y.test <- Y[test]

 




  (2) 로지스틱 회귀모형 적합 (w/ training set)


R의 stats 패키지의 glm() 함수로 training set의 12개 설명변수를 모두 사용하여 로지스틱 회귀모형을 적합하여 보겠습니다. 회귀계수 모수 추정은 최대가능도추정(Maximum Likelihood Estimation)법을 통해서 이루어집니다. 




참고로, GLM 모형의 family 객체에 들어갈 random components 와 link components 에 사용할 수 있는 arguments 는 아래와 같습니다. 


[ Family Objects for Models ]


binomial(link = "logit")  # logistic regression model

poisson(link = "log") # poisson regression model

gaussian(link = "identity") # linear regression model

Gamma(link = "inverse") # gamma regression model

inverse.gaussian(link = "1/mu^2")

quasi(link = "identity", variance = "constant")

quasibinomial(link = "logit")

quasipoisson(link = "log") 




> # train with training set

> glm.fit <- glm(Y ~ ., 

+                data = wdbc_12, 

+                family = binomial(link = "logit")

+                subset = train)

Warning message:

glm.fit: 적합된 확률값들이 0 또는 1 입니다 

> summary(glm.fit)


Call:

glm(formula = Y ~ ., family = binomial(link = "logit"), data = wdbc_12, 

    subset = train)


Deviance Residuals: 

    Min       1Q   Median       3Q      Max  

-1.8232  -0.0662  -0.0074   0.0001   3.8847  


Coefficients:

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

(Intercept)              -36.903201    9.099492  -4.056  0.00005 ***

concavity_se              29.311264   22.989079   1.275 0.202306    

compactness_se          -100.786803   36.755041  -2.742 0.006104 ** 

fractal_dimension_worst   37.094117   43.007209   0.863 0.388407    

symmetry_mean             19.015372   30.163882   0.630 0.528432    

smoothness_mean         -109.619744   74.382548  -1.474 0.140554    

concave.points_se       -123.097155  203.980072  -0.603 0.546192    

texture_mean               0.393147    0.107909   3.643 0.000269 ***

symmetry_worst            13.278414   10.461508   1.269 0.204347    

smoothness_worst         105.722805   35.197258   3.004 0.002667 ** 

perimeter_se               1.033524    0.780904   1.323 0.185670    

area_worst                 0.013239    0.003786   3.496 0.000472 ***

concave.points_mean      104.113327   46.738810   2.228 0.025910 *  

---

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


(Dispersion parameter for binomial family taken to be 1)


    Null deviance: 601.380  on 454  degrees of freedom

Residual deviance:  60.175  on 442  degrees of freedom

AIC: 86.175


Number of Fisher Scoring iterations: 10

 


각 회귀계수의 유의성을 Wald-test()를 통해 해보면 'concave.points_se' 변수가 p-value가 0.54192이므로 유의수준 5% 하에서 [귀무가설(H0): 회귀계수 beta(concave.points_se) = 0]을 채택합니다. 즉, 반응변수에 concave.points_se 변수는 별 영향 혹은 관련이 없다고 판단할 수 있습니다. 따라서 제거를 해도 되겠습니다. 



Wald Test

  • 통계량  
  • 귀무가설(): 회귀계수 
  • 대립가설(): 회귀계수 




  (3) 후진 소거법 (backward elimination)


위의 (2)번 처럼 전체 변수를 포함한 full model (saturated model)에서 시작해서, 가장 유의미하지 않은 변수(가장 큰 p-value 값을 가지는 변수)를 제거하고, 다시 모형을 적합하고, 남은 예측변수가 모두 유의할 때까지 계속 순차적으로 가장 유의미하지 않은 변수를 계속 제거해나가는 후진 소거법(backward elimination)을 사용해서 모델을 적합하겠습니다. 이로써, 모델의 정확도는 어느정도 유지하면서도 해석하기에 용이한 가장 간소한(parsimonious) 모델을 만들어보겠습니다. 



> # Backward Elimination Approach

> # remove 'concave.points_se' variable

> glm.fit.2 <- glm(Y ~ concave.points_mean + area_worst + 

+                    perimeter_se + smoothness_worst + symmetry_worst +

+                    texture_mean + smoothness_mean + symmetry_mean + 

+                    fractal_dimension_worst + compactness_se + concavity_se, 

+                  data = wdbc_12, 

+                  family = binomial(link = "logit")

+                  subset = train)

Warning message:

glm.fit: 적합된 확률값들이 0 또는 1 입니다 

> summary(glm.fit.2)


Call:

glm(formula = Y ~ concave.points_mean + area_worst + perimeter_se + 

    smoothness_worst + symmetry_worst + texture_mean + smoothness_mean + 

    symmetry_mean + fractal_dimension_worst + compactness_se + 

    concavity_se, family = binomial(link = "logit"), data = wdbc_12, 

    subset = train)


Deviance Residuals: 

    Min       1Q   Median       3Q      Max  

-1.8363  -0.0645  -0.0076   0.0001   3.8297  


Coefficients:

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

(Intercept)              -37.182194    8.981711  -4.140 0.0000348 ***

concave.points_mean       87.773808   37.117336   2.365  0.018041 *  

area_worst                 0.013785    0.003734   3.692  0.000222 ***

perimeter_se               0.894502    0.747623   1.196  0.231517    

smoothness_worst         104.169968   35.281001   2.953  0.003151 ** 

symmetry_worst            15.057387   10.048439   1.498  0.134009    

texture_mean               0.385946    0.105254   3.667  0.000246 ***

smoothness_mean         -103.800173   73.953483  -1.404  0.160442    

symmetry_mean             10.952803   26.586484   0.412  0.680362    

fractal_dimension_worst   42.376803   42.322307   1.001  0.316688    

compactness_se           -98.379602   35.688499  -2.757  0.005840 ** 

concavity_se              18.847381   13.973683   1.349  0.177409    

---

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


(Dispersion parameter for binomial family taken to be 1)


    Null deviance: 601.380  on 454  degrees of freedom

Residual deviance:  60.551  on 443  degrees of freedom

AIC: 84.551


Number of Fisher Scoring iterations: 10


> # remove 'symmetry_mean' variable

> glm.fit.3 <- glm(Y ~ concave.points_mean + area_worst + 

+                    perimeter_se + smoothness_worst + symmetry_worst +

+                    texture_mean + smoothness_mean + 

+                    fractal_dimension_worst + compactness_se + concavity_se, 

+                  data = wdbc_12, 

+                  family = binomial(link = "logit")

+                  subset = train)

Warning message:

glm.fit: 적합된 확률값들이 0 또는 1 입니다 

> summary(glm.fit.3)


Call:

glm(formula = Y ~ concave.points_mean + area_worst + perimeter_se + 

    smoothness_worst + symmetry_worst + texture_mean + smoothness_mean + 

    fractal_dimension_worst + compactness_se + concavity_se, 

    family = binomial(link = "logit"), data = wdbc_12, subset = train)


Deviance Residuals: 

    Min       1Q   Median       3Q      Max  

-1.7842  -0.0630  -0.0083   0.0001   3.7588  


Coefficients:

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

(Intercept)             -36.121685   8.504448  -4.247 0.0000216 ***

concave.points_mean      89.488263  37.243081   2.403  0.016269 *  

area_worst                0.013329   0.003532   3.774  0.000161 ***

perimeter_se              0.976273   0.722740   1.351  0.176762    

smoothness_worst        100.865240  33.783402   2.986  0.002830 ** 

symmetry_worst           17.885605   7.503047   2.384  0.017136 *  

texture_mean              0.382870   0.103885   3.686  0.000228 ***

smoothness_mean         -95.335007  70.234774  -1.357  0.174662    

fractal_dimension_worst  40.176580  41.722080   0.963  0.335569    

compactness_se          -97.950160  35.684347  -2.745  0.006053 ** 

concavity_se             19.800333  13.655596   1.450  0.147064    

---

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


(Dispersion parameter for binomial family taken to be 1)


    Null deviance: 601.380  on 454  degrees of freedom

Residual deviance:  60.721  on 444  degrees of freedom

AIC: 82.721


Number of Fisher Scoring iterations: 10


> # remove 'fractal_dimension_worst' variable

> glm.fit.4 <- glm(Y ~ concave.points_mean + area_worst + 

+                    perimeter_se + smoothness_worst + symmetry_worst +

+                    texture_mean + smoothness_mean + 

+                    compactness_se + concavity_se, 

+                  data = wdbc_12, 

+                  family = binomial(link = "logit")

+                  subset = train)

Warning message:

glm.fit: 적합된 확률값들이 0 또는 1 입니다 

> summary(glm.fit.4)


Call:

glm(formula = Y ~ concave.points_mean + area_worst + perimeter_se + 

    smoothness_worst + symmetry_worst + texture_mean + smoothness_mean + 

    compactness_se + concavity_se, family = binomial(link = "logit"), 

    data = wdbc_12, subset = train)


Deviance Residuals: 

    Min       1Q   Median       3Q      Max  

-1.6650  -0.0653  -0.0079   0.0001   3.9170  


Coefficients:

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

(Intercept)          -33.656002    7.809704  -4.310 0.0000164 ***

concave.points_mean   96.137001   36.867939   2.608  0.009118 ** 

area_worst             0.012905    0.003435   3.757  0.000172 ***

perimeter_se           0.676729    0.642490   1.053  0.292208    

smoothness_worst     109.530724   32.239695   3.397  0.000680 ***

symmetry_worst        19.724732    7.543696   2.615  0.008930 ** 

texture_mean           0.381544    0.103567   3.684  0.000230 ***

smoothness_mean     -101.037553   70.839248  -1.426  0.153784    

compactness_se       -80.383972   30.461425  -2.639  0.008318 ** 

concavity_se          20.771502   14.372162   1.445  0.148385    

---

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


(Dispersion parameter for binomial family taken to be 1)


    Null deviance: 601.380  on 454  degrees of freedom

Residual deviance:  61.657  on 445  degrees of freedom

AIC: 81.657


Number of Fisher Scoring iterations: 10


> # remove 'perimeter_se' variable

> glm.fit.5 <- glm(Y ~ concave.points_mean + area_worst + smoothness_worst + 

+                    symmetry_worst + texture_mean + smoothness_mean + 

+                    compactness_se + concavity_se, 

+                  data = wdbc_12, 

+                  family = binomial(link = "logit")

+                  subset = train)

Warning message:

glm.fit: 적합된 확률값들이 0 또는 1 입니다 

> summary(glm.fit.5)


Call:

glm(formula = Y ~ concave.points_mean + area_worst + smoothness_worst + 

    symmetry_worst + texture_mean + smoothness_mean + compactness_se + 

    concavity_se, family = binomial(link = "logit"), data = wdbc_12, 

    subset = train)


Deviance Residuals: 

    Min       1Q   Median       3Q      Max  

-1.7487  -0.0602  -0.0079   0.0002   4.2383  


Coefficients:

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

(Intercept)         -33.156014   7.674110  -4.321 0.0000156 ***

concave.points_mean  98.900172  35.529247   2.784   0.00538 ** 

area_worst            0.013819   0.003339   4.139 0.0000349 ***

smoothness_worst    105.448378  32.444984   3.250   0.00115 ** 

symmetry_worst       17.412883   6.654537   2.617   0.00888 ** 

texture_mean          0.385751   0.102180   3.775   0.00016 ***

smoothness_mean     -89.410469  70.307291  -1.272   0.20348    

compactness_se      -78.009472  29.946520  -2.605   0.00919 ** 

concavity_se         24.607453  13.153360   1.871   0.06137 .  

---

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


(Dispersion parameter for binomial family taken to be 1)


    Null deviance: 601.380  on 454  degrees of freedom

Residual deviance:  62.777  on 446  degrees of freedom

AIC: 80.777


Number of Fisher Scoring iterations: 10


> # remove 'smoothness_mean' variable

> glm.fit.6 <- glm(Y ~ concave.points_mean + area_worst + smoothness_worst + 

+                    symmetry_worst + texture_mean +  

+                    compactness_se + concavity_se, 

+                  data = wdbc_12, 

+                  family = binomial(link = "logit"), 

+                  subset = train)

Warning message:

glm.fit: 적합된 확률값들이 0 또는 1 입니다 

> summary(glm.fit.6)


Call:

glm(formula = Y ~ concave.points_mean + area_worst + smoothness_worst + 

    symmetry_worst + texture_mean + compactness_se + concavity_se, 

    family = binomial(link = "logit"), data = wdbc_12, subset = train)


Deviance Residuals: 

    Min       1Q   Median       3Q      Max  

-1.8607  -0.0640  -0.0080   0.0002   4.0912  


Coefficients:

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

(Intercept)         -39.164619   6.748646  -5.803 0.0000000065 ***

concave.points_mean  73.504634  27.154772   2.707     0.006792 ** 

area_worst            0.015465   0.003248   4.762 0.0000019175 ***

smoothness_worst     78.962941  23.746808   3.325     0.000884 ***

symmetry_worst       17.429475   6.497808   2.682     0.007310 ** 

texture_mean          0.423250   0.099815   4.240 0.0000223171 ***

compactness_se      -79.556564  29.747330  -2.674     0.007486 ** 

concavity_se         27.279384  12.093534   2.256     0.024089 *  

---

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


(Dispersion parameter for binomial family taken to be 1)


    Null deviance: 601.380  on 454  degrees of freedom

Residual deviance:  64.446  on 447  degrees of freedom

AIC: 80.446


Number of Fisher Scoring iterations: 10

 



(2)의 12개 설명변수를 사용한 full model에서 시작해서 (3)의 후진소거법(backward elimination)에 의해 5개의 변수를 하나씩 순차적으로 제거하면서 적합한 모형의 결과를 표로 정리하면 아래와 같습니다. 


자유도(df) 대비 이탈도(deviance)와의 차이가 거의 없으므로 후진소거법에 의해 변수를 제거해도 모델에는 별 영향이 없음을 알 수 있습니다. 그리고 AIC(= -2(로그가능도 - 모형에 있는 모수 개수)) 기준에 의하면 6번째 모형 (concave.points_se, symmetry_mean, fractal_dimension_worst, perimeter_se, smoothness_mean 제거한 모형)이 가장 우수한 모형임을 알 수 있습니다.(AIC 값이 가장 작으므로)



6번째 모델의 '귀무가설 H0: 모든 회귀계수 = 0 (무가치한 모형)'에 대한 가능도비 검정(likelihood ration test) 결과, 


 

> # likelihood ratio test

> LR_statistic = 601.380 - 64.446

> LR_statistic

[1] 536.934

pchisq(LR_statistic, df = (454-447), lower.tail = F)

[1] 9.140339e-112




이므로 적어도 하나의 설명변수는 반응변수를 예측하는데 유의미한 관련이 있다고 볼 수 있다. 



그리고 각 설명변수별 Wald-test 결과 p-value(Pr(>|z|) 가 모두 0.05보다 작으므로, 유의수준 5% 하에서 모두 유의미하다(즉, 반응변수와 관련이 있다)고 판단할 수 있다. 



summary(glm.fit.6)


Call:

glm(formula = Y ~ concave.points_mean + area_worst + smoothness_worst + 

    symmetry_worst + texture_mean + compactness_se + concavity_se, 

    family = binomial(link = "logit"), data = wdbc_12, subset = train)


Deviance Residuals: 

    Min       1Q   Median       3Q      Max  

-1.8607  -0.0640  -0.0080   0.0002   4.0912  


Coefficients:

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

(Intercept)         -39.164619   6.748646  -5.803 0.0000000065 ***

concave.points_mean  73.504634  27.154772   2.707     0.006792 ** 

area_worst            0.015465   0.003248   4.762 0.0000019175 ***

smoothness_worst     78.962941  23.746808   3.325     0.000884 ***

symmetry_worst       17.429475   6.497808   2.682     0.007310 ** 

texture_mean          0.423250   0.099815   4.240 0.0000223171 ***

compactness_se      -79.556564  29.747330  -2.674     0.007486 ** 

concavity_se         27.279384  12.093534   2.256     0.024089 *  

---

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


(Dispersion parameter for binomial family taken to be 1)


    Null deviance: 601.380  on 454  degrees of freedom

Residual deviance:  64.446  on 447  degrees of freedom

AIC: 80.446


Number of Fisher Scoring iterations: 10

 




로지스틱 회귀계수의 95% 신뢰구간은 confint() 함수를 통해 구할 수 있습니다. 



> confint(glm.fit.6)

Waiting for profiling to be done...

                            2.5 %       97.5 %

(Intercept)         -5.483332e+01 -27.92705624

concave.points_mean  2.353047e+01 131.95158315

area_worst           9.906713e-03   0.02282079

smoothness_worst     3.617243e+01 131.06234268

symmetry_worst       6.422548e+00  32.73189772

texture_mean         2.428527e-01   0.64064918

compactness_se      -1.411215e+02 -15.15533223

concavity_se        -1.259598e+01  49.25042906

 




  (4) 모델 성능 평가 (w/ test set) : confusion matrix, ROC curve


모델 성능 평가는 훈련용 데이터셋이 아니라 반드시 테스트셋(test set)을 대상으로 진행해야 합니다! 


(a) confusion matrix를 이용한 분류 정확도, 오분류율을 구해보고, (b) ROC 곡선도 그려보겠습니다. 


(a) Confusion matrix, accuracy, test error rate


cutoff 기준을 이면 이면 으로 예측을 하면, 양성 72개는 모두 양성으로 전부 정확하게 분류를 하였으며, 악성 42개에 대해서는 41개를 악성으로 정확하게 분류하고 1개는 양성으로 오분류 하였습니다. 


따라서 정확도(accuracy) = (TP + TN) / N = (72 + 41)/ 114 = 0.9912281 로서 매우 높게 나왔습니다. 유방암 예측(악성, 양성 분류) 로지스틱 회귀모형이 잘 적합되었네요. 



> # calculate the probabilities on test set

> glm.probs <- predict(glm.fit.6, wdbc_12[test,], type="response")

> glm.probs[1:20]

                6                 8                10                11                12 

0.986917125953939 0.992018931502162 0.998996856414432 0.992792364235021 0.999799557490392 

               13                25                35                39                54 

0.999938548975209 1.000000000000000 0.999981682993826 0.003031968742568 0.999978684326258 

               61                67                72                77                86 

0.000024521469882 0.000434266606978 0.000000006158534 0.002210951131364 0.999999952448192 

               92                98               106               108               115 

0.923496652258054 0.000007994931448 0.998688178459914 0.000576845750543 0.000149114805822 

> # compute the predictions with the threshold of 0.5

> glm.pred <- rep(0, nrow(wdbc_12[test,]))

> glm.pred[glm.probs > .5] = 1

> table(Y.test, glm.pred)

      glm.pred

Y.test  0  1

     0 72  0

     1  1 41

> mean(Y.test == glm.pred) # accuracy

[1] 0.9912281

> mean(Y.test != glm.pred) # test error rate

[1] 0.00877193

 



(b) ROC 곡선 (ROC curve), AUC (The Area Under an ROC Curve)


cutoff 별  'False Positive Rate(1-Specificity)', 'True Positive Rate(Sensitivity)' 값을 연속적으로 계산하여 ROC 곡선을 그려보면 아래와 같습니다. 대각선이 random guessing 했을 때를 의미하며, 대각선으로부터 멀리 떨어져서 좌측 상단으로 곡선이 붙을 수록 더 잘 적합된 모형이라고 해석할 수 있습니다. 


AUC(Area Under an ROC Curve)는 ROC 곡선 아래 부분의 면적의 합이며, 1이면 완벽한 모델, 0.5면 random guessing한 모델의 수준입니다. AUC가 0.9 이상이면 매우 우수한 모형이라고 할 수 있습니다. 


이번 로지스틱회귀모형은 ROC 곡선을 봐서도 그렇고, AUC가 0.9933로서 아주 잘 적합된 모형이네요. 



> # ROC curve

> install.packages("ROCR")

> library(ROCR)

> pr <- prediction(glm.probs, Y.test)

> prf <- performance(pr, measure = "tpr", x.measure = "fpr")

> plot(prf, main="ROC Curve")




> # AUC (The Area Under an ROC Curve)

> auc <- performance(pr, measure = "auc")

> auc <- auc@y.values[[1]]

> auc

[1] 0.9933862

> detach(wdbc)





  (5) 모델 해석


최종 적합된 로지스틱 회귀모형식은 아래와 같습니다. 

(워드 수식입력기에서 '_'를 입력하니 자꾸 '.'으로 자동으로 바뀌는 바람에 변수이름 구분자 '_'가 전부 '.'로 잘못 입력되었습니다....)




위의 최종 모형으로 부터 다른 설명변수가 고정(통제)되었을 때 설명변수 가 한단위 증가할 때 유방암 악성(M)일 확률의 오즈(odds)는 exp() 배만큼 증가한다는 뜻입니다. 


가령, texture_mean 의 회귀계수가0.42 이므로 이를 해석하자면, 다른 설명변수들이 고정되었을 때 texture_mean이 한 단위 증가할 때 유방암 악성일 확률의 오즈(odds)는 exp(0.42)=1.52 배 증가한다는 뜻입니다. 


악성 유방암일 확률은  

식을 사용해서 구할 수 있습니다. (위의 식 beta_hat 에 추정한 회귀계수를 대입해고 x 값에 관측값을 대입해서 풀면 됩니다)


많은 도움이 되었기를 바랍니다. 


그동안 4번에 걸쳐서 로지스틱 회귀모형 분석에 대해 연재를 했는데요, 혹시 포스팅에 잘못된 부분이나 궁금한 부분이 있으면 댓글로 남겨주시면 수정, 보완하도록 하겠습니다. 


이번 포스팅이 도움이 되었다면 아래의 '공감~'를 꾹 눌러주세요. ^^



Posted by R Friend R_Friend

이번 포스팅은 로지스틱 회귀분석을 통한 유방암 예측(분류) 포스팅의 세번째 순서로서 '1차 변수 선택 및 목표변수와 설명변수 간 관계 분석' 차례입니다. 



[로지스틱 회귀분석을 통한 유방암 예측(분류) 포스팅 순서]

  1. WDBC(Wisconsin Diagnostic Breast Cancer) dataset 소개 및 분석 목적과 방향 설정
  2. 탐색적 데이터 분석 및 전처리
  3. 1차 변수 선택 및 목표변수와 설명변수 간 관계 분석
  4. 로지스틱 회귀모형 적합 및 모델 평가, 해석



사실, 이번 포스팅도 '2. 탐색적 데이터 분석 및 전처리'에 포함시켜도 되는데요, 두번째 포스팅이 미친듯이 내용이 많고 포스팅이 길어져서 스크롤 압박이 너무 심했던지라 가독성을 위해서 세번째 포스팅으로 분리를 했습니다. 



  (1) t-test 를 통한 1차 변수 선별


진단 결과(악성 M, 양성 B) 그룹별로 설명변수 간의 차이가 존재하는지 t-test 검정을 해보고, p-value가 0.05 를 초과하는 설명변수는 1차 제거하는 것으로 하겠습니다. 


특히 연속형 설명변수 개수가 상당히 많고, 종속변수가 2개의 범주(class)를 가지는 경우 t-test 를 활용해서 1차로 별 효용성이 없는 설명변수를 screening out 하는데 활용할 수 있습니다. 



> # Variable selection (1st round)

> # Multiple t-tests of explanatory variables between diagnosis

> X_names <- names(data.frame(X_3))

> X_names

 [1] "texture_mean"            "smoothness_mean"         "concave.points_mean"    

 [4] "symmetry_mean"           "fractal_dimension_mean"  "texture_se"             

 [7] "perimeter_se"            "smoothness_se"           "compactness_se"         

[10] "concavity_se"            "concave.points_se"       "symmetry_se"            

[13] "fractal_dimension_se"    "area_worst"              "smoothness_worst"       

[16] "symmetry_worst"          "fractal_dimension_worst"

> t.test_p.value_df <- data.frame() # blank data.frame for saving

> for (i in 1:length(X_names)) {

+   t.test_p.value <- t.test(wdbc[,X_names[i]] ~ wdbc$diagnosis, var.equal = TRUE)$p.value

+   

+   t.test_p.value_df[i,1] <- X_names[i]

+   t.test_p.value_df[i,2] <- t.test_p.value

+ }

> colnames(t.test_p.value_df) <- c("x_var_name", "p.value")

> t.test_p.value_df

                x_var_name       p.value

1             texture_mean  4.058636e-25

2          smoothness_mean  1.051850e-18

3      concave.points_mean 7.101150e-116

4            symmetry_mean  5.733384e-16

5   fractal_dimension_mean  7.599368e-01

6               texture_se  8.433320e-01

7             perimeter_se  1.651905e-47

8            smoothness_se  1.102966e-01

9           compactness_se  9.975995e-13

10            concavity_se  8.260176e-10

11       concave.points_se  3.072309e-24

12             symmetry_se  8.766418e-01

13    fractal_dimension_se  6.307355e-02

14              area_worst  2.828848e-97

15        smoothness_worst  6.575144e-26

16          symmetry_worst  2.951121e-25

17 fractal_dimension_worst  2.316432e-15

 



위의 17개 설명변수별 t-test 결과를 좀더 보기에 편하도록 p.value를 기준으로 작은 것부터 큰 순서대로 dplyr 패키지의 arrange() 함수를 사용하여 정렬해보겠습니다. 



> # sorting by p.value in ascending order

> install.packages("dplyr")

> library(dplyr)

> arrange(t.test_p.value_df, p.value)

                x_var_name       p.value

1      concave.points_mean 7.101150e-116

2               area_worst  2.828848e-97

3             perimeter_se  1.651905e-47

4         smoothness_worst  6.575144e-26

5           symmetry_worst  2.951121e-25

6             texture_mean  4.058636e-25

7        concave.points_se  3.072309e-24

8          smoothness_mean  1.051850e-18

9            symmetry_mean  5.733384e-16

10 fractal_dimension_worst  2.316432e-15

11          compactness_se  9.975995e-13

12            concavity_se  8.260176e-10

13    fractal_dimension_se  6.307355e-02

14           smoothness_se  1.102966e-01

15  fractal_dimension_mean  7.599368e-01

16              texture_se  8.433320e-01

17             symmetry_se  8.766418e-01

 



t-test의 p-value가 0.05 보다 큰 값을 가지는 설명변수인 'symmetry_se', 'texture_se', 'fractal_dimension_mean', 'smoothness_se', 'fractal_dimension_se' 의 5개 설명변수는 1차로 제거하고, 나머지 12개 설명변수만 로지스틱 회귀모형 적합하는데 사용하도록 하겠습니다. (말그대로 1차 선별이구요, 나중에 후진소거법으로 추가로 더 변수 선택할 예정입니다)




위의 1차 선별된 12개 설명변수만을 포함한 데이터프레임 X_4를 만들었습니다. 



> # select x_variables only if p.value < 0.05

> t.test_filtered <- t.test_p.value_df$p.value < 0.05

> X_names_filtered <- X_names[t.test_filtered]

> X_4 <- data.frame(X_3[, X_names_filtered])

> str(X_4)

'data.frame': 569 obs. of  12 variables:

 $ texture_mean           : num  -2.072 -0.353 0.456 0.254 -1.151 ...

 $ smoothness_mean        : num  1.567 -0.826 0.941 3.281 0.28 ...

 $ concave.points_mean    : num  2.53 0.548 2.035 1.45 1.427 ...

 $ symmetry_mean          : num  2.21557 0.00139 0.93886 2.86486 -0.00955 ...

 $ perimeter_se           : num  2.831 0.263 0.85 0.286 1.272 ...

 $ compactness_se         : num  1.3157 -0.6923 0.8143 2.7419 -0.0485 ...

 $ concavity_se           : num  0.723 -0.44 0.213 0.819 0.828 ...

 $ concave.points_se      : num  0.66 0.26 1.42 1.11 1.14 ...

 $ area_worst             : num  2 1.89 1.46 -0.55 1.22 ...

 $ smoothness_worst       : num  1.307 -0.375 0.527 3.391 0.22 ...

 $ symmetry_worst         : num  2.748 -0.244 1.151 6.041 -0.868 ...

 $ fractal_dimension_worst: num  1.935 0.281 0.201 4.931 -0.397 ...





  (2) 목표변수와 설명변수 간 관계 분석 (시각화)


(2-1) 박스 그림 (Box Plot)


다음으로 12개 설명변수에 대해서 진단결과(악성: M, 1 vs. 양성: B, 0) 별로 집단을 분리해서 박스 그림(Box Plot)을 그려서 비교를 해보겠습니다. 


그래프로 보기에 편리하도록 표준화한 데이터셋을 p.value 기준으로 변수의 순서를 정렬한 후에 Y값(diagnosis)과 합쳐서 새로운 데이터셋을 만들었습니다. 



> # sorting by p.value in descending order

> # t.test_p.value_df.sorted <- arrange(t.test_p.value_df[t.test_filtered,], p.value)

> # t.test_p.value_df.sorted

> t.test_p.value_df.sorted_2 <- arrange(t.test_p.value_df[t.test_filtered,], desc(p.value))

> t.test_p.value_df.sorted_2

                x_var_name       p.value

1             concavity_se  8.260176e-10

2           compactness_se  9.975995e-13

3  fractal_dimension_worst  2.316432e-15

4            symmetry_mean  5.733384e-16

5          smoothness_mean  1.051850e-18

6        concave.points_se  3.072309e-24

7             texture_mean  4.058636e-25

8           symmetry_worst  2.951121e-25

9         smoothness_worst  6.575144e-26

10            perimeter_se  1.651905e-47

11              area_worst  2.828848e-97

12     concave.points_mean 7.101150e-116

> x_names_sorted <- t.test_p.value_df.sorted_2$x_var_name

> x_names_sorted

 [1] "concavity_se"            "compactness_se"          "fractal_dimension_worst"

 [4] "symmetry_mean"           "smoothness_mean"         "concave.points_se"      

 [7] "texture_mean"            "symmetry_worst"          "smoothness_worst"       

[10] "perimeter_se"            "area_worst"              "concave.points_mean"    

> X_5 <- X_4[x_names_sorted] # rearrange column order for plotting below 

> head(X_5,2)

  concavity_se compactness_se fractal_dimension_worst symmetry_mean smoothness_mean

1    0.7233897      1.3157039               1.9353117   2.215565542       1.5670875

2   -0.4403926     -0.6923171               0.2809428   0.001391139      -0.8262354

  concave.points_se texture_mean symmetry_worst smoothness_worst perimeter_se area_worst

1         0.6602390   -2.0715123      2.7482041        1.3065367    2.8305403   1.999478

2         0.2599334   -0.3533215     -0.2436753       -0.3752817    0.2630955   1.888827

  concave.points_mean

1           2.5302489

2           0.5476623

> #-----

> # combine Y and X

> wdbc_2 <- data.frame(Y, X_5)

> str(wdbc_2)

'data.frame': 569 obs. of  13 variables:

 $ Y                      : num  1 1 1 1 1 1 1 1 1 1 ...

 $ concavity_se           : num  0.723 -0.44 0.213 0.819 0.828 ...

 $ compactness_se         : num  1.3157 -0.6923 0.8143 2.7419 -0.0485 ...

 $ fractal_dimension_worst: num  1.935 0.281 0.201 4.931 -0.397 ...

 $ symmetry_mean          : num  2.21557 0.00139 0.93886 2.86486 -0.00955 ...

 $ smoothness_mean        : num  1.567 -0.826 0.941 3.281 0.28 ...

 $ concave.points_se      : num  0.66 0.26 1.42 1.11 1.14 ...

 $ texture_mean           : num  -2.072 -0.353 0.456 0.254 -1.151 ...

 $ symmetry_worst         : num  2.748 -0.244 1.151 6.041 -0.868 ...

 $ smoothness_worst       : num  1.307 -0.375 0.527 3.391 0.22 ...

 $ perimeter_se           : num  2.831 0.263 0.85 0.286 1.272 ...

 $ area_worst             : num  2 1.89 1.46 -0.55 1.22 ...

 $ concave.points_mean    : num  2.53 0.548 2.035 1.45 1.427 ...

 



그 다음으로 reshape2 패키지의 melt() 함수를 사용해서 데이터를 세로로 길게 재구조화한 다음에, ggplot2패키지의 ggplot() + geom_boxplot() 함수를 사용하여 박스 그림을 그렸습니다. 



> #-----

> # Box plot of X per Y(M, B)

> install.packages("reshape2")

> library(reshape2)

> wdbc_2_melt <- melt(wdbc_2, id.var = "Y")

> install.packages("ggplot2")

> library(ggplot2)

> ggplot(data = wdbc_2_melt[wdbc_2_melt$value < 3,], aes(x=variable, y=value)) + 

+   geom_boxplot(aes(fill=as.factor(Y))) +

+   theme_bw() + # white background

+   coord_flip() # flip the x & y-axis

 





이 박스그림은 위의 t-test 결과 p-value 가 작았던 설명변수 순서대로 위에서 부터 아래로 그려진 것인데요, 역시나 p-value가 작을 수록 박스그림에서 보면 진단결과 Malignant(악성)인 그룹과 Benign(양성) 간의 차이가 크게 나고 있음을 눈으로 재확인할 수 있습니다. 





(2-2) 산점도 그림 (Scatter Plot)


다음으로 p-value가 가장 작게 나왔던 상위 변수들 중에서 일부인 'concave.points_mean', 'area_worst', 'texture_mean'를 조합해서 사용해서 예시로 산점도를 그려보았습니다. 이때 진단결과(악성: 'M', 1, vs. 양성: 'B', 0) 별로 색깔을 다르게 해서 산점도를 그렸습니다. 


아래 예시의 2개 산점도를 봐서는 로지스틱 회귀모형으로 분류모델을 적합해도 잘 될 것으로 예상이 되네요. 



> # scatter plot of x=concave.points_mean, y=area_worst

> ggplot(data=wdbc_2, aes(x=concave.points_mean, y=area_worst, colour=as.factor(Y), alpha=0.5)) +

+   geom_point(shape=19, size=3) +

+   ggtitle("Scatter Plot of concave.points_mean & area_worst by Y")




> ggplot(data=wdbc_2, aes(x=area_worst, y=texture_mean, colour=as.factor(Y), alpha=0.5)) +

+   geom_point(shape=19, size=3) +

+   ggtitle("Scatter Plot of area_worst & texture_mean by Y")




이상으로 세번째 포스팅인 't-test를 활용한 1차 변수 선별 및 목표변수와 설명변수간 관계 (탐색적) 분석'을 마치겠습니다. 


다음번 포스팅은 마지막인 '로지스틱 회귀모형 적합 및 모델 평가, 해석'에 대해서 다루겠습니다. 


이번 포스팅이 도움이 되었다면 아래의 '공감~'를 꾸욱 눌러주세요. ^^



Posted by R Friend R_Friend