지난번 WDBC(Wisconsin Diagnostic Breast Cancer) dataset 소개 및 분석 목적과 방향 설정 포스팅에 이어서, 이번 포스팅은 두번째 순서로 'WDBC 데이터셋에 대한 탐색적 데이터 분석과 전처리'에 대해서 알아보겠습니다. 



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

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



지난번 첫번째 포스팅에서 데이터 셋 가져와서 DataFrame으로 만들었을 때 사용했던 R코드는 아래와 같습니다. 



rm(list=ls())

options(scipen=30)


# data loading: WDBC (Wisconsin Diagnostic Breast Cancer)

library(data.table)

library(curl)


url <- c("https://archive.ics.uci.edu/ml/machine-learning-databases/breast-cancer-wisconsin/wdbc.data")

wdbc <- data.frame(fread(url))

str(wdbc)


# column names

colnames(wdbc) <- c("id", "diagnosis", "radius_mean", "texture_mean", 

                    "perimeter_mean", "area_mean", "smoothness_mean", "compactness_mean",

                    "concavity_mean", "concave.points_mean", "symmetry_mean", "fractal_dimension_mean", 

                    "radius_se", "texture_se", "perimeter_se", "area_se", 

                    "smoothness_se", "compactness_se", "concavity_se", "concave.points_se", 

                    "symmetry_se", "fractal_dimension_se", "radius_worst", "texture_worst", 

                    "perimeter_worst", "area_worst", "smoothness_worst", "compactness_worst", 

                    "concavity_worst", "concave.points_worst", "symmetry_worst", "fractal_dimension_worst")


str(wdbc)

head(wdbc, 2)

 



이렇게 불러온 WDBC DataFrame에 대해서 


(1) 결측값 확인 및 처리

(2) 중복 데이터 확인 및 처리

(3) 목표변수 범주/계급 구성 분포 확인 및 처리

(4) 설명변수 간 다중공선성 확인 및 처리

(5) 표준화, 척도 변환


의 순서로 탐색적 데이터 분석(EDA: Exploratory Data Analysis) 및 데이터 전처리(Data Preprocessing)를 해보겠습니다. 



 (1) 결측값 (Missing Value) 확인 및 처리


is.na() 함수를 사용해서 결측값 여부를 확인 결과 다행히도 결측값이 없네요.  데이터셋 설명된 사이트에 보면 처음에는 결측값이 있었다고 하는데요, 지금 사용하는 데이터셋은 결측값을 다 제거한 자료네요. 


혹시 결측값이 있다면 http://rfriend.tistory.com/34 포스팅을 참고하여 분석 목적에 맞게 결측값을 처리하시기 바랍니다. 



> attach(wdbc)

The following objects are masked from wdbc (pos = 3):


    area_mean, area_se, area_worst, compactness_mean, compactness_se, compactness_worst,

    concave.points_mean, concave.points_se, concave.points_worst, concavity_mean,

    concavity_se, concavity_worst, diagnosis, fractal_dimension_mean,

    fractal_dimension_se, fractal_dimension_worst, id, perimeter_mean, perimeter_se,

    perimeter_worst, radius_mean, radius_se, radius_worst, smoothness_mean,

    smoothness_se, smoothness_worst, symmetry_mean, symmetry_se, symmetry_worst,

    texture_mean, texture_se, texture_worst


> # check missing value

> colSums(is.na(wdbc)) # no missing value, good

                     id               diagnosis             radius_mean            texture_mean 

                      0                       0                       0                       0 

         perimeter_mean               area_mean         smoothness_mean        compactness_mean 

                      0                       0                       0                       0 

         concavity_mean     concave.points_mean           symmetry_mean  fractal_dimension_mean 

                      0                       0                       0                       0 

              radius_se              texture_se            perimeter_se                 area_se 

                      0                       0                       0                       0 

          smoothness_se          compactness_se            concavity_se       concave.points_se 

                      0                       0                       0                       0 

            symmetry_se    fractal_dimension_se            radius_worst           texture_worst 

                      0                       0                       0                       0 

        perimeter_worst              area_worst        smoothness_worst       compactness_worst 

                      0                       0                       0                       0 

        concavity_worst    concave.points_worst          symmetry_worst fractal_dimension_worst 

                      0                       0                       0                       0

> sum(is.na(wdbc))

[1] 0





  (2) 중복 데이터 (Duplicated Data) 확인 및 처리


다음으로 duplicated() 함수를 사용해서 중복 데이터가 있는지 확인해보았는데요, 다행히도 중복 데이터는 없네요. ^^


혹시 중복 데이터가 있다면 유일한 값을 선별하는 방법은 http://rfriend.tistory.com/165 포스팅을 참고하시기 바랍니다. 



> # check duplicated data

> sum(duplicated(wdbc)) # no duplicated data, good~!

[1] 0

 




  (3) 목표변수 범주/계급 구성 분포(class distribution) 확인 및 처리


세번째로 목표변수(반응변수, 종속변수) diagnosis 의 범주(category, class)의 분포를 살펴보겠습니다. 

table() 함수로 각 계급별 빈도 수 집계, margin.table() 함수로 총 계 계산, prop.table() 함수로 비율을 구하였습니다. 


다행히도 양성(Benign)과 악성(Malignant) 비율이 62.7% vs. 37.2% 로서, 어느 한쪽으로 심하게 불균형이 아니고 양쪽으로 적정하게 배분이 되어 있네요. (데이터를 공개한 교수님들이 참 친절하십니다!)



> # check Y: diagnosis

> # The diagnosis of breast tissues (M = malignant, B = benign)

> table(diagnosis); cat("total :", margin.table(table(diagnosis)))

diagnosis

  B   M 

357 212 

total : 569

> prop.table(table(diagnosis)) # Benign 62.7% vs. Malignant 37.2%

diagnosis

        B         M 

0.6274165 0.3725835

 


만약, 우리가 관심있는 악성(M)의 비율이 1% 미만인 심하게 불균형한 자료(imbalanced dataset)인 경우에는 관측치 빈도수가 클 경우 majority undersampling (정보 손실의 단점 있음)을 하거나, 아니면 SMOTE(Synthetic Minority Over-Sampling TEchnique) 알고리즘을 사용하여 minority oversampling 을 하거나, 혹은 모델 알고리즘에서 가중치를 목표변수 계급의 구성비를 기준으로 조정하는 방법을 사용하면 됩니다. 




  (4) 설명변수 간 다중공선성(Multicollinearity) 확인 및 처리


네번째로, 설명변수 간 강한 상관관계가 있는 다중공선성(Multicolleniarity)이 존재하는지를 


 - 산점도 (Scatter plot)

 - 상관계수 (Correlation Coefficient)

 - 분산팽창지수 (VIF: Variance Inflation Factor)


를 사용하여 확인해보겠습니다. 회귀모형에서는 설명변수간 독립성(independace)을 가정하므로, 독립성 가정 만족 여부를 만족하지 않을 경우 조치하기 위함입니다. 


PerformanceAnalytics 패키지의 chart.Correlation() 함수를 사용하면 한꺼번에 편리하게 각 변수의 히스토그램, 변수들 간의 산점도와 상관계수를 그려볼 수 있습니다. 


코딩하기에 편하도록 설명변수만 따로 떼어서 X 라는 데이터프레임을 만들었습니다. 그리고 설명변수가 30개씩이나 되어 한꺼번에 그래프를 그릴 경우 그래프 결과나 너무나 작고 조밀하게 나와서 보기에 불편하므로, 편의상 'mean' 측정 10개 변수, 'se(standard error' 측정 10개 변수, 'worst' 측정 10개 변수의 3개 그룹으로 나누어서 그래프를 그려보았습니다. (변수 간 상관성이 높게 나온 결과가 3개 그룹간 유사합니다). 


분석 결과 매우 높은 상관관계(상관계수 0.9 이상)를 보이는 설명변수들이 여러개 존재하므로, 다중공선성을 강하게 의심할 수 있습니다. 



> # split Y, X

> Y <- ifelse(wdbc$diagnosis == 'M', 1, 0)

> X <- wdbc[,c(3:32)]

> names(X)

 [1] "radius_mean"             "texture_mean"            "perimeter_mean"         

 [4] "area_mean"               "smoothness_mean"         "compactness_mean"       

 [7] "concavity_mean"          "concave.points_mean"     "symmetry_mean"          

[10] "fractal_dimension_mean"  "radius_se"               "texture_se"             

[13] "perimeter_se"            "area_se"                 "smoothness_se"          

[16] "compactness_se"          "concavity_se"            "concave.points_se"      

[19] "symmetry_se"             "fractal_dimension_se"    "radius_worst"           

[22] "texture_worst"           "perimeter_worst"         "area_worst"             

[25] "smoothness_worst"        "compactness_worst"       "concavity_worst"        

[28] "concave.points_worst"    "symmetry_worst"          "fractal_dimension_worst"

 

> # distribution and correlation check among input variables of WDBC

> install.packages("PerformanceAnalytics")

> library(PerformanceAnalytics)


> chart.Correlation(X[,c(1:10)], histogram=TRUE, col="grey10", pch=1) # MEAN




> chart.Correlation(X[,c(11:20)], histogram=TRUE, col="grey10", pch=1) # SE




> chart.Correlation(X[,c(21:30)], histogram=TRUE, col="grey10", pch=1) # WORST




GGally 패키지의 ggcorr() 함수를 사용하여 30개의 모든 설명변수간 상관계수를 구하고, -1~1의 상관계수에 따라 색을 달리하는 히트맵(heatmap)을 그려보았습니다.  역시나 여러개의 설명변수들이 서로 상관계수 0.9이상의 매우 강한 상관관계를 가지고 있음을 확인할 수 있습니다. 


> # heatmap plot of correlation coefficients of WDBC

> install.packages("GGally")

> library(GGally)

> ggcorr(X, name="corr", label=T)

 



설명변수들 간에 강한 상관성을 가지는 다중공선성(Multicolleniarity)가 존재하면 추정한 회귀계수의 분산이 매우 커지게 되어 추정한 회귀계수를 신뢰하기 힘들게 됩니다. 다시 말하면, 다중공선성이 있는 변수들을 사용해서 회귀계수를 추정하면, 원래 유의미하게 나와야 할 회귀계수가 검정을 해보면 유의미하지 않게 나올 수도 있으며, 반대로 유의미하지 않은 설명변수가 회귀계수 검정 결과 유의미하게 나오는 경우도 유발할 수 있습니다. 


게다가 설명변수들간 강한 상관관계가 존재할 경우 해석하는데도 문제가 생길 수 있습니다. 회귀모형의 경우 다른 설명변수를 고정(fix, control)한 상태에서 해당 설명변수가 한 단위 증가할 때 회귀계수 만큼 종속변수(목표변수)가 변화한다고 해석을 하는데요, 만약 다중공선성이 존재한다면 '다른 설명변수를 고정한다'는 설명변수의 독립성 가정이 안맞기 때문에 해석이 곤란해집니다. 따라서 다중공선성이 의심되면 처리를 해주는게 필요합니다. 


일반적으로 k개의 설명변수별 분산팽창지수(VIF: Variance Inflation Factor)를 구했을 때 가장 큰 VIF 값이 5 이상이면 다중공선성이 있다고 보며, VIF 값이 10 이상이며 다중공선성이 매우 심각하다고 평가합니다. 





다중공선성이 확인되면 이를 해결하기 위한 방법으로 (1) 상관계수가 가장 높은 변수를 제거(remove the highly correlated Xj variable), (2) 주성분분석(PCA), 요인분석(factor analysis), VAE(Variable Auto Encoder) 등의 알고리즘을 이용한 차원 축소 (dimension reduction), (3) 모수 추정 시 최소자승법(Least Squares Method) 대신에 Regularization (penalty)를 부여하는 Ridge regression, LASSO, PLS(Partial Least Squares Regression) 등을 사용하는 방법이 있습니다. 



지난번 포스팅에서 분석 방향 설정 부분에서 의사 선생님들이 진단하는데 사용하는 모델인 만큼 '모델 해석력(interpretability)'이 매우 중요할 것으로 예상이 된다고 했으므로, 이번 포스팅에서는 (1) 분산팽창지수(VIF) 10 이상인 설명변수를 제거하는 방법을 사용하겠습니다. 


R의 fmsb 패키지의 VIF() 함수를 사용해서 모든 설명변수의 분산팽창지수를 구한 후에 가장 큰 값이 10 이상이 경우 순차적으로 제거하겠습니다. 먼저, 예시로 첫번째와 두번째 설명변수인 radius_mean, texture_mean에 대한 분산팽창지수를 계산해본 것인데요, radius_mean 변수의 VIF가 3,806으로서 미친듯이 높음을 알 수 있습니다. (나머지 29의 설명변수와 매우 강한 선형상관관계를 가지고 있다는 뜻)



> # multicolleniarity check

> install.packages("fmsb")

> library(fmsb)

> VIF(lm(radius_mean ~ .,data=X))

[1] 3806.115

> VIF(lm(texture_mean ~ .,data=X))

[1] 11.88405

 



문제는 설명변수가 30개씩이나 되고, 30개 변수의 VIF를 모두 구해서 제일 큰 VIF 값의 설명변수를 제거하고, 다시 29개의 설명변수의 VIF를 모두 구해서 제일 큰 VIF 값의 설명변수를 제거하고.... 이런 단순 반복작업을 VIF 값이 10 이상인 설명변수가 없을 때 까지 반복해야만 한다는 점입니다. 


이처럼 반복작업인 경우 사용자 정의 함수(User Defined Function)을 짜놓고 사용하면 편리하고 시간도 줄일 수 있어서 좋습니다. 



[분산팽창지수를 구하고 VIF 10 이상인 변수 중 가장 큰 값을 순차적으로 제거하는 R 사용자 정의함수]

 

# Multi-colleniarity check and remove the highly correlated variables step by step

# UDF of stepwise VIF function with preallocated vectors

# code source: https://beckmw.wordpress.com/2013/02/05/collinearity-and-stepwise-vif-selection/


vif_func <- function(in_frame,thresh=10, trace=F,...){

  

  require(fmsb)

  

  if(class(in_frame) != 'data.frame') in_frame<-data.frame(in_frame)

  

  #get initial vif value for all comparisons of variables

  vif_init <- vector('list', length = ncol(in_frame))

  names(vif_init) <- names(in_frame)

  var_names <- names(in_frame)

  

  for(val in var_names){

    regressors <- var_names[-which(var_names == val)]

    form <- paste(regressors, collapse = '+')

    form_in <- formula(paste(val,' ~ .'))

    vif_init[[val]] <- VIF(lm(form_in,data=in_frame,...))

  }

  vif_max<-max(unlist(vif_init))

  

  if(vif_max < thresh){

    if(trace==T){ #print output of each iteration

      prmatrix(vif_init,collab=c('var','vif'),rowlab=rep('', times = nrow(vif_init) ),quote=F)

      cat('\n')

      cat(paste('All variables have VIF < ', thresh,', max VIF ',round(vif_max,2), sep=''),'\n\n')

    }

    return(names(in_frame))

  }

  else{

    

    in_dat<-in_frame

    

    #backwards selection of explanatory variables, stops when all VIF values are below 'thresh'

    while(vif_max >= thresh){

      

      vif_vals <- vector('list', length = ncol(in_dat))

      names(vif_vals) <- names(in_dat)

      var_names <- names(in_dat)

      

      for(val in var_names){

        regressors <- var_names[-which(var_names == val)]

        form <- paste(regressors, collapse = '+')

        form_in <- formula(paste(val,' ~ .'))

        vif_add <- VIF(lm(form_in,data=in_dat,...))

        vif_vals[[val]] <- vif_add

      }

      

      max_row <- which.max(vif_vals)

      #max_row <- which( as.vector(vif_vals) == max(as.vector(vif_vals)) )

      

      vif_max<-vif_vals[max_row]

      

      if(vif_max<thresh) break

      

      if(trace==T){ #print output of each iteration

        vif_vals <- do.call('rbind', vif_vals)

        vif_vals

        prmatrix(vif_vals,collab='vif',rowlab=row.names(vif_vals),quote=F)

        cat('\n')

        cat('removed: ', names(vif_max),unlist(vif_max),'\n\n')

        flush.console()

      }

      in_dat<-in_dat[,!names(in_dat) %in% names(vif_max)]

    }

    return(names(in_dat))

  }

}




위의 사용자정의함수 vif_func() 를 먼저 실행시키고, 다음으로 아래처럼 설명변수 X DataFrame과 VIF 기준(threshold)을 10으로 설정하고, 순차적인 제거(remove) 결과를 프린트하도록 해서 vif_func(X, thresh=10, trace=T) 함수를 실행시키면 아래와 같습니다. 



> # run vif_function

> X_independent <- vif_func(X, thresh=10, trace=T)

                                vif

radius_mean             3806.115296

texture_mean              11.884048

perimeter_mean          3786.400419

area_mean                347.878657

smoothness_mean            8.194282

compactness_mean          50.505168

concavity_mean            70.767720

concave.points_mean       60.041733

symmetry_mean              4.220656

fractal_dimension_mean    15.756977

radius_se                 75.462027

texture_se                 4.205423

perimeter_se              70.359695

area_se                   41.163091

smoothness_se              4.027923

compactness_se            15.366324

concavity_se              15.694833

concave.points_se         11.520796

symmetry_se                5.175426

fractal_dimension_se       9.717987

radius_worst             799.105946

texture_worst             18.569966

perimeter_worst          405.023336

area_worst               337.221924

smoothness_worst          10.923061

compactness_worst         36.982755

concavity_worst           31.970723

concave.points_worst      36.763714

symmetry_worst             9.520570

fractal_dimension_worst   18.861533


removed:  radius_mean 3806.115 


                               vif

texture_mean             11.882933

perimeter_mean          541.712931

area_mean               317.093211

smoothness_mean           7.990641

compactness_mean         38.106611

concavity_mean           65.978202

concave.points_mean      60.025840

symmetry_mean             4.203501

fractal_dimension_mean   15.677673

radius_se                75.101495

texture_se                4.185513

perimeter_se             67.720819

area_se                  41.089343

smoothness_se             4.017499

compactness_se           15.341790

concavity_se             15.234133

concave.points_se        11.399633

symmetry_se               5.175369

fractal_dimension_se      9.699518

radius_worst            616.350861

texture_worst            18.539292

perimeter_worst         375.408537

area_worst              304.471896

smoothness_worst         10.727206

compactness_worst        36.053404

concavity_worst          31.968539

concave.points_worst     36.763168

symmetry_worst            9.511243

fractal_dimension_worst  18.841897


removed:  radius_worst 616.3509 


                               vif

texture_mean             11.759131

perimeter_mean          325.641312

area_mean               237.012095

smoothness_mean           7.988003

compactness_mean         36.681620

concavity_mean           64.836935

concave.points_mean      60.019062

symmetry_mean             4.123603

fractal_dimension_mean   15.670406

radius_se                38.637579

texture_se                4.132025

perimeter_se             59.062677

area_se                  33.911923

smoothness_se             4.010296

compactness_se           15.304014

concavity_se             15.002055

concave.points_se        11.218541

symmetry_se               5.156085

fractal_dimension_se      9.542616

texture_worst            18.191599

perimeter_worst         308.052048

area_worst              168.343121

smoothness_worst         10.679641

compactness_worst        35.779767

concavity_worst          31.942417

concave.points_worst     35.761242

symmetry_worst            9.312564

fractal_dimension_worst  18.445566


removed:  perimeter_mean 325.6413 


                               vif

texture_mean             11.714252

area_mean                34.491349

smoothness_mean           7.964156

compactness_mean         31.979571

concavity_mean           64.655174

concave.points_mean      59.967015

symmetry_mean             4.123603

fractal_dimension_mean   14.921612

radius_se                36.056151

texture_se                4.092556

perimeter_se             42.980382

area_se                  32.570748

smoothness_se             3.914161

compactness_se           15.283194

concavity_se             14.769198

concave.points_se        10.464462

symmetry_se               5.128175

fractal_dimension_se      9.542575

texture_worst            18.112512

perimeter_worst         123.257811

area_worst               72.764912

smoothness_worst         10.648133

compactness_worst        34.263137

concavity_worst          31.681663

concave.points_worst     35.231031

symmetry_worst            9.268771

fractal_dimension_worst  18.287262


removed:  perimeter_worst 123.2578 


                              vif

texture_mean            11.679833

area_mean               28.534200

smoothness_mean          7.909212

compactness_mean        28.746302

concavity_mean          64.654796

concave.points_mean     59.816820

symmetry_mean            4.071436

fractal_dimension_mean  12.724264

radius_se               36.045576

texture_se               4.040107

perimeter_se            31.225949

area_se                 20.995394

smoothness_se            3.894739

compactness_se          15.199363

concavity_se            14.766025

concave.points_se       10.344938

symmetry_se              5.007681

fractal_dimension_se     9.302515

texture_worst           18.004692

area_worst              23.311066

smoothness_worst        10.619439

compactness_worst       34.253186

concavity_worst         31.669493

concave.points_worst    34.141124

symmetry_worst           9.077526

fractal_dimension_worst 18.285365


removed:  concavity_mean 64.6548 


                              vif

texture_mean            11.679763

area_mean               28.512892

smoothness_mean          7.543056

compactness_mean        26.110203

concave.points_mean     28.499831

symmetry_mean            4.064239

fractal_dimension_mean  12.668596

radius_se               35.617518

texture_se               4.034866

perimeter_se            31.178650

area_se                 19.985188

smoothness_se            3.872144

compactness_se          14.858964

concavity_se            11.587995

concave.points_se       10.013827

symmetry_se              5.005381

fractal_dimension_se     9.248401

texture_worst           18.003124

area_worst              23.026283

smoothness_worst        10.598388

compactness_worst       33.972256

concavity_worst         21.188494

concave.points_worst    32.115706

symmetry_worst           9.068073

fractal_dimension_worst 18.090939


removed:  radius_se 35.61752 


                              vif

texture_mean            11.465264

area_mean               25.427695

smoothness_mean          7.482881

compactness_mean        26.043226

concave.points_mean     27.801240

symmetry_mean            4.047744

fractal_dimension_mean  12.232750

texture_se               4.011486

perimeter_se            16.007813

area_se                 16.951023

smoothness_se            3.861783

compactness_se          14.577396

concavity_se            11.499061

concave.points_se        9.999462

symmetry_se              5.003384

fractal_dimension_se     8.800984

texture_worst           17.724662

area_worst              20.617761

smoothness_worst        10.595373

compactness_worst       33.960639

concavity_worst         21.056095

concave.points_worst    32.088040

symmetry_worst           9.065905

fractal_dimension_worst 18.084650


removed:  compactness_worst 33.96064 


                              vif

texture_mean            11.448852

area_mean               25.389376

smoothness_mean          7.479515

compactness_mean        19.208231

concave.points_mean     25.697888

symmetry_mean            3.982308

fractal_dimension_mean  10.961595

texture_se               3.998307

perimeter_se            15.960835

area_se                 16.935889

smoothness_se            3.859669

compactness_se           9.974677

concavity_se            10.850219

concave.points_se        9.805676

symmetry_se              4.941233

fractal_dimension_se     7.983689

texture_worst           17.701635

area_worst              20.613295

smoothness_worst        10.586935

concavity_worst         18.432076

concave.points_worst    30.596655

symmetry_worst           8.754474

fractal_dimension_worst 13.187594


removed:  concave.points_worst 30.59666 


                              vif

texture_mean            11.278555

area_mean               25.387830

smoothness_mean          7.162886

compactness_mean        19.175385

concave.points_mean     19.091402

symmetry_mean            3.918815

fractal_dimension_mean  10.902634

texture_se               3.937299

perimeter_se            15.808730

area_se                 16.917891

smoothness_se            3.637606

compactness_se           9.956527

concavity_se             9.775933

concave.points_se        5.347299

symmetry_se              4.900803

fractal_dimension_se     7.965699

texture_worst           17.427935

area_worst              20.406468

smoothness_worst         9.741783

concavity_worst         16.192763

symmetry_worst           8.435382

fractal_dimension_worst 13.047801


removed:  area_mean 25.38783 


                              vif

texture_mean            11.179202

smoothness_mean          7.162712

compactness_mean        18.843208

concave.points_mean     15.619381

symmetry_mean            3.895936

fractal_dimension_mean   9.707446

texture_se               3.937174

perimeter_se            15.619268

area_se                 16.655447

smoothness_se            3.632008

compactness_se           9.936443

concavity_se             9.705569

concave.points_se        5.250584

symmetry_se              4.872228

fractal_dimension_se     7.946733

texture_worst           17.236427

area_worst              10.626847

smoothness_worst         9.608528

concavity_worst         16.109962

symmetry_worst           8.409532

fractal_dimension_worst 13.023306


removed:  compactness_mean 18.84321 


                              vif

texture_mean            11.134313

smoothness_mean          6.970849

concave.points_mean     11.753066

symmetry_mean            3.829642

fractal_dimension_mean   7.907186

texture_se               3.890957

perimeter_se            15.333308

area_se                 16.345495

smoothness_se            3.552541

compactness_se           6.363339

concavity_se             9.367267

concave.points_se        5.245367

symmetry_se              4.871870

fractal_dimension_se     7.584276

texture_worst           17.232376

area_worst              10.602010

smoothness_worst         9.606389

concavity_worst         15.700019

symmetry_worst           8.401090

fractal_dimension_worst 13.023120


removed:  texture_worst 17.23238 


                              vif

texture_mean             1.715846

smoothness_mean          6.795612

concave.points_mean     11.715292

symmetry_mean            3.654734

fractal_dimension_mean   7.890069

texture_se               2.033874

perimeter_se            15.281161

area_se                 16.333806

smoothness_se            3.384881

compactness_se           6.337432

concavity_se             9.364521

concave.points_se        5.235966

symmetry_se              4.312472

fractal_dimension_se     7.575192

area_worst              10.540176

smoothness_worst         8.644833

concavity_worst         15.699140

symmetry_worst           7.294569

fractal_dimension_worst 13.021119


removed:  area_se 16.33381 


                              vif

texture_mean             1.709993

smoothness_mean          6.701262

concave.points_mean     11.653729

symmetry_mean            3.651771

fractal_dimension_mean   7.750052

texture_se               2.009042

perimeter_se             4.317808

smoothness_se            3.338723

compactness_se           6.317986

concavity_se             8.849322

concave.points_se        4.645375

symmetry_se              4.312339

fractal_dimension_se     7.575158

area_worst               8.677078

smoothness_worst         8.642994

concavity_worst         15.510661

symmetry_worst           7.265658

fractal_dimension_worst 12.938791


removed:  concavity_worst 15.51066 


> X_independent

 [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"

 



이렇게 VIF 10 이상인 설명변수를 순차적으로 제거해서 처음에 30개의 설명변수가 -> 17개의 설명변수로 줄어들었습니다. 




남은 17개 설명변수만을 가진 X_2 데이터프레임을 만들고, 상관계수 히트맵을 다시 그려보았습니다. 일부 변수가 상관계수 0.8인 경우가 있기는 합니다만, 0.9이상의 매우 강한 상관관계를 가진 설명변수는 없네요.  



> # explanatory/independant variables after VIF test

> X_2 <- X[, X_independent]

> str(X_2)

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

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

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

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

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

 $ fractal_dimension_mean : num  0.0787 0.0567 0.06 0.0974 0.0588 ...

 $ texture_se             : num  0.905 0.734 0.787 1.156 0.781 ...

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

 $ smoothness_se          : num  0.0064 0.00522 0.00615 0.00911 0.01149 ...

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

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

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

 $ symmetry_se            : num  0.03 0.0139 0.0225 0.0596 0.0176 ...

 $ fractal_dimension_se   : num  0.00619 0.00353 0.00457 0.00921 0.00511 ...

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

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

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

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

> # correlation heatmap plot again

> ggcorr(X_2, name="corr", label=T)



 




  (5) 설명변수 표준화, 척도 변환 (standardization, rescal)


지난번 첫번째 포스팅에서 데이터셋 설명할 때 설명변수들간의 측정 척도(scale)이 서로 다르다고 했는데요, 나중에 그래프 그려서 비교하기에 유용하도록 척도를 평균이 '0', 표준편차가 '1'인 새로운 척도로 표준화(standardization) 하도록 하겠습니다.  




scale() 함수를 이용하면 간단하게 표준화를 할 수 있습니다. summary() 함수로 확인해보니 평균이 '0'으로 중심이 바뀌었으며, max 값이 들쭉 날쭉 한걸로 봐서는 변수별로 outlier들이 꽤 있는 것처럼 보이네요. 로지스틱 회귀모형을 적합할 계획이므로 별도로 이상치(outlier) 처리는 다루지 않겠습니다. 



> # Standardization

> X_3 <- scale(X_2)

> summary(X_3)

  texture_mean     smoothness_mean    concave.points_mean symmetry_mean      fractal_dimension_mean

 Min.   :-2.2273   Min.   :-3.10935   Min.   :-1.2607     Min.   :-2.74171   Min.   :-1.8183       

 1st Qu.:-0.7253   1st Qu.:-0.71034   1st Qu.:-0.7373     1st Qu.:-0.70262   1st Qu.:-0.7220       

 Median :-0.1045   Median :-0.03486   Median :-0.3974     Median :-0.07156   Median :-0.1781       

 Mean   : 0.0000   Mean   : 0.00000   Mean   : 0.0000     Mean   : 0.00000   Mean   : 0.0000       

 3rd Qu.: 0.5837   3rd Qu.: 0.63564   3rd Qu.: 0.6464     3rd Qu.: 0.53031   3rd Qu.: 0.4706       

 Max.   : 4.6478   Max.   : 4.76672   Max.   : 3.9245     Max.   : 4.48081   Max.   : 4.9066       

   texture_se       perimeter_se     smoothness_se     compactness_se     concavity_se    

 Min.   :-1.5529   Min.   :-1.0431   Min.   :-1.7745   Min.   :-1.2970   Min.   :-1.0566  

 1st Qu.:-0.6942   1st Qu.:-0.6232   1st Qu.:-0.6235   1st Qu.:-0.6923   1st Qu.:-0.5567  

 Median :-0.1973   Median :-0.2864   Median :-0.2201   Median :-0.2808   Median :-0.1989  

 Mean   : 0.0000   Mean   : 0.0000   Mean   : 0.0000   Mean   : 0.0000   Mean   : 0.0000  

 3rd Qu.: 0.4661   3rd Qu.: 0.2428   3rd Qu.: 0.3680   3rd Qu.: 0.3893   3rd Qu.: 0.3365  

 Max.   : 6.6494   Max.   : 9.4537   Max.   : 8.0229   Max.   : 6.1381   Max.   :12.0621  

 concave.points_se  symmetry_se      fractal_dimension_se   area_worst      smoothness_worst 

 Min.   :-1.9118   Min.   :-1.5315   Min.   :-1.0960      Min.   :-1.2213   Min.   :-2.6803  

 1st Qu.:-0.6739   1st Qu.:-0.6511   1st Qu.:-0.5846      1st Qu.:-0.6416   1st Qu.:-0.6906  

 Median :-0.1404   Median :-0.2192   Median :-0.2297      Median :-0.3409   Median :-0.0468  

 Mean   : 0.0000   Mean   : 0.0000   Mean   : 0.0000      Mean   : 0.0000   Mean   : 0.0000  

 3rd Qu.: 0.4722   3rd Qu.: 0.3554   3rd Qu.: 0.2884      3rd Qu.: 0.3573   3rd Qu.: 0.5970  

 Max.   : 6.6438   Max.   : 7.0657   Max.   : 9.8429      Max.   : 5.9250   Max.   : 3.9519  

 symmetry_worst    fractal_dimension_worst

 Min.   :-2.1591   Min.   :-1.6004        

 1st Qu.:-0.6413   1st Qu.:-0.6913        

 Median :-0.1273   Median :-0.2163        

 Mean   : 0.0000   Mean   : 0.0000        

 3rd Qu.: 0.4497   3rd Qu.: 0.4504        

 Max.   : 6.0407   Max.   : 6.8408

 



이상으로 탐색적 데이터 분석 및 데이터 전처리 포스팅을 마치겠습니다. 

다음번 포스팅에서는 세번째로 '1차 변수 선택 및 목표변수와 설명변수 간 관계 분석'에 대해서 다루어보겠습니다. 


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



Posted by R Friend R_Friend

이번 포스팅부터 4번에 나누어서 로지스틱 회귀분석을 통한 유방암 예측(분류)를 소개하겠습니다. 


보통 통계학 수업 듣다보면 모델링 및 해석에 집중하고 탐색적 데이터 분석이나 데이터 전처리는 언급이 없거나 건너뛰는 경우가 많은데요, 이번 포스팅에서는 가급적이면 실제 데이터 분석을 하는데 있어서 거치게 되는 데이터 로딩, 탐색적 데이터 분석, 전처리, 시각화 등의 절차까지 모두 포함하여 상세하게 설명하려고 노력하겠습니다. 



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

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



  (1) WDBC(Wisconsin Diagnostic Breast Cancer) Dataset 소개


이번 분석에 사용할 데이터는 미국 중북부 위스콘신 주 남부 Madson에 있는 Wisconsin 대학병원 Dr.William H. Wolberg 교수님, 컴퓨터 공학과 W.Nick Street 교수님, Olvi L.Mangasarian 교수님이 1995년에 공개한 자료입니다. 

(for more information: https://archive.ics.uci.edu/ml/machine-learning-databases/breast-cancer-wisconsin/breast-cancer-wisconsin.names)


가는 주사바늘을 사용하여 낭종(주모니 혹은 내용물이나 멍울)의 세포나 조직을 떼어내서 조직검사한 이미지를 디지털로 변환하여 측정한 30개의 설명변수와, 환자 ID, 진단 결과(악성 M=malignant, 양성 B=benign) 의 변수를 포함하여, 총 32개의 변수로 구성되어 있고, 569명의 환자에 대해 조사한 데이터셋입니다. 


데이터셋이 공개되어 있는 주소는 https://archive.ics.uci.edu/ml/machine-learning-databases/breast-cancer-wisconsin/wdbc.data 입니다. 


* image source: https://pathology.jhu.edu/breast/my-results/types-of-breast-cancer/



외부 공개된 데이터셋을 data.table 패키지의 fread() 함수를 사용해서 불러온 후 R의 DataFrame 객체로 만들어보겠습니다. 



> rm(list=ls()) # clear all

> options(scipen=30)

> # data loading: WDBC (Wisconsin Diagnostic Breast Cancer)

> library(data.table)

data.table 1.11.8  Latest news: r-datatable.com

> library(curl)

> url <- c("https://archive.ics.uci.edu/ml/machine-learning-databases/breast-cancer-wisconsin/wdbc.data")

> wdbc <- data.frame(fread(url))

  % Total    % Received % Xferd  Average Speed   Time    Time     Time  Current

                                 Dload  Upload   Total   Spent    Left  Speed

  0     0    0     0    0     0      0      0 --:--:-- --:--:-- --:--:--     0  0  121k    0     0    0     0      0      0 --:--:-- --:--:-- --:--:--     0100  121k  100  121k    0     0  81060      0  0:00:01  0:00:01 --:--:-- 81060

> str(wdbc)

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

 $ V1 : int  842302 842517 84300903 84348301 84358402 843786 844359 84458202 844981 84501001 ...

 $ V2 : chr  "M" "M" "M" "M" ...

 $ V3 : num  18 20.6 19.7 11.4 20.3 ...

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

 $ V5 : num  122.8 132.9 130 77.6 135.1 ...

 $ V6 : num  1001 1326 1203 386 1297 ...

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

 $ V8 : num  0.2776 0.0786 0.1599 0.2839 0.1328 ...

 $ V9 : num  0.3001 0.0869 0.1974 0.2414 0.198 ...

 $ V10: num  0.1471 0.0702 0.1279 0.1052 0.1043 ...

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

 $ V12: num  0.0787 0.0567 0.06 0.0974 0.0588 ...

 $ V13: num  1.095 0.543 0.746 0.496 0.757 ...

 $ V14: num  0.905 0.734 0.787 1.156 0.781 ...

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

 $ V16: num  153.4 74.1 94 27.2 94.4 ...

 $ V17: num  0.0064 0.00522 0.00615 0.00911 0.01149 ...

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

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

 $ V20: num  0.0159 0.0134 0.0206 0.0187 0.0188 ...

 $ V21: num  0.03 0.0139 0.0225 0.0596 0.0176 ...

 $ V22: num  0.00619 0.00353 0.00457 0.00921 0.00511 ...

 $ V23: num  25.4 25 23.6 14.9 22.5 ...

 $ V24: num  17.3 23.4 25.5 26.5 16.7 ...

 $ V25: num  184.6 158.8 152.5 98.9 152.2 ...

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

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

 $ V28: num  0.666 0.187 0.424 0.866 0.205 ...

 $ V29: num  0.712 0.242 0.45 0.687 0.4 ...

 $ V30: num  0.265 0.186 0.243 0.258 0.163 ...

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

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




V1 은 환자 ID, V2는 진단 결과(악성 M, 양성 B), V3 ~ V32 까지 30개의 설명변수입니다. 설명변수로는 radius, texture, perimeter, area, smoothness, compactness, concavity, concave points, symmetry, fractal dimension 의 10개 지표에 대해서 mean(평균), se(standard error), worst(가장 나쁜 혹은 큰 측정치) 를 각각 측정(즉 10개 지표 * 3개 측정 기준 = 30개 설명변수)한 자료 입니다. 


(data info: https://archive.ics.uci.edu/ml/machine-learning-databases/breast-cancer-wisconsin/wdbc.names)




colnames() 함수로 변수명을 변경해보겠습니다. 


 

> # column names

> colnames(wdbc) <- c("id", "diagnosis", "radius_mean", "texture_mean", 

+                     "perimeter_mean", "area_mean", "smoothness_mean", "compactness_mean",

+                     "concavity_mean", "concave.points_mean", "symmetry_mean", "fractal_dimension_mean", 

+                     "radius_se", "texture_se", "perimeter_se", "area_se", 

+                     "smoothness_se", "compactness_se", "concavity_se", "concave.points_se", 

+                     "symmetry_se", "fractal_dimension_se", "radius_worst", "texture_worst", 

+                     "perimeter_worst", "area_worst", "smoothness_worst", "compactness_worst", 

+                     "concavity_worst", "concave.points_worst", "symmetry_worst", "fractal_dimension_worst")



str() 함수로 데이터셋 구조와 변수에 대해서 살펴보니 목표변수 diagnosis는 M(악성)과 B(양성) 의 두 개 범주(category)를 가진 범주형 데이터이며, 30개 설명변수는 모두 숫자형(numeric), 연속형 변수입니다. 각 설명변수별로 측정 단위가 다르므로 (길이, 면적 등) 나중에 표준화를 해서 분석하겠습니다. 



str(wdbc)

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

 $ id                     : int  842302 842517 84300903 84348301 84358402 843786 844359 84458202 844981 84501001 ...

 $ diagnosis              : chr  "M" "M" "M" "M" ...

 $ radius_mean            : num  18 20.6 19.7 11.4 20.3 ...

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

 $ perimeter_mean         : num  122.8 132.9 130 77.6 135.1 ...

 $ area_mean              : num  1001 1326 1203 386 1297 ...

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

 $ compactness_mean       : num  0.2776 0.0786 0.1599 0.2839 0.1328 ...

 $ concavity_mean         : num  0.3001 0.0869 0.1974 0.2414 0.198 ...

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

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

 $ fractal_dimension_mean : num  0.0787 0.0567 0.06 0.0974 0.0588 ...

 $ radius_se              : num  1.095 0.543 0.746 0.496 0.757 ...

 $ texture_se             : num  0.905 0.734 0.787 1.156 0.781 ...

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

 $ area_se                : num  153.4 74.1 94 27.2 94.4 ...

 $ smoothness_se          : num  0.0064 0.00522 0.00615 0.00911 0.01149 ...

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

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

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

 $ symmetry_se            : num  0.03 0.0139 0.0225 0.0596 0.0176 ...

 $ fractal_dimension_se   : num  0.00619 0.00353 0.00457 0.00921 0.00511 ...

 $ radius_worst           : num  25.4 25 23.6 14.9 22.5 ...

 $ texture_worst          : num  17.3 23.4 25.5 26.5 16.7 ...

 $ perimeter_worst        : num  184.6 158.8 152.5 98.9 152.2 ...

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

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

 $ compactness_worst      : num  0.666 0.187 0.424 0.866 0.205 ...

 $ concavity_worst        : num  0.712 0.242 0.45 0.687 0.4 ...

 $ concave.points_worst   : num  0.265 0.186 0.243 0.258 0.163 ...

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

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


> head(wdbc, 2)

      id diagnosis radius_mean texture_mean perimeter_mean area_mean smoothness_mean compactness_mean

1 842302         M       17.99        10.38          122.8      1001         0.11840          0.27760

2 842517         M       20.57        17.77          132.9      1326         0.08474          0.07864

  concavity_mean concave.points_mean symmetry_mean fractal_dimension_mean radius_se texture_se

1         0.3001             0.14710        0.2419                0.07871    1.0950     0.9053

2         0.0869             0.07017        0.1812                0.05667    0.5435     0.7339

  perimeter_se area_se smoothness_se compactness_se concavity_se concave.points_se symmetry_se

1        8.589  153.40      0.006399        0.04904      0.05373           0.01587     0.03003

2        3.398   74.08      0.005225        0.01308      0.01860           0.01340     0.01389

  fractal_dimension_se radius_worst texture_worst perimeter_worst area_worst smoothness_worst

1             0.006193        25.38         17.33           184.6       2019           0.1622

2             0.003532        24.99         23.41           158.8       1956           0.1238

  compactness_worst concavity_worst concave.points_worst symmetry_worst fractal_dimension_worst

1            0.6656          0.7119               0.2654         0.4601                 0.11890

2            0.1866          0.2416               0.1860         0.2750                 0.08902

 




  (2) 분석 목적 및 분석 방향 설정


WDBC 데이터셋의 diagnosis 진단 변수는 목표변수(반응변수, 종속변수, Output 변수)로서 2개의 class (악성 Malignant, 양성 Benign)를 가진 범주형 자료로서, 30개의 연속형 설명변수를 사용해서 2개의 class (악성, 양성 여부)를 분류(classification) 하는 문제입니다. 




이진 분류에 사용할 수 있는 분석 알고리즘은 무척 많습니다. 가령, KNN, Naive Bayes, Logistic Regression, SVM, Decision Tree, Random Forest, GBM, Deep Neural Network 등 많습니다. 


다만, 이번 데이터가 의료 데이터인지라 의사 선생님이 진단 결과에 대한 해석과 이해가 용이한 분석 모델이 활용도 측면에서 좋을 것으로 예상하여 로지스틱 회귀모형(Logistic Regression)을 선택하였습니다. 즉, 30개의 설명변수를 사용하여 유방암 악성(Malignant)에 속할 0~1 사이의 확률값을 계산하여 진단하는 분류/예측 모델을 만들어보겠습니다. (당연히 다른 분석 알고리즘도 얼마든지 사용 가능합니다 ^^)


참고로 일반화 선형회귀모형(GLM: Generalized Linear Model)은 반응변수가 정규분포를 따르는 연속형변수가 아닐 경우 & 고정효과 만을 다루는 경우에 사용하는데요, 특히 이항변수를 반응변수로 가지며 Logit link를 사용하는 GLM 이 바로 로지스틱 회귀모형이 되겠습니다. 

(참고: https://en.wikipedia.org/wiki/Generalized_linear_model)


[ 일반화 선형회귀모형 (GLM: Generalized Linear Model) 구분 ]





다음번 포스팅에서는 WDBC 데이터에 대한 탐색적 데이터 분석과 데이터 전처리에 대해서 다루겠습니다. 


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



Posted by R Friend R_Friend

이번 포스팅에서는 R로 쇼핑몰 웹사이트의 텍스트를 가져와서 최저가 Top 10 리스트를 선별해서 DataFrame으로 만든 후에, 이를 엑셀로 내보내기를 해보겠습니다. 


아래 코드는 제 블로그에 질문을 남겨주신 분께서 거의 다 짜셨구요, 저는 엑셀로 내보내는 부분의 에러를 바로 잡는 방법을 안내해드렸었습니다. 


textreadr, rvest, xlsx 등의 R 패키지와 사용법에 대해서 블로그의 본문에 정식으로 한번 더 소개해드리면 더욱 많은 분들께서 검색해서 이용하시기에 편리할 듯 해서 본문에 다시 한번 포스팅합니다. 



1. R 패키지 설치 및 로딩



# clear all

rm(list=ls())


# install and load libries

install.packages("textreadr")

install.packages("rvest")

install.packages("xlsx")


library(textreadr)

library(rvest)

library(xlsx)

 




2. 각 제품별로 엑셀 sheet 를 분리해서 최저가 Top 10 업체 결과 내보내기



urlbase <- "https://search.shopping.naver.com/search/all.nhn?query="

low <- c("에어팟","코카콜라","건전지","삼다수")


df <- data.frame()


# for loop statement

for (i in low){

  url <- paste(urlbase, i ,"&cat_id=&frm=NVSHATC",sep="")

  nv <- read_html(url, encoding = "UTF-8")

  st <- html_nodes(nv, '.mall_name') %>% html_text()

  st <- st[1:10]

  price3 <- html_nodes(nv, '._lowPriceByMall') %>% html_nodes('.price') %>% html_text()

  price3 <- gsub('최저가','',price3)

  price3 <- price3[1:10] # 최저가 Top 10

  df <- data.frame(가게명=st, 가격=price3)

  

  # append df to df_all

  df <- rbind(df, data.frame(가게명=url, 가격="r"))

  

  # adding sheetName

  write.xlsx(df, file="C:/Users/admin/Documents/df_sheet.xlsx", 

             sheetName=i, # excel sheet per each product or brand

             col.names=TRUE

             append=TRUE)

 

  # check progress

  print(paste0(i, " has been completed."))

}

 



위의 코드를 실행한 엑셀 파일은 아래와 같이 각 sheet 별로 제품이 구분이 되어있습니다. 






3. 모든 제품에 대해서 하나의 엑셀 sheet에 모두 모아서 최저가 Top 10 결과 내보내기



urlbase <- "https://search.shopping.naver.com/search/all.nhn?query="

low <- c("에어팟","코카콜라","건전지","삼다수")


df_all <- data.frame()


# for loop statement

for (i in low){

  url <- paste(urlbase, i ,"&cat_id=&frm=NVSHATC",sep="")

  nv <- read_html(url, encoding = "UTF-8")

  st <- html_nodes(nv, '.mall_name') %>% html_text()

  st <- st[1:10]

  price3 <- html_nodes(nv, '._lowPriceByMall') %>% html_nodes('.price') %>% html_text()

  price3 <- gsub('최저가','',price3)

  price3 <- price3[1:10]

  df <- data.frame(item=i, 가게명=st, 가격=price3)

  

  # append url info to df_all

  df <- rbind(df, data.frame(item=i, 가게명=url, 가격="r"))

  

  # combining df to df_all one by one

  df_all <- rbind(df_all, df)

  

  # check progress

  print(paste0(i, " has been completed."))

}


# exporting df_all to excel

write.xlsx(df_all, 

           file="C:/Users/admin/Documents/df_all.xlsx", 

           col.names=TRUE)

 




위의 코드를 실행한 엑셀 파일 결과는 아래와 같이 하나의 sheet에 각 제품(item) 변수로 구분이 되어서 한꺼번에 모아서 엑셀로 내보내기가 되었음을 알 수 있습니다. 





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

Posted by R Friend R_Friend

이번 포스팅에서는 카이제곱 적합도 검정(Chi-squared goodness of fit test)를 활용하여 관측치가 특정 분포를 따르는지 아닌지를 검정하는 방법을 소개하겠습니다.

 

다음과 같은 문제가 있다고 해봅시다.

 

 

[문제]  시간이 지남에 따라 일어난 어떤 특정 사건의 발생건수에 대한 결과가 다음과 같이 주어졌다고 하자. 그 때 특정 사건의 발생 회수가 포아송(λ) 분포를 따르는지에 대한 검정을 유의수준 5%에서 실시하시오.

 

 

 발생건수

0회

1회

2회 

3회 

4회 

5회 

6회 

7회 

8회 

9회 이상 

 관찰회수

6

20

36 

50 

52 

40 

25 

12 

 

 

(풀이)

특정 관측치들이 어떤 이론적 분포를 따르고 있는지를 검정하는 데에는 적합도 검정(Chi-squared goodness of fit test) 를 이용합니다.

 

검정해야 할 귀무가설(H0)와 대립가설(H1)을 세워보면,

  • 귀무가설 : 특정 사건의 발생 회수가 포아송 분포(Poisson(λ))를 따른다.
                     (즉, 관측값의 도수와 가정한 분포로부터의 기대관측도수가 같다)
  • 대립가설 : 특정 사건의 발생 회수가 포아송 분포를 따르지 않는다. (Not )

 

이며,

카이제곱 검정 통계량은  입니다. (이 문제에서는 i=0, 1, ..., 9)

이때 k 는 범주(class)의 개수를 나타내며, 는 도수분포의 각 구간에 있는 관측도수,

는 도수분포의 각 구간의 기대도수를 나타냅니다.

 

 

주어진 데이터로부터 가정된 포아송 분포의 평균(λ)의 추정치을 구하면, =3.844가 됩니다.

 

 = (0*6+1*20+2*36+3*50*4*52+5*40+6*25+7*12+8*4+9*5)/250 = 3.844

 

R 로 위의 관측치 데이터셋을 입력하고 가정된 포아송 분포의 평균 추정치를 구해보도록 하겠습니다.

 

 

> #===========================================
> # Chi-squared Goodness-of-fit Test
> # Poisson Distribution, 5% significant level
> #===========================================
> rm(list=ls()) # clear all
>
> # Class of Events
> classnum <- c(0:9)
>
> # Observed Frequency
> obsnum <- c(6, 20, 36, 50, 52, 40, 25, 12, 4, 5)
>
> # Sum of Observed Frequency
> obssum <- sum(obsnum); obssum
[1] 250
>
> # Calculate Lambda
> lamb <- sum(classnum*obsnum) / sum(obsnum); lamb
[1] 3.844 

 

 

 

포아송분포의 모수 λ에 대한 추정치로 3.844를 사용하여 각 class 별 가설적인 확률 (i=0, 1, ..., 9)를 구해보겠습니다.

 

에 i=0, 1, ..., 9 를 대입해서 풀어보면 아래의 표와 같이 가설적인 확률 를 구할 수 있습니다.

 

R 의 dpois() 함수를 사용하여 확률 구해봤는데요, sum(hyp_prob)=0.994로서 전체 합이 1이 안되고 0.994 이네요.

 

> # Calculate the hypothesised probabilities
> hyp_prob = round(dpois(classnum, lamb), 3); hyp_prob
 [1] 0.021 0.082 0.158 0.203 0.195 0.150 0.096 0.053 0.025 0.011
> sum(hyp_prob) # not 1
[1] 0.994

 

확률의 공리적 정의 상 표본공간의 P(S)=1 이어야 하므로 모든 class별 확률을 전부 더했을 때 1이 되도록 보정을 해줄 필요가 있습니다. 만약 보정을 안한 상태에서 바로 카이제곱 검정을 실시하면 아래와 같이 "확률값들의 총계는 반드시 1 이어야 합니다"와 같은 에러 메시지가 뜹니다.

 

> # Chi-squared goodness-of-fit test
> chisq.gof.test <- chisq.test(obsnum, p=hyp_prob) # error
Error in chisq.test(obsnum, p = hyp_prob) : 
  확률값들의 총계는 반드시 1 이어야 합니다 

 

 

R의 indexing을 사용하여 가장 마지막인 '9회 이상' class의 확률 를 조정하여 전체 class별 확률의 합이 1이 되도록 조정 해보겠습니다.

 

 

> # Change the last hypothesised probability to make sum(hyp_prob)=1
> hyp_prob[length(hyp_prob)] <- 1-sum(hyp_prob[1:(length(hyp_prob)-1)]); hyp_prob
[1] 0.021 0.082 0.158 0.203 0.195 0.150 0.096 0.053 0.025 0.017
> sum(hyp_prob) # 1
[1] 1 

 

 

 

발생확률 Pi를 구했으므로 전체 발생회수 합인 250과 각 class별 가정된 포아송분포로 부터의 확률을 곱한 250*로 부터 기대도수(Expected frequency)를 구할 수 있습니다.

 

 발생건수

0회

1회

2회 

3회 

4회 

5회 

6회 

7회 

8회 

9회 이상 

 계

 발생확률

0.021

0.082

0.158

0.203

0.195

0.150

0.096

0.053

0.025

0.017

 1

기대도수

(Ei=Pi*250)

5.25

20.50

39.50

50.75

48.75

37.50

24.00

13.25

6.25

4.25

250

 관측도수

(Oi)

6

20

36

50

52

40

25

12

4

5

250

 

 

이제 카이제곱 적합도 검정(Chi-squared goodness of fit test)을 할 준비가 되었네요. 가정된 포아송분포(λ=3.844)로 부터의 기대도수(Expected frequency)와 250개 샘플 데이터로 부터 실제 관측한 도수(Observed frequency) 간에 차이가 있다고 볼 수 있는 것인지, 관측치가 포아송분포를 따른다고 볼 수 있는 것인지 가설을 5% 유의수준 하에서 검정해보겠습니다.

 

R의 chisq.test() 함수를 사용하여 카이제곱 적합도 검정을 수행한 결과는 아래와 같습니다.

 

 

> # Chi-squared goodness-of-fit test
> chisq.gof.test <- chisq.test(obsnum, p=hyp_prob)
Warning message:
In chisq.test(obsnum, p = hyp_prob) :
카이제곱 approximation은 정확하지 않을수도 있습니다
> chisq.gof.test

Chi-squared test for given probabilities

data:  obsnum
X-squared = 1.9258, df = 9, p-value = 0.9926 

 

 

검정통계량인 X-squared 통계량이 1.9258이고 P-value 가 0.9926 이므로 유의수준 5% 하에서 관측값이 포아송 분포를 따른다는 귀무가설을 채택(Accept Null Hyperthesis)합니다.

 

 

R로 카이제곱 적합도 검정을 한 결과를 chisq.gof.test 라는 이름의 객체로 저장을 했는데요, names() 함수로 확인을 해보면 chisq.gof.test 객체 리스트 안에 여러 분석 결과 통계량들이 들어있음을 알 수 있습니다.

 

> names(chisq.gof.test)
[1] "statistic" "parameter" "p.value"   "method"    "data.name"
[6] "observed"  "expected"  "residuals" "stdres"  

 

 

이들 중에서 index을 해서 관측도수(Observed frequency)와 기대도수(Expected frequency)를 확인해보겠습니다.

 

 

> chisq.gof.test$observed
[1]  6 20 36 50 52 40 25 12  4  5
> chisq.gof.test$expected
[1]  5.25 20.50 39.50 50.75 48.75 37.50 24.00 13.25  6.25  4.25
> round(hyp_prob, 3)
[1] 0.021 0.082 0.158 0.203 0.195 0.150 0.096 0.053 0.025 0.017 

 

 

 

카이제곱 검정을 할 때 각 class의 기대도수가 너무 작을 경우 검정이 유효하지 않을 수가 있습니다. 위에 R분석 결과에서 빨간색으로 "카이제곱 approximation은 정확하지 않을 수도 있습니다"라는 warning message가 나온 이유가 '9회 이상' class의 기대도수가 4.25로서 작기 때문입니다. 최소 기대도수에 대해서는 일반적인 합의는 없으며, 보통 3, 4, 5개가 종종 사용되며, 3~5개 이하일 경우 인접 class와 병합해서 재분석을 하기도 합니다.

 

이번 문제에서는 최소 기대도수가 4.25이고, 또 P-value가 0.9926으로서 매우 높은 수준으로 귀무가설을 채택하였으므로 '9회 이상' class를 '8회'와 합쳐서 재분석하는 것은 필요없어 보입니다. (병합해서 재분석해도 결과는 동일하게 귀무가설 채택함)

 

마지막으로, 각 class별 관측도수와 기대도수를 그래프로 그려서 시각화로 비교를 해보겠습니다.

 

> # Plot of Poisson Distribution: Observed vs. Expected
> # Observed Frequency
> plot(c(0:9), chisq.gof.test$observed, 
+      main="Poisson Distribution: Observed vs. Expected", 
+      type='b', 
+      pch=0,
+      col='blue',
+      xlab="Number of Events",
+      ylab="Frequency",
+      ylim=c(0,55))
> 
> # Dual Y-axes plot (secondary Y-axis)
> par(new=T)
> 
> # Expected frequency
> plot(c(0:9), chisq.gof.test$expected, 
+      type='b', 
+      pch=1,
+      col='red', 
+      xlab="", 
+      ylab="",
+      ylim=c(0,55))
> 
> legend(x=6.5, y=50, 
+        c("Observed", "Expected"), 
+        pch=c(0,1), 
+        col=c('blue', 'red'))

 

 

 

 

 

Poisson(λ=3.844) 로부터의 각 class별 기대도수(Expected frequency)와 250개 샘플로 부터의 관측도수(Observed frequency)가 매우 유사한 분포를 띠고 있음을 눈으로도 확인할 수 있습니다.

 

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

 

Posted by R Friend R_Friend

이번 포스팅에서는 R ggplot2 로 y축이 2개 있는 이중 축 그래프 그리는 방법을 소개하겠습니다.

 

엑셀로는 이중 축 그래프 그리기가 그리 어렵지 않은데요, R의 ggplot 도 sec.axis 을 사용하면 두번째 y축을 추가할 수 가 있습니다.

 

먼저 예제로 사용할 데이터셋을 만들어보겠습니다.

 

 

#=========================================
# Dual y-axes plot using ggplot2 sec.axis
#=========================================

rm(list=ls()) # clear all

# make a DataFrame
data_val <- matrix(
  c(-4.2922960, 0.00000000000389755, 0.00000000268579,
  -3.9102123, 0.00000000000939471, 0.00000003082191,
  -3.5281286, 0.00000000002264455, 0.00000027610486,
  -3.1460449, 0.00000000005458189, 0.00000193070482,
  -2.7639612, 0.00000000013156254, 0.00001053865139,
  -2.3818776, 0.00000000031711433, 0.00004490363779,
  -1.9997939, 0.00000000076436268, 0.00014935003987,
  -1.6177102, 0.00000000184239690, 0.00038775418083,
  -1.2356265, 0.00000000444085801, 0.00078584150863,
  -0.8535428, 0.00000001070411038, 0.00124319933681,
  -0.4714591, 0.00000002580086522, 0.00153523161114,
  -0.0893754, 0.00000006218962756, 0.00147990680737,
  0.2927083, 0.00000014989999897, 0.00111358189390,
  0.6747920, 0.00000036131440584, 0.00065408965954,
  1.0568757, 0.00000087090114242, 0.00029990225777,
  1.4389594, 0.00000209919260108, 0.00010733701573,
  1.8210431, 0.00000505982314747, 0.00002998793957,
  2.2031268, 0.00001219600184343, 0.00000653989939,
  2.5852105, 0.00002939662278856, 0.00000111332724,
  3.5852105, 0.00029392734366929, 0.00000000000000,
  4.5852105, 0.00293538878457633, 0.00000000000000,
  5.5852105, 0.02896916461589227, 0.00000000000000,
  6.5852105, 0.25470155892952129, 0.00000000000000,
  7.5852105, 0.94711869936516890, 0.00000000000000,
  8.5852105, 0.99999999999982903, 0.00000000000000,
  9.5852105, 1.00000000000000000, 0.0000000000000),
  nrow=26,
  byrow=TRUE)

# convert a matrix into a dataframe
df <- data.frame(data_val)

# column name
colnames(df) <- c("x", "y1", "y2")

 

> head(df)
          x           y1           y2
1 -4.292296 3.897550e-12 2.685790e-09
2 -3.910212 9.394710e-12 3.082191e-08
3 -3.528129 2.264455e-11 2.761049e-07
4 -3.146045 5.458189e-11 1.930705e-06
5 -2.763961 1.315625e-10 1.053865e-05
6 -2.381878 3.171143e-10 4.490364e-05

 

 

 

두 개 y축의 요약 정보를 보니 최대값의 차이가 매우 크다는 것을 알 수 있습니다. y1 축과 y2 축의 비율을 보니 약 651배 차이가 나네요. 이런 상태에서 그래프를 그리면 y2 축의 그래프가 바닥에 쫙 붙어서 그려지므로 두 축 간의 패턴, 특징, 구조를 비교하기가 어렵게 됩니다. 따라서 y2 축의 데이터를 y1축과 비교하기 쉽도록 y2에 'y1과 y2의 최대값의 비율(max_ratio)'를 곱해서 그래프를 그려보겠습니다.

 

 

> # max ratio b/w 2 axes
> summary(df)
       x                 y1                  y2          
Min.   :-4.2923   Min.   :0.0000000   Min.   :0.000e+00 
1st Qu.:-1.9043   1st Qu.:0.0000000   1st Qu.:7.000e-10 
Median : 0.4838   Median :0.0000003   Median :8.539e-06 
Mean   : 1.1492   Mean   :0.1243873   Mean   :3.020e-04 
3rd Qu.: 3.3352   3rd Qu.:0.0002278   3rd Qu.:3.658e-04 
Max.   : 9.5852   Max.   :1.0000000   Max.   :1.535e-03 
> max_ratio <- max(df$y1)/max(df$y2); max_ratio
[1] 651.367

 

 

 

y1 데이터로는 선 그래프를 그리고, y2 데이터로는 막대그래프를 그려보겠습니다.

이때 y2 에 해당하는 두번째 축의 막대그래프를 그릴 때는

 - (1) max(y1)/max(y2) 로 계산한 두 축의 최대값의 비율인 651을 y2에 곱했으며

       (geom_bar(aes(y = y2*max_ratio) 

 - (2) scale_y_continuous(sec.axis = sec_axis(~.*maxratio, name="y2") 를 사용해서 오른쪽의 두번째 축을 추가하였습니다.

 

 

> # Dual y-axes plot using ggplot2 sec.axis
> library(ggplot2)
> 
> g <- ggplot(df, aes(x = x))
>   g <- g + geom_line(aes(y = y1), colour = "red", size = 2)
>   
>   # adding the relative result5.data.A, 
> # transformed to match roughly the range of the B
> g <- g + geom_bar(aes(y = y2*max_ratio),
> fill = "blue",
> stat = "identity")
> > # adding secondary axis > g <- g + scale_y_continuous(sec.axis = sec_axis(~.*max_ratio, name="y2*651")) > > g

 

 

 

해석하기에 편리하도록 데이터의 소스가 무엇인지를 나타내는 텍스트와 화살표를 추가해보겠습니다.

 

 

> # adding text > g <- g + annotate("text", x = -2, y = 0.75, colour="black", label="y2*651") > g <- g + annotate("text", x = 5, y = 0.5, colour="black", label="y1") > > # adding arrow > library(grid) > g <- g + annotate("segment", x = -1.8, y = 0.75, + xend = -1, yend = 0.6, colour="blue", arrow=arrow()) > g <- g + annotate("segment", x = 5.2, y = 0.5, + xend = 6.7, yend = 0.4, colour="red", arrow=arrow()) > g

 

 

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

 

Posted by R Friend R_Friend

이번 포스팅에서는 


(1) R Environment 상에 존재하는 여러개의 DataFrame 중에서 이름이 특정 조건을 만족하는 DataFrame을 선별하여 

 : ls(pattern = "xx")


(2) 하나의 List 로 묶는 방법, 

 : mget()


(3) List로 부터 특정 DataFrame을 Indexing 하는 방법

 : list[[1]], or  list[["name"]]


에 대해서 소개하겠습니다. 



먼저, 이번 포스팅의 예제의 특성 상 Environment 와 Console을 깨끗하게 청소하는 것부터 시작하겠습니다.  rm(list=ls()) 로 Environment에 있는 모든 객체를 삭제하고, cat("\014")로 Console을 백지상태로 만들어주었습니다. 



# To make an environment, console clear 

rm(list=ls()) # clearing environment

cat("\014") # clearing Console




객체 이름에 'data_'를 공통으로 포함한 3개의 DataFrame 예제와, 객체 이름에 'file_'을 포함한 3개의 DataFrame 예제를 만들어보겠습니다. 



# To make several sample DataFrames

data_1 <- data.frame(var1 = c(1, 2), var2 = c(3, 4))

data_2 <- data.frame(var1 = c('a', 'b'), var2 = c('c', 'd'))

data_3 <- data.frame(var1 = c(TRUE, TRUE), var2 = c(FALSE, TRUE))

file_1 <- data.frame(var1 = c(1, 2), var2 = c(3, 4))

file_2 <- data.frame(var1 = c('a', 'b'), var2 = c('c', 'd'))

file_3 <- data.frame(var1 = c(TRUE, TRUE), var2 = c(FALSE, TRUE))

  [DataFrame samples]

 



Environment에 생성된 객체를 확인해보려면 ls() 함수를 사용하면 됩니다. 그리고 특정 문자열을 포함한 객체만을 선별해서 찾아보려면 ls(pattern = "xx") 처럼 pattern = "xx" 를 추가하면 됩니다. 



> # To list up all objects in an environment

> ls()

[1] "data_1" "data_2" "data_3" "file_1" "file_2" "file_3"

> # To list up objects which have a matching 'pattern' in a DataFrame name

> ls(pattern = "data_")

[1] "data_1" "data_2" "data_3"

> ls(pattern = "file_")

[1] "file_1" "file_2" "file_3"

 




다음으로, mget() 함수를 사용하여 ls(pattern = "data_"), ls(pattern = "file_"로 각각 선별한 DataFrame 객체들을 2개의 List로 각각 묶어서 생성해보겠습니다.  



> # To combine all DataFrame into a List using mget()

> list_df_data <- mget(ls(pattern = "data_"))

> list_df_data

$data_1

  var1 var2

1    1    3

2    2    4


$data_2

  var1 var2

1    a    c

2    b    d


$data_3

  var1  var2

1 TRUE FALSE

2 TRUE  TRUE


> list_df_file <- mget(ls(pattern = "file_"))

> list_df_file

$file_1

  var1 var2

1    1    3

2    2    4


$file_2

  var1 var2

1    a    c

2    b    d


$file_3

  var1  var2

1 TRUE FALSE

2 TRUE  TRUE






마지막으로, List에 묶인 DataFrame을 [[ ]] 을 사용하여 위치(숫자) 혹은 이름으로 Indexing 하는 방법을 아래에 소개합니다. 



> # To index one of object in a List using [[ ]]

> list_df_data[[1]] # using 'number'

  var1 var2

1    1    3

2    2    4

> list_df_data["data_1"] # using 'name'

$data_1

  var1 var2

1    1    3

2    2    4 



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

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



Posted by R Friend R_Friend

이번 포스팅은 페이스북의 R User Group 에서 아이디어 토론 주제로 올라왔길레 한번 짜봤던 코드입니다.  


'0'과 '1'로 가지는 긴 벡터를 처음부터 끝까지 쭉 훓어가면서 


(1) '0'이 연속으로 3번 나오면 그 구간을 기준으로 나누어서 

=> (2) 나누어진 구간을 새로운 벡터로 생성하기


입니다.  


'0'이 연속으로 나오는 회수는 분석가가 필요로 하는 회수로 지정할 수 있도록 매개변수(argument)로 지정해서 프로그래밍해보겠습니다. 


간단한 예제 벡터로서, '0'과 '1'을 원소로 해서 30개의 무작위수(이항분포 랜덤 샘플링, 0이 80%, 1이 20%)로 구성된 벡터를 생성해보겠습니다. 


 

> rm(list=ls()) # clean all


# Sample vector


> set.seed(123) # for reproducibility

> vec_raw <- sample(c(0, 1), size=30, replace=TRUE, prob=(c(0.8, 0.2)))

> vec_raw

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

> vec <- vec_raw # copy for manipulation




R 프로그래밍을 하고자 하는 일의 아웃풋 이미지는 아래와 같습니다. 

'0'이 연속 세번 나오는 구간을 노란색으로 표시를 했습니다. 이 구간을 기준으로 원래의 벡터를 구분해서 나눈 후에(split), 각각을 새로운 벡터로 생성하여 분리해주는 프로그램을 짜보는 것이 이번 포스팅의 주제입니다. 




아래에 wihle() 반복문, if(), else if() 조건문, assign() 함수를 사용해서 위의 요건을 수행하는 코드를 짜보았습니다. 원래 벡터에서 제일 왼쪽에 있는 숫자가 '0'이 아니면 하나씩 빼서 빈 벡터에 차곡차곡 쌓아 놓구요, '0'이 3개 나란히 있는 '000'이 나타나면 그동안 쌓아두었던 벡터를 잘라내서 다른 이름의 벡터로 저장하고, '000'도 그 다음 이름의 벡터로 저장하면서 원래 벡터에서 '000'을 빼놓도록 하는 반복문입니다. assign() 함수는 저장하는 객체의 이름에 paste()를 사용해서 변화하는 숫자를 붙여주고자 할 때 사용하는 함수인데요, 자세한 사용법 예제는 여기(=> http://rfriend.tistory.com/108)를 참고하세요. 


아래 프로그램에서 'bin_range'에 할당하는 숫자를 변경하면 원하는 회수만큼 '0'이 반복되었을 때 구간을 분할하도록 조정할 수 있습니다.  가령, '0'이 연속 10번 나오면 분할하고 싶으면 bin_range <- 10  이라고 할당해주면 됩니다. 



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

# Vector Bin Split by successive '0's criteria

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


# Setting

vec_tmp <- c() # null vector

bin_range <- 3 # successive '0's criterion

vec_idx <- 1

vec_num <- 1


# Zero_Bin_Splitter

while (length(vec) > 0){

  if (sum(vec[1:bin_range], na.rm=T) != 0){

    vec_tmp[vec_idx] <- vec[1]

    vec <- vec[-1]

    vec_idx <- vec_idx + 1

  } else if (is.null(vec_tmp)){

    assign(paste0("vec_split_", vec_num), vec[1:bin_range])

    vec <- vec[-(1:bin_range)]

    vec_idx <- 1 # initialization

    vec_num <- vec_num + 1

  } else {

    assign(paste0("vec_split_", vec_num), vec_tmp)

    vec_tmp <- c() # initialization

    vec_num <- vec_num + 1

    assign(paste0("vec_split_", vec_num), vec[1:bin_range])

    vec <- vec[-(1:bin_range)]

    vec_idx <- 1 # initialization

    vec_num <- vec_num + 1

  }

}


# delete temp vector and arguments

rm(vec, vec_tmp, bin_range, vec_idx, vec_num)




아래는 위의 프로그램을 실행시켰을 때의 결과입니다.  원래 의도했던대로 잘 수행되었네요. 




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

아래의 코드는 페이스북 R User Group의 회원이신 June Young Lee 님께서 data.table 라이브러리를 사용해서 짜신 것입니다.  코드도 간결할 뿐만 아니라, 위에 제가 짠 코드보다 월등히 빠른 실행속도를 자랑합니다. 대용량 데이터를 빠른 속도로 조작, 처리하려면 data.table 만한게 없지요. 아래 코드는 저도 공부하려고 옮겨 놓았습니다. 


똑같은 일을 하는데 있어서도 누가 어떤 로직으로 무슨 패키지, 함수를 사용해서 짜느냐에 따라서 복잡도나 연산속도가 이렇게 크게 차이가 날 수 있구나 하는 좋은 예제가 될 것 같습니다. 이런게 프로그래밍의 묘미이겠지요? ^^  June Young Lee 님께 엄지척! ^^b



# by June Young Lee


library(data.table)


# 0이 세번 연속으로 들어 있는 아무 벡터 생성

vec <- c(0,0,0,3,10,2,3,0,4,0,0,0,1,2,50,4,0,0,32,1,0,0,0,1,1)


# data.table함수인 rleid이용하여 연속인 행들을 구분하는 idx1 컬럼 생성 

dt <- data.table(v1=vec, idx1=rleid(vec))


# idx1 기준으로 갯수를 세어 N 컬럼 생성 

dt[,N:=.N, by=idx1]


# 0이 세번 연속으로 나온 그룹(=N컬럼변수가 3인)을 idx2로 체크하기 

dt[v1==0&N==3,idx2:=1L]


# split 하기 위한 기준 벡터를 마찬가지로 rleid함수를 이용하여 idx3로 생성


dt[,idx3:=rleid(idx2)]


# data.table의 split 함수 적용하고, lapply 적용하여 v1 칼럼만 추출 

res <- lapply(split(dt, by="idx3"), function(x) x[,v1])

res

# $`1`

# [1] 0 0 0

# $`2`

# [1] 3 10 2 3 0 4

# $`3`

# [1] 0 0 0

# $`4`

# [1] 1 2 50 4 0 0 32 1

# $`5`

# [1] 0 0 0

# $`6`

# [1] 1 1


# 0이 3번 연속으로 나온 놈들을 남겨두도록 작성한 이유는, 

# 위치확인 & 확인용입니다. 

# 필요없으면, 아래 정도의 코드를 추가하면 되겠네요.


res2 <- res[!unlist(lapply(res, function(x) length(x)==3&sum(x)==0))]

res2

# $`2`

# [1] 3 10 2 3 0 4

# $`4`

# [1] 1 2 50 4 0 0 32 1

# $`6`

# [1] 1 1

 


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


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



Posted by R Friend R_Friend

페이스북의 KRUG (Korean R User Group) 에서 2017.11.18일 주말 퀴즈로 아래의 시각화를 푸는 문제가 있어서 재미로 풀어보았습니다. 


처음엔 금방 코드 짤 것으로 예상했는데요, 출제자분께서 중간 중간 장애물을 심어놓으셔서 제 예상보다 좀 더 걸렸네요. 


문제 풀면 KRUG 운영자분께서 스타벅스 기프티콘 주신다고 하는데요, 기대되네요. ^___^



# KRUG's plot quiz, as of 18th NOV. 2017


library(ggplot2)

#install.packages("dplyr")

library(dplyr)

str(mtcars)


# making 'model' variable

mtcars$model <- rownames(mtcars)


# ordering by cyl, hp

mtcars_ord <- arrange(mtcars[,c('model', 'cyl', 'hp')], cyl, desc(hp))


# marking an inflection point within cluster : threshold >= 40

mtcars_cyl <- c(4, 6, 8)

mtcars_ord_2 <- data.frame()


for (i in 1:length(mtcars_cyl)) {

  

  mtcars_tmp <- subset(mtcars_ord, subset = (cyl == mtcars_cyl[i]))

  

  for (j in 1:nrow(mtcars_tmp)) {

    if (j != nrow(mtcars_tmp) & mtcars_tmp$hp[j] - mtcars_tmp$hp[j+1] >= 40) {

      mtcars_tmp$hp_outlier[j] = '1_outlier'

    } else {

      mtcars_tmp$hp_outlier[j] = '2_normal'

    }

  }

  

  mtcars_ord_2 <- rbind(mtcars_ord_2, mtcars_tmp)

  rm(mtcars_tmp)

}


# converting cyl variable type from numeric to factor

mtcars_ord_2$cyl_cd <- paste0("cyl :", mtcars_ord_2$cyl)


model_order <- mtcars_ord_2$model[order(mtcars_ord_2$cyl_cd, 

                                        mtcars_ord_2$hp, 

                                        decreasing = FALSE)]


mtcars_ord_2$model <- factor(mtcars_ord_2$model, levels = model_order)



# drawing cleveland dot plot

ggplot(mtcars_ord_2, aes(x = hp, y = model)) +

  geom_point(size = 2, aes(colour = hp_outlier)) +

  scale_colour_manual(values = c("red", "black")) + 

  theme_bw() +

  facet_grid(. ~ cyl_cd, scales = "free_y", space = "free_y") +

  xlim(0, max(mtcars_ord_2$hp)) +

  geom_hline(yintercept = nrow(mtcars_ord_2[mtcars_ord_2$cyl == 4,]) + 0.5, 

             colour = "black", linetype = "dashed", size = 0.5) +

  geom_hline(yintercept = nrow(mtcars_ord_2) - 

               nrow(mtcars_ord_2[mtcars_ord_2$cyl == 8,]) + 0.5, 

             colour = "black", linetype = "dashed", size = 0.5) +

  theme(legend.position = 'none')

 




R Korea - KRUG(Korean R User Group) 에 문제 출제해주신 정우준님께서 나중에 정답지 올려주신 코드도 아래에 같이 공유합니다. 제가 짠 코드보다 한결 간결하네요. 



library(data.table)
library(dplyr)
library(ggplot2)

mtcars %>%
mutate(car.name=rownames(.)) %>%
arrange(cyl, hp) %>%
mutate(order.key=1:n()) -> data

data %>%
ggplot(aes(x=hp, y=reorder(car.name, order.key))) +
geom_point(
colour=case_when(
data$car.name %in% c('Ferrari Dino','Maserati Bora') ~ 'red', 
TRUE ~ 'black')) +
geom_hline(yintercept = 11.5, linetype='dashed') +
geom_hline(yintercept = 18.5, linetype='dashed') +
facet_wrap(~ cyl, labeller = label_both) +
scale_x_continuous(limits=c(0,max(data$hp))) +
theme_bw() +
theme(axis.title.y=element_blank())

 






Posted by R Friend R_Friend

문자열이나 숫자를 특정 형식으로 길이를 지정해주면 데이터를 출력했을 때 깨끗하게 정리가 되어 보이기 때문에 가독성이 좋아집니다. 


혹은 데이터가 특정 형식(format)으로 DB에 이미 지정이 되어 있어서 데이터 간 병합이나 join을 하기 위해 특정 형식으로 데이터를 표준화 해주어야 할 경우가 있습니다. 


이번 포스팅에서는 {base} package의 sprintf() 함수를 사용해서 


 - (1) 문자열을 매개변수 width 길이로 만들고, 빈 자리는 '0'으로 채우기 : sprintf("%05d", var)


 - (2) 소수점 숫자(numeric)의 자리수를 지정해주기 : sprintf(".5f", var)


하는 방법에 대해서 알아보겠습니다. 


이번 포스팅의 함수 sprintf()는 데이터 전처리할 때 종종 사용하는 편이예요. 



 (1) 문자열을 특정 길이로 만들고, 빈 자리수만큼 '0'을 채우기 : sprintf("%05d", var)





1자리, 2자리, 3자리, 4자리를 가진 데이터를 가지고 예제로 사용할 간단한 DataFrame을 만들어보겠습니다. 



> # making a sample DataFrame

> df <- data.frame(var1 = c(1, 11, 111, 1111))

> df

  var1

1    1

2   11

3  111

4 1111

 




위의 예제 데이터셋을 가지고, 칼럼 var1의 데이터를 


- '1 자리수를 가진 문자열'로 만들되, '1자리수가 안되면 모자라는 자리수 만큼'0'으로 채우기'

- '2 자리수를 가진 문자열'로 만들되, '2자리수가 안되면 모자라는 자리수 만큼'0'으로 채우기'

- '3 자리수를 가진 문자열'로 만들되, '3자리수가 안되면 모자라는 자리수 만큼'0'으로 채우기'

- '4 자리수를 가진 문자열'로 만들되, '4자리수가 안되면 모자라는 자리수 만큼'0'으로 채우기'

- '5 자리수를 가진 문자열'로 만들되, '5자리수가 안되면 모자라는 자리수 만큼'0'으로 채우기'


를 해보겠습니다. 


만약 매개변수 자리수보다 데이터의 길이가 더 크다면 '0'이 채워지지는 않습니다.  아래의 예제의 결과를 View(df)로 해서 보면 원래의 변수 var1과 sprintf() 함수를 사용해서 만든 var1_01d 변수의 데이터 출력 형식이 다른 것을 알 수 있습니다. 


그리고 class 함수로 데이터 형식을 살펴보니 원래 변수 var1은 숫자형(numeric)이지만 sprintf() 함수로 만든 새로운 변수는 요인형(factor)의 문자열로 바뀌어 있음을 알 수 있습니다. 



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

> # (1) sprintf(%03d, var) : Format number as fixed width, with leading zeros

> df <- transform(df, 

+                 var1_01d = sprintf("%01d", var1), 

+                 var1_02d = sprintf("%02d", var1), 

+                 var1_03d = sprintf("%03d", var1), 

+                 var1_04d = sprintf("%04d", var1), 

+                 var1_05d = sprintf("%05d", var1))

> df

  var1 var1_01d var1_02d var1_03d var1_04d var1_05d

1    1        1       01      001     0001    00001

2   11       11       11      011     0011    00011

3  111      111      111      111     0111    00111

4 1111     1111     1111     1111     1111    01111



> View(df)



> sapply(df, class)

     var1  var1_01d  var1_02d  var1_03d  var1_04d  var1_05d

"numeric"  "factor"  "factor"  "factor"  "factor"  "factor"




 (2) 소수점 숫자(numeric)의 자리수를 지정해주기 : sprintf(".5f", var)


무리수인 자연상수 e의 소수점 10째 자리까지의 수를 대상으로 sprintf("%.5f", e) 함수를 사용해서 소수점의 자리수를 설정해보겠습니다. "%.숫자f"의 숫자 만큼 소수점을 표시해주는데요, 반올림을 해서 표시해줍니다. 아래의 예제를 보시면 금방 이해할 수 있을 것입니다. 



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

> # (3) sprintf("%.5f", x) : formatting decimal point, 

> e <- c(2.7182818284) # mathematical constant, the base of the natural logarithm

> sprintf("%.0f", e)

[1] "3"

> sprintf("%.1f", e)

[1] "2.7"

> sprintf("%.2f", e)

[1] "2.72"

> sprintf("%.3f", e)

[1] "2.718"

> sprintf("%.5f", e)

[1] "2.71828"

> sprintf("%.10f", e)

[1] "2.7182818284"

 




아래의 예시는 sprintf("%숫자.f, e)로 '숫자' 부분에 매개변수로 정수 부분의 자리수를 지정해주는 예시입니다. 소수점의 자리도 모두 포함해서 '숫자' 부분 매개변수만큼의 길이로 표시 형식을 맞추어줍니다. 



> e <- c(2.7182818284) # mathematical constant, the base of the natural logarithm

> sprintf("%1.1f", e)

[1] "2.7"

> sprintf("%2.1f", e)

[1] "2.7"

> sprintf("%3.1f", e)

[1] "2.7"

> sprintf("%5.1f", e)

[1] "  2.7"

> sprintf("%10.1f", e)

[1] "       2.7"

>  



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

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





Posted by R Friend R_Friend

이번 포스팅에서는 R ggplot2 패키지의 coord_fixed(ratio = number) 옵션을 사용해서 그래프의 크기 (가로, 세로 비율)를 조정하는 방법을 소개하겠습니다. 


1 ~ 4의 정수 좌표를 가지는 X 축과 

1 ~ 4의 정수 좌표를 가지는 Y 축을 가지는 

가상의 데이터를 사용해서 예를 들어보겠습니다. 


X축이 Y축이 1:1 비율인 그래프가 되겠습니다. 



[ 가상 데이터셋 생성 ]



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

# Rssizing the plot

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


install.packages("ggplot2")

library(ggplot2)


# X, Y coordinate and segments(A, B category) data set

my.df <- data.frame(XCoord = c(1, 1, 1, 2, 2, 2, 3, 3, 4, 4), 

                    YCoord = c(1, 2, 4, 1, 3, 4, 2, 4, 1, 4), 

                    Seg = c("A", "A", "A", "A", "B", "A", "B", "B", "B", "B"))





X와 Y좌표별로 Seg. 변수의 "A"는 검정색, "B"는 흰색으로 색깔을 지정해서 Heatmap을 그려보겠습니다. 

X좌표와 Y좌표가 1~4 범위를 동일하게 가지고 있으므로 크기에 대한 설정없이 디폴트 세팅으로 그래프를 그리면 아래 처럼 정사각형 1:1 비율로 그래프가 그려집니다. 



[ 그림 1 ] 원본 그래프 (Original Plot) : 가로축과 세로축이 1:1



# Original plot

ggplot(my.df, aes(x=XCoord, y=YCoord, fill=Seg)) +

  geom_tile(colour="gray80") +

  scale_fill_manual(values = c("black", "white")) +

  ggtitle("Original Heatmap by X4*Y4 size")






[ 그림 2 ] 가로축과 세로축의 비율을 1:2 로 설정하기 : coord_fixed(ratio = 2)


[그림 1]에서 원본 이미지가 가로축과 세로축이 1:1 비율의 정사각형 그래프였는데요, 이를 가로:세로 비율을 1:2로 세로가 가로의 2배 비율인 그래프(가로:세로 = 1:2)로 바꾸어 주려면 coord_fixed(ratio = 2) 를 설정해주면 됩니다. 



# Resized plot using coord_fixed(ratio = number)

ggplot(my.df, aes(x=XCoord, y=YCoord, fill=Seg)) +

  geom_tile(colour="gray80") +

  scale_fill_manual(values = c("black", "white")) +

  coord_fixed(ratio = 2) +

  ggtitle("Heatmap : Resized X:Y axis size with 1:2 ratio")







[ 그림 3 ] 가로축과 세로축의 비율을 2:1 로 설정하기 : coord_fixed(ratio = 0.5)


가로와 세로축 비율이 1:1인 원본 이미지 [그림 1] 을 가로가 세로축의 2배인 그래프 (가로:세로 = 2:1)로 바꾸고 싶다면 coord_fixed(ratio = 0.5) 로 설정해주면 됩니다. 



# Resized plot using coord_fixed(ratio = number)

ggplot(my.df, aes(x=XCoord, y=YCoord, fill=Seg)) +

  geom_tile(colour="gray80") +

  scale_fill_manual(values = c("black", "white")) +

  coord_fixed(ratio = 0.5) +

  ggtitle("Heatmap : Resized X:Y axis size with 2:1 ratio")









아래는 EBImage 패키지의 resize() 함수를 사용해서 png 이미지 파일로 출력할 때 이미지 크기를 (a) 특정 가로, 세로 크기로 설정해주는 방법과 (b) 비율로 설정해주는 방법입니다. 

R code는 stackoverflow 의 답변 중에서 aoles 님께서 달아놓은 것인데요, 코드 그대로 인용해서 소개합니다. ( * R code 출처 : https://stackoverflow.com/questions/35786744/resizing-image-in-r )


#==============
# Image resize using EBImage package

# installing EBImage package
source("http://bioconductor.org/biocLite.R")
biocLite("EBImage")


# resizing image using EBImage package's resize() function

library("EBImage")

x <- readImage(system.file("images", "sample-color.png", package="EBImage"))


# width and height of the original image
dim(x)[1:2]


# scale to a specific width and height
y <- resize(x, w = 200, h = 100)


# scale by 50%; the height is determined automatically so that
# the aspect ratio is preserved
y <- resize(x, dim(x)[1]/2)


# show the scaled image
display(y)


# extract the pixel array
z <- imageData(y)


# or
z <- as.array(y)

 



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

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



Posted by R Friend R_Friend