자기상관계수(Autocorrelation coefficients), 자기상관그림(Autocorrelation Plot) (a) 데이터가 시간에 의존하는 것 없이 무작위성 (randomness, no time dependence in the data) 을 띠는지 확인하는데 사용됩니다. 그리고 (b) 시계열분석에서 Box-Jenkins ARIMA (Autoregressive Integrated Moving Average) 모형 식별 단계(model identification stage)에서 AR 차수를 구하는데도 사용됩니다. 


이번 포스팅에서는 R을 사용해서 


(1) 자기상관계수(Autocorrelation coefficients)를 계산하고, 

(2) 자기상관계수의 95% 신뢰구간(95% confidence interval)을 계산하고, 

(3) 자기상관그림(Autocorrelation Plot)을 그리는 방법


을 소개하겠습니다. 


본론에 들어가기 전에 먼저, 상관계수(correlation coefficients)와 자기상관계수(autocorrelation coefficients)를 비교해보겠습니다. 


상관계수와 자기상관계수의 유사한 점은 두 연속형 변수 간의 관계를 분석하다는 점입니다. 기본 개념은 공분산을 표준편차로 나누어서 표준화해주며, (자기)상관계수가 -1 ~ 1 사이의 값을 가지게 됩니다. 


상관계수와 자기상관계수가 서로 다른 점은, 상관계수는 특정 동일 시점을 횡단면으로 해서 Y와 다른 X1, X2, ... 변수들 간의 관계를 분석합니다. 반면에 자기상관계수는 동일한 변수(Yt, Yt-1, Yt-2, ...)의 서로 다른 시간 차이 (time lag) 를 두고 관계를 분석하는 것입니다. 


기존에 cross-sectional 관점의 Y와 X 변수들 간의 상관관계 분석에 많이 익숙해져 있는 경우에, 시계열 분석을 공부할 때 보면 자기 자신의 시간 차이에 따른 자기상관관계 관점이 처음엔 헷갈리고 잘 이해가 안가기도 합니다. 관점이 바뀐것일 뿐 어려운 개념은 아니므로 정확하게 이해하고 가면 좋겠습니다. 




  (1) 자기상관계수(Autocorrelation coefficients) 구하고, 자기상관그림 그리기


자기상관계수는 아래의 공식과 같이  로 계산합니다. 



간단한 예로서, 1 ~ 50까지의 시간 t에 대해서 싸인곡선 형태의 주기적인 파동을 띠는 값에 정규확률분포 N(0, 1) 에서 추출한 난수를 더하여 생성한 Y 데이터셋에 대해서 R로 위의 공식을 이용해서 각 시차(time lag)별로 자기상관계수를 구해보겠습니다. 



# sine wave with noise

set.seed(123)

noise <- 0.5*rnorm(50, mean=0, sd=1)

t <- seq(1, 50, 1)

y <- sin(t) + noise


plot(y, type="b")





자기상관계수를 구하고 자기상관그림을 그리는 가장 간단한 방법은 R stats 패키지에 들어있는 acf() 함수를 사용하는 것입니다. 



z <- acf(y, type=c("correlation"), plot=TRUE)


# autocorrelation coefficients per time lag k

round(z$acf, 4)


         [,1]

 [1,]  1.0000

 [2,]  0.3753

 [3,] -0.3244

 [4,] -0.5600

 [5,] -0.4216

 [6,]  0.2267

 [7,]  0.5729

 [8,]  0.3551

 [9,] -0.0952

[10,] -0.5586

[11,] -0.4582

[12,]  0.0586

[13,]  0.3776

[14,]  0.4330

[15,]  0.0961

[16,] -0.4588

[17,] -0.4759


# autocorrelation plot



좀 복잡하기는 합니다만, 위의 자기상관계수를 구하는 공식을 이용해서 R로 직접 사용자 정의 함수를 만들고, for loop 반복문을 이용해서 시차(time lag) k 별로 자기상관계수를 구할 수도 있습니다. 자기공분산을 분산으로 나누어서 표준화해주었습니다. 시차 0(time lag 0) 일 경우에는 자기자신과의 상관관계이므로 자기상관계수는 '1'이 되며, 표준화를 해주었기 때문에 자기상관계수는 -1 <= autocorr(Y, Yt-k) <= 1  사이의 값을 가집니다. 



# User Defined Function of Autocorrelation

acf_func <- function(y, lag_k){

  # y: input vector

  # lag_k : Lag order of k

  

  N = length(y) # total number of observations

  y_bar = mean(y)

  

  # Variance

  var = sum((y - y_bar)^2) / N

  

  # Autocovariance

  auto_cov = sum((y[1:(N-lag_k)] - y_bar) * (y[(1+lag_k):(N)] - y_bar)) / N

  

  # Autocorrelation coefficient = Autocovariance / Variance

  r_k = auto_cov / var

  

  return(r_k)

}



# Compute Autocorrelation per lag (from 0 to 9 in this case)

acf <- data.frame()

for (k in 0:(length(y)-1)){

  acf_k <- round(acf_func(y, lag_k = k), 4)

  acf[k+1, 'lag'] = k

  acf[k+1, 'ACF'] = acf_k

}


> print(acf)

   lag     ACF

1    0  1.0000

2    1  0.3753

3    2 -0.3244

4    3 -0.5600

5    4 -0.4216

6    5  0.2267

7    6  0.5729

8    7  0.3551

9    8 -0.0952

10   9 -0.5586

11  10 -0.4582

12  11  0.0586

13  12  0.3776

14  13  0.4330

15  14  0.0961

16  15 -0.4588

17  16 -0.4759

18  17 -0.0594

19  18  0.2982

20  19  0.3968

21  20  0.1227

22  21 -0.2544

23  22 -0.3902

24  23 -0.1981

25  24  0.1891

26  25  0.3701

27  26  0.1360

28  27 -0.1339

29  28 -0.2167

30  29 -0.1641

31  30  0.0842

32  31  0.2594

33  32  0.1837

34  33  0.0139

35  34 -0.1695

36  35 -0.1782

37  36 -0.0363

38  37  0.1157

39  38  0.1754

40  39  0.0201

41  40 -0.1428

42  41 -0.1014

43  42  0.0000

44  43  0.0750

45  44  0.0641

46  45 -0.0031

47  46 -0.0354

48  47 -0.0405

49  48 -0.0177

50  49 -0.0054

 



이번 예의 경우 시간의 흐름에 따라 싸인 파동 (sine wave) 형태를 보이는 데이터이므로, 당연하게도 시차를 두고 주기적으로 큰 양(+)의 자기상관, 큰 음(-)의 자기 상관을 보이고 있습니다. 




  (2) 자기상관계수의 95% 신뢰구간(95% confidence interval) 계산


자기상관계수의 95% 신뢰구간은 아래의 공식으로 간단하게 구할 수 있습니다. 


자기상관계수의 95% 신뢰구간 = 


이때 N은 샘플 개수, z 는 표준정규분포의 누적분포함수, 는 유의수준 5% 입니다. 


이미 익숙하게 알고 있다시피 유의수준 =0.05 에서의 z = 1.96 입니다. 이번 예에서의 관측치 개수 N = 50 입니다. 따라서 자기상관계수의 95% 신뢰구간은  이 됩니다. 



> # z quantile of 95% confidence level

> qnorm(0.975, mean=0, sd=1, lower.tail=TRUE)

[1] 1.959964

> qnorm(0.025, mean=0, sd=1, lower.tail=TRUE)

[1] -1.959964

> # sample size

> N <- length(y)

> print(N)

[1] 50

> # 95% confidence interval of autocorrelation

> qnorm(0.975, mean=0, sd=1, lower.tail=TRUE)/ sqrt(N)

[1] 0.2771808

> qnorm(0.025, mean=0, sd=1, lower.tail=TRUE)/ sqrt(N)

[1] -0.2771808




위의 acf() 함수로 그린 자기상관그림(autocorrelation plot)에서 ACF 축 +0.277 과 -0.277 위치에 수평으로 그어진 점선이 바로 95% 신뢰구간의 상, 하한 선입니다. 


위의 자기상관그림에서 보면 95% 신뢰구간의 상, 하한 점선을 주기적으로 벗어나므로 이 Y 데이터셋은 무작위적이지 않으며 (not random), 시계열적인 자기상관(autocorrelation)이 존재한다고 판단할 수 있습니다. 


만약, 무작위적인 데이터셋에 대해서 자기상관계수를 계산하고 자기상관그림을 그린다면, 시차(time lag)에 따른 자기상관계수가 모두 95% 신뢰구간 상, 하한 점선의 안쪽에 위치할 것입니다. 


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

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



728x90
반응형
Posted by Rfriend
,

이번 포스팅에서는 STEFANY COXE, STEPHEN G. WEST, AND LEONA S. AIKEN, 2009, "The Analysis of Count Data: A Gentle Introduction to Poisson Regression and Its Alternatives" 논문을 요약해서 한글로 번역한 내용을 소개하겠습니다. 

 

빈도(count data) 데이터에 대한 분석을 할 때 어떤 분포를 가정하고, 어떤 회귀모형을 적합하는 것이 좋을지 분석 방향을 잡는데 도움이 될 것입니다. 

 

 

 

 


  왜 포아송 회귀모형이 필요한가(Why Poisson Regression)?

 

특정 시간 동안 발생한 사건의 건수에 대한 도수 자료(count data, 음수가 아닌 정수)를 목표변수, 특히 평균이 10 미만일 경우, 최소제곱법(Ordinary Least Squares) 회귀모형을 적합하면 표준 오차와 유의수준이 편향되는 문제가 발생한다. OLS 회귀모형은 오차의 조건부 정규성(conditional normality), 등분산성(homoscedasiticity), 독립성(independence)을 가정한다. 하지만 평균이 작은 도수 데이터는 0 이상의 정수만 있고, 작은 계급에 많은 관측치가 몰려있으며, 우측으로 길게 치우친 분포를 띠어 회귀모형의 가정을 위배하고 음수값도 결과값으로 반환하는 문제가 있다. 종속변수가 도수 데이터이면서 오차가 정규분포, 등분산이 아닌 경우 포아송 회귀모형을 사용한다.

 

 

  • 포아송 회귀모형은 무엇인가? (Overview of Poisson Regression)

포아송 회귀모형은 일반화선형모형(Generalized Linear Model, GLM)의 한 종류로서,

 

-       Random component: Poisson Distribution, 

 

 

-       Systematic component: Mixed linear relationship

 

-       Link function: Log Link (ln)

 

 

 

표현 예

 

해석

 이며, X1을 한 단위 증가시켰을 때 효과를 보면 

  이므로, 다른 변수들을 통제했을 때 X1이 한단위 증가하면 도수의 추정치

 만큼 승법적으로 변화(multiplicative changes)한다.

 

 

선형회귀모형이 최소제곱법(OLS)로 모수를 추정한다면, 포아송 회귀모형은 최대가능도추정(MLE, Maximum Likelihood Estimation)을 통해 모수를 추정한다.

 

 

선형회귀모형은 모델에 의해 설명되는 잔차 제곱합 SS (Sum of Squares)를 총 잔차제곱합(TSS, Total Sum of Squares)로 나눈 값인  R^2로 모델 적합도(Goodness of Fit of the model) 평가하지만, 포아송 회귀모형은 설명변수의 추가에 따른 이탈도의 감소 비율(proportional reduction in deviance)로 표현되는 pseudo-R^2를 사용하여 완전모델(perfect model)에 얼마나 가깝게 적합이 되었는지를 평가한다

 

 

선형회귀모형은 F-test를 통해 모델의 유의성을 검정하지만, 포아송 회귀모형은 base modelcomplete model 간의 이탈도 차이에 대한 Chi-square test를 통해 모델 유의성을 검정한다 

선형회귀모형은 모형 타당성 평가 (Assessing model adequacy)를 위해 잔차도(residual plot)를 그려서 확인해보는 반면, 포아송 회귀모형은 이탈도 잔차도(deviance residuals vs. predicted values)를 그려보거나 실측치와 예측치를 비교해보는 것이다.

 

 

* [참고] 포아송 분포 : https://rfriend.tistory.com/101

 

 

  과대산포(Over-dispersion) 문제는 무엇이며어떻게 해결할 수 있나?

 

포아송분포는 평균과 분산이 동일(equidispersion)하다고 가정하므로, 분산이 평균보다 큰 과대산포 (Overdispersion) 시에는 포아송 회귀모형을 사용할 수 없다. 과대산포 시에 적용할 수 있는 대안 모델로는 과대산포 포아송 회귀(Overdispersed Poisson Regression), 음이항회귀(Negative Binomial Regression)가 있다.

 

 

-   Overdispersed Poisson Regression: 과대산포를 모델링하기 위해 두번째 모수인 조건부 분산 overdispersion scaling parameter

를 추정하며, 과대산포 포아송 회귀모형의 오차 분포는 

 를 가진다.

 

 

-   Negative Binomial Regression Models: 음이항회귀모형은 평균에는 영향이 없으면서 과대산포를 유발하는, 설명이 되지 않는 추가적인 가변성이 있다고 가정한다. 음이항 회귀모형은 똑같은 설명변수 값을 가지는 관측치가 다른 평균 모수를 가지고 포아송 회귀모형에 적합될 수 있도록 해준다. 오차 함수는 포아송 분포와 감마분포의 혼합 분포이다.

 

 

과대산포 검정은 가능도비 검정(Likelihood ratio test) 또는 score 검정을 사용한다.

 

모델 선택을 위해서는 AIC(Akaike Information Criterion), BIC(Bayesian Information Criterion) 을 사용하며, AIC, BIC의 값이 작은 모델을 선택하는데, 샘플 크기가 작을 때는 간단한 모델을 선택하고, 샘플 크기가 클 때는 좀더 복잡한 모델을 선택한다.

 

 

 

  과대영(Excess Zeros) 문제는 무엇이며어떻게 해결할 수 있나?

 

특정 평균을 모수로 가지는 포아송 분포에서 나타나는 ‘0’보다 많은 ‘0’을 가지는 샘플의 경우 과대영(excess zeros)이 발생했다고 말한다. 과대영은 (a) 특정 행동을 나타내지 않은 그룹의 포함 (: 비흡연자 그룹), (b) 특정 행동을 나타내지만 샘플링 계획에서 제외 (: 병원 입원한 사람만 연구) 하는 경우에 발생할 수 있다. 과대영일 경우 ZIP 회귀모형을 대안으로 사용할 수 있다.

 

-   Zero-Inflated Poisson(ZIP) Regression: 두 부분으로 나뉘는데, (1) 로지스틱 회귀모형을 사용하여 항상 ‘0’인 집단(structural zeros) 또는 항상 ‘0’은 아니지만 조사 시점에 ‘0’이라고 응답한 집단에 속할 확률을 구하고, (2) 포아송 회귀 또는 음이항 회귀모형으로 structural zeros를 포함하지 않은 나머지 관측치를 추정하는 모델을 적합한다.

 

 


  음이항분포(Negative Binomial Distribution) 

 

* Reference : https://en.wikipedia.org/wiki/Negative_binomial_distribution


 
확률 이론과 통계에서, 음 이항 분포는 특정 (무작위) 실패 횟수 (r로 표시)가 발생하기 전에 독립적이고 동일하게 분산 된 베르누이 시행 순서에서 성공 횟수의 이산 확률 분포이다. 예를 들어 1을 실패로, 1이 아닌 모든 것을 성공으로 정의하고 1이 세 번째로 나타날 때까지 반복적으로 던지면 (r = 3 번의 실패), 출현한 1이 아닌 수(non-1)의 확률 분포는 음 이항 분포가 된다.

 

 파스칼 분포 (Blaise Pascal 이후) Polya 분포 (George Pólya의 경우)는 음 이항 분포의 특수한 경우이다. 엔지니어, 기후 학자 및 다른 사람들 사이의 협약은 정수 값 정지 시간 매개 변수 r의 경우 "음 이항" 또는 "파스칼"을 사용하고, 실수 값의 경우에는 "Polya"를 사용한다. 토네이도 발생과 같은 "전염성 있는" 이산 이벤트가 발생하는 경우 Polya 분포를 사용하여 Poisson과 달리 평균 및 분산을 다르게 하여 Poisson 분포보다 더 정확한 모델을 제공 할 수 있다. "전염성 있는" 사건은 양의 공분산으로 인해 사건이 독립적이었던 경우보다 양의 상관관계가 있어 Poisson 분포에 의한 분산보다 더 큰 분산을 일으킨다. Poisson 분포는 평균과 분산이 동일하다고 가정하기 때문에 평균보다 분산이 큰 분포(Overdispersed)의 경우 적절하지 않으며, 2개의 모수를 가져서 과대산포분포를 적합할 수 있는 음 이항 분포를 사용할 수 있다.

 

 일련의 독립적 인 베르누이 재판이 있다고 가정하자. 따라서 각 임상 시험은 "성공" "실패"라는 두 가지 결과를 나타낸다. 각 시험에서 성공 확률은 p이고 실패 확률은 (1 - p)이다. 미리 정의 된 실패 횟수 r이 발생할 때까지의 서열을 관찰한다. 그러면 우리가 본 성공의 무작위 수 X는 음 이항 (또는 파스칼) 분포를 가질 것이다

 

 



 R을 활용한 Poisson model, Negative binomial model, Zero-Inflated model

 

Poisson model, Negative binomial model, Zero-Inflated model 의 유형, 분포, 모수 추정 방법, R패키지와 함수를 표로 정리하면 아래 [1]을 참고하세요. 
 

 

 

[1] Overview of count regression models

 

유형
분포
추정
방법
Description
R package & function
GLM
Poisson
ML
Poisson Regression: classical GLM, estimated by maximum likelihood(ML)
{stats} 패키지의 glm()
NB
ML
NB Regression: exteded GLM, estimated by ML including additional shape parameter
{stats} 패키지의 glm(), {MASS} 패키지의 glm.nb() 함수
Zero-augmented
Poisson
ML
Zero-Inflated Poission(ZIP)
{pscl} 패키지의 zeroinfl() 함수

 

각 회귀모형별로 R 패키지와 함수의 활용 방법은 아래와 같다.

-       (1) R을 활용한 Poisson model : glm() in the {stats} package

 



glm(formula, data, subset, na.action, weights, offset,


    family = poisson, start = NULL, control = glm.control(...),


    model = TRUE, y = TRUE, x = FALSE, ...)


 

 

-       (2) R을 활용한 Negative binomial model
   :
glm() in the {stats} package

 



glm(formula, data, subset, na.action, weights, offset,


    family = negative.binomial, start = NULL, control = glm.control(...),


    model = TRUE, y = TRUE, x = FALSE, ...)


 

 

   : glm.nb() in the {MASS} package

 



glm.nb(formula, data, weights, subset, na.action,


       start = NULL, etastart, mustart,


       control = glm.control(...), method = "glm.fit",


       model = TRUE, x = FALSE, y = TRUE, contrasts = NULL, ...,


       init.theta, link = log)


 

-       (3) R을 활용한 Zero-inflated regression model: zeroinfl() in the {pscl} package

 



zeroinfl(formula, data, subset, na.action, weights, offset,


dist = "poisson", link = "logit", control = zeroinfl.control(...),


model = TRUE, y = TRUE, x = FALSE, ...)


 

 

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

 

 

728x90
반응형
Posted by Rfriend
,

이번 포스팅에서는 마지막 네번째 포스팅으로서, 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 ratio 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번에 걸쳐서 로지스틱 회귀모형 분석에 대해 연재를 했는데요, 혹시 포스팅에 잘못된 부분이나 궁금한 부분이 있으면 댓글로 남겨주시면 수정, 보완하도록 하겠습니다. 


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



728x90
반응형
Posted by Rfriend
,

이번 포스팅은 로지스틱 회귀분석을 통한 유방암 예측(분류) 포스팅의 세번째 순서로서 '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차 변수 선별 및 목표변수와 설명변수간 관계 (탐색적) 분석'을 마치겠습니다. 


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


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



728x90
반응형
Posted by Rfriend
,

지난번 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)


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


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, VIF 10 이상인 설명변수들 중에서 가장 큰 VIF 변수 제거 -> 나머지 모든 변수에 대해 VIF 계산 -> VIF 10 이상인 설명변수들 중에서 가장 큰 VIF 변수 제거 -> 나머지 변수들 VIF 계산 -> 제거 .... 의 과정을 반복함), 


(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-collinearity 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))

  }

}





== 참고로 Python으로 다중공선성 처리하는 사용자 정의 함수는 아래 참고하세요. ==



# Remove multicollinarity recursively using Python

from statsmodels.stats.outliers_influence import variance_inflation_factor


def X_filter_multicollinearity(X, thresh=5.0):

    from datetime import datetime

    start_tm = datetime.now()

    

    variables = list(range(X.shape[1]))

    dropped = True

    while dropped:

        dropped = False

        vif = [variance_inflation_factor(X.iloc[:, variables].values, ix) 

               for ix in range(X.iloc[:, variables].shape[1])]

        

        maxloc = vif.index(max(vif))

        

        if max(vif) > thresh:

            print('==> [Dropped variable] : ' + X.iloc[:, variables].columns[maxloc])

            del variables[maxloc]

            

            if len(variables) > 1:

                dropped = True


    print('[Remaining variables] :')

    print(X.columns[variables])

    print('[Elapsed time] :', str(datetime.now() - start_tm))

    

    return variables


# run the UDF

X_remained_idx = X_filter_multicollinearity(X)

print('X index after filtering multicollinearity:', X_remained_idx)

 




위의 사용자정의함수 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차 변수 선택 및 목표변수와 설명변수 간 관계 분석'에 대해서 다루어보겠습니다. 


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



728x90
반응형
Posted by Rfriend
,

이번 포스팅부터 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 데이터에 대한 탐색적 데이터 분석과 데이터 전처리에 대해서 다루겠습니다. 


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



728x90
반응형
Posted by Rfriend
,

이번 포스팅에서는 카이제곱 적합도 검정(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)가 매우 유사한 분포를 띠고 있음을 눈으로도 확인할 수 있습니다.

 

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

 

728x90
반응형
Posted by Rfriend
,

관측값이 질적 자료(qualitative data) 또는 어떤 속성에 따라 분류되어 범주(category)에 속하는 도수(frequency)로 주어질 경우 이를 범주형 자료(categorical data) 라고 합니다. 범주형 자료의 예로는 학력(초등졸, 중등졸, 고등졸, 대졸, 대학원졸), 연수익(극빈, 하, 중, 상, 극상) 등이 있습니다.

 

지난 포스팅에서는 범주형 자료분석 중에서 (1) 적합도 검정, (2) 독립성 검정 (test of independence)에 대해서 소개하였다면, 이번 포스팅에서는 (3) 동질성 검정(test of homogeneity)을 다루도록 하겠습니다.

 

범주형 자료 분석 유형별 간략한 소개는 아래와 같습니다.

 

 

 

 

[ 범주형 자료 분석 (categorical data test) ]

 

(1) 적합도 검정(goodness of fit test) : 관측값들이 어떤 이론적 분포를 따르고 있는지를 검정. 한 개의 요인을 대상으로 함 

 

(2) 독립성 검정(test of independence) : 서로 다른 요인들에 의해 분할되어 있는 경우 그 요인들이 관찰값에 영향을 주고 있는지 아닌지, 요인들이 서로 연관이 있는지 없는지를 검정. 두 개의 요인을 대상으로 함.

 

(3) 동질성 검정(test of homogeneity) : 관측값들이 정해진 범주 내에서 서로 비슷하게 나타나고 있는지를 검정. 속성 A, B를 가진 부모집단(subpopulation) 각각으로부터 정해진 표본의 크기만큼 자료를 추출하는 경우에 분할표에서 부모집단의 비율이 동일한가를 검정. 두 개의 요인을 대상으로 함.  

 

 

 

 

 

 

 

독립성 검정과 동질성 검정의 이해를 돕기 위해 서로 비교를 해보자면,

 

(1) 독립성 검정이 두 변수 X와 Y가 서로 독립인지 아닌지에 대한 판단이라면, 동질성 검정은 r개의 행과 c개의 열을 가진 두 변수 X와 Y로부터 작성된 분할표의 각 열분포에서 행들이 균일한 값을 가지는지 즉, 각 열에서 행들의 동질성(homegeneity)를 검정하는 것입니다. 두 검정법의 이런 차이점은 개념상의 차이일 뿐이며 검정을 하는 방법은 카이제곱 검정을 이용해서 동일합니다.

 

(2) 독립성 검정은 하나의 모집단에서 표본을 무작위로 추출한 후 추출된 표본을 두가지 속성(변수)에 따라 분류합니다.  반면에 동질성 검정은 부모집단(subpopulation)을 먼저 설정한 후 각 부모집단으로부터 정해진 표본의 크기만큼 무작위로 추출하여 분할표에서 부모집단의 비율이 동일한가를 검정하게 됩니다. 가령, 소득수준에 따라 지지 정당이 동일한지 여부를 검정한다고 할 때, 우선 소득수준을 부모집단으로 설정하고, 각 소득수준별로 정해진 크기의 표본을 무작위로 추출하는 식입니다.

 

 

r개의 행과 c개의 열을 가진 두 변수 X와 Y로부터 작성된 r*c 분할표를 이용한 동질성 검정을 위한 데이터셋 구조는 아래와 같습니다.

 

 

[ 동질성 검정 자료 구조 (dataset for test of homogeneity) ]

 

 

 

동질성 검정을 위한 가설과 검정통계량, 검정방법은 아래와 같습니다. 검정통계량 X^2 는 귀무가설 H0가 사실일 때 근사적으로 자유도 (r-1)(c-1)인 카이제곱 분포를 따르는 것으로 알려져있습니다.

 

 

[ 동질성 검정 가설 및 검정 통계량, 검정 방법 ]

 

(1) 가설

 

  - 귀무가설 H0 : p1j = p2j = ... = prj,   j = 1, ..., c

 

  - 대립가설 H1 : H0가 아니다

 

 

(2) 검정 통계량

 

 

(3) 검정 방법

 

 

 

이제 아래의 문제를 R의 chisq.test() 함수를 이용해서 풀어보도록 하겠습니다.

 

 

(문제) 초등학교 1학년 남학생 100명과 여학생 200명을 무작위로 추출하여 TV 프로그램 선호도를 조사하였다. 유의수준 α 0.05 에서 남학생의 TV 프로그램 선호도와 여학생의 TV프로그램 선호도가 동일한지 검정하여라.

 

 

선호 TV 프로그램 

row total 

 뽀로로

짱구는 못말려 

로봇카 폴리 

 남학생

 50

30

20 

100

 여학생

 50

80

70 

200

 column total

100

110

90

300

 

 

 

 

[가설]
- H0 : 남학생과 여학생별로 TV 선호도는 동일하다
         (p1j = p2j,   j = 뽀로로, 짱구는 못말려, 로봇카 폴리)

- H1 : H0가 아니다

 

 

아래 분석에 사용한 R chisq.test() 함수는 이전 포스팅의 독립성 검정과 동일합니다.

 

 

> ##---------------------------------------------------------------------
> ## categorical data analysis - (3) test of homogeneity : chisq.test()
> ##---------------------------------------------------------------------
> 
> ##-------------
> # (a) textbook problem
> 
> ## data key-in
> # data key-in way 1 : rbind()
> row_1 <- c(50, 30, 20)
> row_2 <- c(50, 80, 70)
> 
> data_rbind <- rbind(row_1, row_2)
> data_rbind
      [,1] [,2] [,3]
row_1   50   30   20
row_2   50   80   70
> 
> 
> # data key-in way 2 : matrix()
> raw_data <- c(50, 30, 20, 50, 80, 70)
> data_matrix <- matrix(raw_data, byrow=TRUE, nrow=2)
> data_matrix
     [,1] [,2] [,3]
[1,]   50   30   20
[2,]   50   80   70
> 
> 
> # giving names to the rows and columns of the data table : dimnames()
> dimnames(data_matrix) <- list("Gender" = c("Boys", "Girls"), 
+                               "TV_Preferences" = c("Pororo", "JJangGu", "RobotCar"))
> 
> data_matrix
       TV_Preferences
Gender  Pororo JJangGu RobotCar
  Boys      50      30       20
  Girls     50      80       70
> 
> 
> ## exploratory data analysis
> # marginal distribution : addmargins()
> addmargins(data_matrix)
       TV_Preferences
Gender  Pororo JJangGu RobotCar Sum
  Boys      50      30       20 100
  Girls     50      80       70 200
  Sum      100     110       90 300
> 
> 
> # proportional distribution : prop.table()
> prop.table(data_matrix)
       TV_Preferences
Gender     Pororo   JJangGu   RobotCar
  Boys  0.1666667 0.1000000 0.06666667
  Girls 0.1666667 0.2666667 0.23333333
> 
> addmargins(prop.table(data_matrix))
       TV_Preferences
Gender     Pororo   JJangGu   RobotCar       Sum
  Boys  0.1666667 0.1000000 0.06666667 0.3333333
  Girls 0.1666667 0.2666667 0.23333333 0.6666667
  Sum   0.3333333 0.3666667 0.30000000 1.0000000
> 
> 
> # bar plot : barplot()
> barplot(t(data_matrix), beside=TRUE, legend=TRUE, 
+         ylim=c(0, 120), 
+         ylab="Observed frequencies in sample", 
+         main="TV viewing preferences by gender")
> 

 

> 
> ## chisquared test : chisq.test()
> chisq.test(data_matrix)

	Pearson's Chi-squared test

data:  data_matrix
X-squared = 19.318, df = 2, p-value = 6.384e-05

> 
> 
> # indexing statistics of chisq.test()
> chisq.test_output_2 <- chisq.test(data_matrix)
> 
> chisq.test_output_2$observed # observed frequency
       TV_Preferences
Gender  Pororo JJangGu RobotCar
  Boys      50      30       20
  Girls     50      80       70
> chisq.test_output_2$expected # expected frequeycy
       TV_Preferences
Gender    Pororo  JJangGu RobotCar
  Boys  33.33333 36.66667       30
  Girls 66.66667 73.33333       60
> chisq.test_output_2$residuals # residual between observed and expected frequecy
       TV_Preferences
Gender     Pororo    JJangGu  RobotCar
  Boys   2.886751 -1.1009638 -1.825742
  Girls -2.041241  0.7784989  1.290994
> 
> chisq.test_output_2$statistic # chi-squared statistics
X-squared 
 19.31818 
> chisq.test_output_2$parameter # degrees of freedom
df 
 2 
> chisq.test_output_2$p.value # P-value
[1] 6.384253e-05

 

 

 

위의 분석 결과를 해석해보자면, 카이제곱 통계량 값이 19.318이 나왔고 P-value가 6.384e-05 로서 유의수준 α 0.05 보다 훨씬 작기때문에 귀무가설 H0 를 기각하고 대립가설 H1을 채택하여 "남학생/여학생별 선호하는 TV프로그램은 동일하지 않다"고 판단할 수 있겠습니다.

 

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

 

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

 

728x90
반응형
Posted by Rfriend
,

지난 포스팅에서는 범주형 자료분석 중에서 (1) 적합도 검정에 대해서 소개하였다면, 이번 포스팅에서는 (2) 독립성 검정 (test of independence)에 대해서 알아보겠으며, 다음번 포스팅에서는 (3) 동질성 검정을 다루도록 하겠습니다.

 

 

범주형 자료 분석 유형별 간략한 소개는 아래와 같습니다.

 

 

 

[ 범주형 자료 분석 (categorical data test) ]

 

(1) 적합도 검정(goodness of fit test) : 관측값들이 어떤 이론적 분포를 따르고 있는지를 검정. 한 개의 요인을 대상으로 함 

 

(2) 독립성 검정(test of independence) : 서로 다른 요인들에 의해 분할되어 있는 경우 그 요인들이 관찰값에 영향을 주고 있는지 아닌지, 요인들이 서로 연관이 있는지 없는지를 검정. 두 개의 요인을 대상으로 함.

 

(3) 동질성 검정(test of homogeneity) : 관측값들이 정해진 범주 내에서 서로 비슷하게 나타나고 있는지를 검정. 속성 A, B를 가진 부모집단(subpopulation) 각각으로부터 정해진 표본의 크기만큼 자료를 추출하는 경우에 분할표에서 부모집단의 비율이 동일한가를 검정. 두 개의 요인을 대상으로 함.

 

 

 

 

 

 

 

 

 

독립성 검정(test of independence)은 두 개의 범주형 변수/요인(2 factors)이 서로 연관성이 있는지, 상관이 있는지, 독립적인지를 카이제곱 검정(chisquared test)을 통해 통계적으로 판단하는 방법입니다.

 

가령, 학력(초등졸, 중등졸, 고등졸, 대졸, 대학원졸)이라는 범주형 변수(variable X)/요인(factor 1)와 연소득(하, 중, 상)이라는 범주형 변수(variable Y)/요인(factor 2) 간에 서로 관련성이 있는 것인지 아니면 관련이 없이 독립적인지를 판단하는 것과 같은 문제에 독립성 검정을 사용합니다.

 

참고로, 두 변수가 양적변수(qualitative variable)인 경우 두 변수 간 상관관계 분석을 위해서는 공분산 분석, 상관계수 분석, 회귀분석 등을 활용합니다.

 

범주형 자료분석의 경우 두 변수의 관련성을 보려면 분할표를 만들어서 카이제곱 검정을 하게 되는데요, 두 변수간의 관련성을 양적변수 분석할 때처럼 숫자로 얼마나 관련성이 큰지를 알 수 있는 통계량을 제공하지는 않습니다.  (범주형 자료분석에서 두 변수 간 관련성 측도로 파이계수, 속성계수, 크레머 V 등의 통계량이 있는데요, 이번 포스팅에서는 이에 대한 설명은 건너뛰겠습니다.)

 

 

자료를 분류하는 두 변수를 x와 Y라고 하고, 변수 X는 m개, 변수 Y는 n개의 범주(혹은 계급 class, 혹은 요인 수준 factor level)를 가진다고 했을 때 관측도수 Oij 는 m개와 n개의 층으로 이루어진 아래와 같은 표로 정리할 수 있습니다.  이를 m*n 분할표 (m*n contingency table)이라고 부릅니다.

 

 

[ 독립성 검정 자료 구조 (dataset for test of independence) ]

 

 

 

 

독립성 검정에는 카이제곱 X^2 통계량을 사용하는데요, 귀무가설 H0 가 사실일 때 자유도 (m-1)(n-1)인 카이제곱분포에 근사하는 것으로 알려져 있습니다. (☞ 카이제곱 분포(Chi-squared distribution) 참고 포스팅

 

검정통계량 카이제곱 X^2 은 각 범주의 기대도수가 5 이상인 경우에 사용하는 것이 바람직하며, 기대도수가 5 미만인 경우에는 주의를 요합니다. (5보다 작으면 인접 범주와 합치는 것도 방법)

 

기본 원리는, 관측도수 O11, O21, ..., Omn 이 기대도수 E11, E21, ..., Emn 과 차이가 없다면 검정통계량 X0^2 값이 '0'이 되고, 반대로 관측도수와 기대도수가 차이가 크다면 검정통계량 값 또한 커지게 된다는 것입니다.

 

 

 

 

 

[ 독립성 검정(test of independence) 가설 및 검정 통계량, 검정 방법 ]

 

(1) 가설 (hypothesis)

 

  - 귀무가설 H0 : 두 변수 X와 Y는 서로 독립이다 (관련성이 없다)
                        ( pij = pim * pnj,   i = 1, 2, .., m,   j = 1, 2, ..., n )

 

  - 대립가설 H1 : 두 변수 X와 Y는 서로 독립이 아니다 (관련성이 있다)

 

 

 

(2) 검정 통계량 (chisquared test statistics)

 

 

 

(3) 검정 방법 (test method)

 

 

  • (a) chisq.test() of data from text problem

 

아래 문제에 대해서 R의 chisq.test() 함수를 사용해서 풀어보도록 하겠습니다. 

 

 

(문제)  학급 (class 1, class 2, class 3)과 수학 성적 (math score High, Middle, Low, Fail) 간의 관련성이 있는지를 조사한 아래의 분할표를 사용하여 유의수준 α 0.05 로 검정하여라.  

 

 

학급과 수학성적 분할표 (contingency table of class & math score)

 

 

score High 

score Middle 

score Low 

Fail 

Class 1

7

13

9

12

 Class 2

13

21

10

19

 Class 3

11

18

12

13

 

 

 

먼저 데이터 입력 및 탐색적 분석을 위한 R 함수입니다.

 

 

> ##---------------------------------------------------------------------
> ## categorical data analysis - (2) test of independence : chisq.test()
> ##---------------------------------------------------------------------
> 
> ##-------------
> # (a) textbook problem
> 
> ## data key-in
> # data key-in way 1 : rbind()
> row_1 <- c(7, 13, 9, 12)
> row_2 <- c(13, 21, 10, 19)
> row_3 <- c(11, 18, 12, 13)
> 
> data_rbind <- rbind(row_1, row_2, row_3)
> data_rbind
      [,1] [,2] [,3] [,4]
row_1    7   13    9   12
row_2   13   21   10   19
row_3   11   18   12   13
> 
> 
> # data key-in way 2 : matrix()
> raw_data <- c(7, 13, 9, 12, 13, 21, 10, 19, 11, 18, 12, 13)
> data_matrix <- matrix(raw_data, byrow=TRUE, nrow=3)
> data_matrix
     [,1] [,2] [,3] [,4]
[1,]    7   13    9   12
[2,]   13   21   10   19
[3,]   11   18   12   13
> 
> 
> # giving names to the rows and columns of the data table : dimnames()
> dimnames(data_matrix) <- list("Class" = c("Class_1", "Class_2", "Class_3"), 
+                               "Score" = c("Score_H", "Score_M", "Score_L", "Fail"))
> 
> data_matrix
         Score
Class     Score_H Score_M Score_L Fail
  Class_1       7      13       9   12
  Class_2      13      21      10   19
  Class_3      11      18      12   13
> 
> 
> ## exploratory data analysis
> # marginal distribution : addmargins()
> addmargins(data_matrix)
         Score
Class     Score_H Score_M Score_L Fail Sum
  Class_1       7      13       9   12  41
  Class_2      13      21      10   19  63
  Class_3      11      18      12   13  54
  Sum          31      52      31   44 158
> 
> 
> # proportional distribution : prop.table()
> prop.table(data_matrix)
         Score
Class        Score_H    Score_M    Score_L       Fail
  Class_1 0.04430380 0.08227848 0.05696203 0.07594937
  Class_2 0.08227848 0.13291139 0.06329114 0.12025316
  Class_3 0.06962025 0.11392405 0.07594937 0.08227848
> 
> addmargins(prop.table(data_matrix))
         Score
Class        Score_H    Score_M    Score_L       Fail       Sum
  Class_1 0.04430380 0.08227848 0.05696203 0.07594937 0.2594937
  Class_2 0.08227848 0.13291139 0.06329114 0.12025316 0.3987342
  Class_3 0.06962025 0.11392405 0.07594937 0.08227848 0.3417722
  Sum     0.19620253 0.32911392 0.19620253 0.27848101 1.0000000
> 
> 
> # bar plot : barplot()
> barplot(t(data_matrix), beside=TRUE, legend=TRUE, 
+         ylim=c(0, 30), 
+         ylab="Observed frequencies in sample", 
+         main="Frequency of math score by class")
> 

 

 

 

 

 

다음으로 카이제곱 검정 및 통계량 indexing 방법입니다.

 

 

> 
> ## chisquared test : chisq.test()
> chisq.test(data_matrix)

	Pearson's Chi-squared test

data:  data_matrix
X-squared = 1.3859, df = 6, p-value = 0.9667

> 
> 
> # indexing statistics of chisq.test()
> chisq.test_output_2 <- chisq.test(data_matrix)
> 
> chisq.test_output_2$observed # observed frequency
         Score
Class     Score_H Score_M Score_L Fail
  Class_1       7      13       9   12
  Class_2      13      21      10   19
  Class_3      11      18      12   13
> chisq.test_output_2$expected # expected frequeycy
         Score
Class       Score_H  Score_M   Score_L     Fail
  Class_1  8.044304 13.49367  8.044304 11.41772
  Class_2 12.360759 20.73418 12.360759 17.54430
  Class_3 10.594937 17.77215 10.594937 15.03797
> chisq.test_output_2$residuals # residual between observed and expected frequecy
         Score
Class        Score_H     Score_M    Score_L       Fail
  Class_1 -0.3681990 -0.13439170  0.3369579  0.1723221
  Class_2  0.1818200  0.05837794 -0.6714739  0.3475383
  Class_3  0.1244439  0.05404747  0.4316649 -0.5255380
> 
> chisq.test_output_2$statistic # chi-squared statistics
X-squared 
 1.385926 
> chisq.test_output_2$parameter # degrees of freedom
df 
 6 
> chisq.test_output_2$p.value # P-value
[1] 0.966709

  

 

 

위 분석결과를 보면 P-value가 0.966709 로서 유의수준 α 0.05보다 크므로 귀무가설 H0를 채택하여 "학급과 수학성적 간에는 서로 관련성이 없다. 즉, 독립적이다"고 판단할 수 있겠습니다.

 

 

 

이상으로 독립성 검정(test of independence)을 카이제곱 검정 기법을 사용해서 하는 방법을 소개하였습니다.

 

 

 

 

다음번 포스팅에서는 (3) 동질성 검정(test of homegeneity)에 대해서 알아보도록 하겠습니다.

 

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

 

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

 

 

 

 

728x90
반응형
Posted by Rfriend
,

관측값이 질적 자료(qualitative data) 또는 어떤 속성에 따라 분류되어 범주(category)에 속하는 도수(frequency)로 주어질 경우 이를 범주형 자료(categorical data) 라고 합니다. 범주형 자료의 예로는 학력(초등졸, 중등졸, 고등졸, 대졸, 대학원졸), 연수익(극빈, 하, 중, 상, 극상) 등이 있습니다.

 

앞서의 포스팅에서는 종속변수가 연속형 자료(continuous data)인 경우에 사용하는 검정 방법으로 t-Test와 ANOVA에 대해서 소개하였습니다.

 

이번 포스팅부터는 종속변수가 범주형 자료(categorical data)인 경우에 사용하는 분석기법으로 카이제곱 검정(Chi-Squared Test)에 대해서 알아보도록 하겠습니다.

 

범주형 자료 분석은 크게 적합도 검정(goodness f fit test), 독립성 검정(test of independence), 동질성 검정(test of homogeneity)의 3가지로 분류할 수 있으며, 이번 포스팅에서는 (1) 적합도 검정에 대해서 알아보도록 하겠습니다.

 

 

 

[ 범주형 자료 분석 (categorical data test) ]

 

(1) 적합도 검정(goodness of fit test) : 관측값들이 어떤 이론적 분포를 따르고 있는지를 검정. 한 개의 요인을 대상으로 함 

 

(2) 독립성 검정(test of independence) : 서로 다른 요인들에 의해 분할되어 있는 경우 그 요인들이 관찰값에 영향을 주고 있는지 아닌지, 요인들이 서로 연관이 있는지 없는지를 검정. 두 개의 요인을 대상으로 함.

 

(3) 동질성 검정(test of homogeneity) : 관측값들이 정해진 범주 내에서 서로 비슷하게 나타나고 있는지를 검정. 속성 A, B를 가진 부모집단(subpopulation) 각각으로부터 정해진 표본의 크기만큼 자료를 추출하는 경우에 분할표에서 부모집단의 비율이 동일한가를 검정. 두 개의 요인을 대상으로 함.

 

 

 

 

 

 

적합도 검정(goodness of fit test)은 k개의 범주 (혹은 계급)을 가지는 한 개의 요인(factor)에 대해서 어떤 이론적 분포를 따르고 있는지를 검정하는 방법입니다. 

 

기본 원리는, 도수분포의 각 구간에 있는 관측도수를 O1, O2, ..., Ok 라 하고, 각 범주 ( 혹은 계급)가 일어날 확률을 p1, p2, ..., pk 라고 할 때 기대되는 관측도수 E1, E2, ..., Ek 를 계산하여 실제 관측도수와 기대 관측도수의 차이를 카이제곱 검정 통계량(Chi-squared statistics)을 활용하여 가정한 확률모형에 적합한지를 평가하게 됩니다. 만약 귀무가설 H0가 맞다면 관측도수와 기대도수가 별 차이가 없을 것이므로 검정통계량 X0^2 값이 작을 것이며, 반대로 대립가설 H1이 맞다면 관측도수와 기대도수의 차이가 클 것이므로 검정통계량 X0^2 값이 커질 것입니다.

 

 

 

[ 적합도 검정 가설 및 검정 통계량, 검정 방법 ]

(1) 가설

 

  - 귀무가설 H0 : 관측값의 도수와 가정한 이론 도수(기대 관측도수)가 동일하다
                        ( p1 = p10, p2 = p20, ..., pk = pko )

 

  - 대립가설 H1 : 적어도 하나의 범주 (혹은 계급)의 도수가 가정한 이론 도수(기대 관측도수)와 다르다

                        (적어도 하나의 pi는 가정된 값 pi0과 다르다)

 

 

(2) 검정 통계량

 

 

 

(3) 검정 방법

 

 

 

 

  • (a) chisq.test() of data from text problem

 

아래 문제에 대해서 R의 chisq.test() 함수를 사용해서 풀어보도록 하겠습니다.

 

 

 

(문제)  유전학자 멘델은 콩 교배에 대한 유전의 이론적 모형으로서 잡종비율을 A : B : C = 2 : 3 : 5 라고 주장하였다.  이 이론의 진위를 가리기 위해 두 콩 종자의 교배로 나타난 100개의 콩을 조사하였더니 A형 19개, B형 41개, C형 40개였다.  이러한 관찰값을 얻었을 때 멘델 유전학자의 이론이 맞다고 할 수 있는지를 유의수준 α = 0.05 에서 검정하여라.

 

 

 

> ##---------------------------------------------------------------------
> ## categorical data analysis - (1) goodness of fit test : chisq.test()
> ##---------------------------------------------------------------------
> 
> obs <- c(19, 41, 40)
> null.probs <- c(2/10, 3/10, 5/10)
> 
> chisq.test(obs, p=null.probs)

	Chi-squared test for given probabilities

data:  obs
X-squared = 6.0833, df = 2, p-value = 0.04776

 

 

위 분석결과를 보면 P-value가 0.04776 이므로 유의수준 α 0.05보다 작으므로 귀무가설 H0를 기각하고 대립가설 H1을 채택하여 "멘델이 주장한 콩의 잡종비율 이론적 분포는 적합하지 않다"고 판단할 수 있겠습니다.

 

 

 참고로, R로 통계분석을 하면 콘솔 창에 보여지는 내용 말고도 실제로 다양한 통계량이 계산이 되어 list 형태로 메모리상에 가지고 있으며 단지 눈에 보이지 않을 뿐인데요, indexing 기법을 활용하면 chisq.test() 함수로 카이제곱 검정 실행한 후에 다양한 통계량들을 선별해서 볼 수도 있고, 통계량을 다른 분석 혹은 애플리케이션에 input으로 넣어 재활용할 수도 있습니다.  R의 큰 장점 중의 하나이니 팁으로 알아두면 좋겠습니다. 주요 통계량 몇 개를 아래에 소개합니다.

 

 

> # To see results of chisquared test
> chisq.test_output_1 <- chisq.test(obs, p=null.probs)
> 
> chisq.test_output_1$observed # observed frequency
[1] 19 41 40
> chisq.test_output_1$expected # expected frequeycy
[1] 20 30 50
> chisq.test_output_1$residuals # residual between observed and expected frequecy
[1] -0.2236068  2.0083160 -1.4142136
> 
> chisq.test_output_1$statistic # chi-squared statistics
X-squared 
 6.083333 
> chisq.test_output_1$parameter # degrees of freedom
df 
 2 
> chisq.test_output_1$p.value # P-value
[1] 0.04775523

 

 

 

참고로 하나더, 검정통계량 X^2는 귀무가설 H0가 참이라는 가정 하에 근사적으로 자유도가 k-1 인 카이제곱분포를 따르는 것으로 알려져 있습니다.  (☞ 카이제곱 분포(Chi-squared distribution) 참고 포스팅)

 

카이제곱 검정에 사용하는 카이제곱 분포는 범주형 자료의 도수 추정에 사용되는데요, 이때 도수가 너무 작으면 "카이제곱 approximation)이 정확하지 않을 수도 있습니다" 라는 경고메시지가 뜹니다.

 

 

> # Warning message when there are not sufficient frequencies > # R will issue a warning message if any of the EFs fall below 5

> obs_2 <- c(5, 5) > null.probs_2 <- c(0.3, 0.7) > > chisq.test(obs_2, p=null.probs_2) Chi-squared test for given probabilities data: obs_2 X-squared = 1.9048, df = 1, p-value = 0.1675 Warning message: In chisq.test(obs_2, p = null.probs_2) : 카이제곱 approximation은 정확하지 않을수도 있습니다

 

 

 

 

  • (b) chisqtest() of data from a table object

위의 문제는 텍스트 문제로 도수와 확률이 주어지면 'x'와 'p'를 직접 입력하였는데요, 데이터가 Table object 형태로 주어졌을 때 카이제곱 검정으로 적합도 검정하는 방법을 소개하겠습니다.

 

R에 내장된 HairEyeColor table 데이터셋에 있는 Hair 요인 변수를 대상으로, Hair의 요인 수준(factor levels, 혹은 계급 class) 별로 생물학자가 주장하기를 확률이 Black 20%, Brown 50%, Red 10%, Blond 20% 라고 하는데요, 유의수준 α 0.05 로 검정을 해보겠습니다.

 

> ##------------
> # chisq.test() of data from a table object
> str(HairEyeColor) 
 table [1:4, 1:4, 1:2] 32 53 10 3 11 50 10 30 10 25 ...
 - attr(*, "dimnames")=List of 3
  ..$ Hair: chr [1:4] "Black" "Brown" "Red" "Blond"
  ..$ Eye : chr [1:4] "Brown" "Blue" "Hazel" "Green"
  ..$ Sex : chr [1:2] "Male" "Female"
> 
> HairEyeColor 
, , Sex = Male

       Eye
Hair    Brown Blue Hazel Green
  Black    32   11    10     3
  Brown    53   50    25    15
  Red      10   10     7     7
  Blond     3   30     5     8

, , Sex = Female

       Eye
Hair    Brown Blue Hazel Green
  Black    36    9     5     2
  Brown    66   34    29    14
  Red      16    7     7     7
  Blond     4   64     5     8

> 
> dimnames(HairEyeColor)
$Hair
[1] "Black" "Brown" "Red"   "Blond"

$Eye
[1] "Brown" "Blue"  "Hazel" "Green"

$Sex
[1] "Male"   "Female"

> 
> margin.table(HairEyeColor, 1) # Hair
Hair
Black Brown   Red Blond 
  108   286    71   127 
> margin.table(HairEyeColor, 2) # Eye
Eye
Brown  Blue Hazel Green 
  220   215    93    64 
> margin.table(HairEyeColor, 3) # Sex
Sex
  Male Female 
   279    313 
> 
> # vector of observed frequencies and probabilities
> Hair_Freq <- c(margin.table(HairEyeColor, 1))
> Hair_Freq
Black Brown   Red Blond 
  108   286    71   127 
> Hair_Prob <- c(0.2, 0.5, 0.1, 0.2)
> 
> chisq.test(x=Hair_Freq, p=Hair_Prob)

	Chi-squared test for given probabilities

data:  Hair_Freq
X-squared = 4.228, df = 3, p-value = 0.2379

 

 

 

  • (c) chisq.test() of data from Data Frame

 

데이터가 텍스트, Table object가 아니고 Data Frame일 경우에 카이제급 검정하는 방법도 소개해드리겠습니다.  MASS 패키지에 내장된 Cars93 Data Frame 을 이용하며, 자동차종(Type)의 이론상 분포 Compact 20%,  Large 10%, Midsize 20%, Small 20%, Sporty  20%, Van 10% 에 대해서 유의수준 α 0.05로 검정해보겠습니다. 

data(data frame, package="xxx"), table() 함수를 이용합니다.

 

> ##------------
> # chisq.test() of data from data frame
> data(Cars93, package="MASS")
> head(Cars93)
  Manufacturer   Model    Type Min.Price Price Max.Price MPG.city MPG.highway            AirBags
1        Acura Integra   Small      12.9  15.9      18.8       25          31               None
2        Acura  Legend Midsize      29.2  33.9      38.7       18          25 Driver & Passenger
3         Audi      90 Compact      25.9  29.1      32.3       20          26        Driver only
4         Audi     100 Midsize      30.8  37.7      44.6       19          26 Driver & Passenger
5          BMW    535i Midsize      23.7  30.0      36.2       22          30        Driver only
6        Buick Century Midsize      14.2  15.7      17.3       22          31        Driver only
  DriveTrain Cylinders EngineSize Horsepower  RPM Rev.per.mile Man.trans.avail Fuel.tank.capacity
1      Front         4        1.8        140 6300         2890             Yes               13.2
2      Front         6        3.2        200 5500         2335             Yes               18.0
3      Front         6        2.8        172 5500         2280             Yes               16.9
4      Front         6        2.8        172 5500         2535             Yes               21.1
5       Rear         4        3.5        208 5700         2545             Yes               21.1
6      Front         4        2.2        110 5200         2565              No               16.4
  Passengers Length Wheelbase Width Turn.circle Rear.seat.room Luggage.room Weight  Origin
1          5    177       102    68          37           26.5           11   2705 non-USA
2          5    195       115    71          38           30.0           15   3560 non-USA
3          5    180       102    67          37           28.0           14   3375 non-USA
4          6    193       106    70          37           31.0           17   3405 non-USA
5          4    186       109    69          39           27.0           13   3640 non-USA
6          6    189       105    69          41           28.0           16   2880     USA
           Make
1 Acura Integra
2  Acura Legend
3       Audi 90
4      Audi 100
5      BMW 535i
6 Buick Century
> str(Cars93)
'data.frame':	93 obs. of  27 variables:
 $ Manufacturer      : Factor w/ 32 levels "Acura","Audi",..: 1 1 2 2 3 4 4 4 4 5 ...
 $ Model             : Factor w/ 93 levels "100","190E","240",..: 49 56 9 1 6 24 54 74 73 35 ...
 $ Type              : Factor w/ 6 levels "Compact","Large",..: 4 3 1 3 3 3 2 2 3 2 ...
 $ Min.Price         : num  12.9 29.2 25.9 30.8 23.7 14.2 19.9 22.6 26.3 33 ...
 $ Price             : num  15.9 33.9 29.1 37.7 30 15.7 20.8 23.7 26.3 34.7 ...
 $ Max.Price         : num  18.8 38.7 32.3 44.6 36.2 17.3 21.7 24.9 26.3 36.3 ...
 $ MPG.city          : int  25 18 20 19 22 22 19 16 19 16 ...
 $ MPG.highway       : int  31 25 26 26 30 31 28 25 27 25 ...
 $ AirBags           : Factor w/ 3 levels "Driver & Passenger",..: 3 1 2 1 2 2 2 2 2 2 ...
 $ DriveTrain        : Factor w/ 3 levels "4WD","Front",..: 2 2 2 2 3 2 2 3 2 2 ...
 $ Cylinders         : Factor w/ 6 levels "3","4","5","6",..: 2 4 4 4 2 2 4 4 4 5 ...
 $ EngineSize        : num  1.8 3.2 2.8 2.8 3.5 2.2 3.8 5.7 3.8 4.9 ...
 $ Horsepower        : int  140 200 172 172 208 110 170 180 170 200 ...
 $ RPM               : int  6300 5500 5500 5500 5700 5200 4800 4000 4800 4100 ...
 $ Rev.per.mile      : int  2890 2335 2280 2535 2545 2565 1570 1320 1690 1510 ...
 $ Man.trans.avail   : Factor w/ 2 levels "No","Yes": 2 2 2 2 2 1 1 1 1 1 ...
 $ Fuel.tank.capacity: num  13.2 18 16.9 21.1 21.1 16.4 18 23 18.8 18 ...
 $ Passengers        : int  5 5 5 6 4 6 6 6 5 6 ...
 $ Length            : int  177 195 180 193 186 189 200 216 198 206 ...
 $ Wheelbase         : int  102 115 102 106 109 105 111 116 108 114 ...
 $ Width             : int  68 71 67 70 69 69 74 78 73 73 ...
 $ Turn.circle       : int  37 38 37 37 39 41 42 45 41 43 ...
 $ Rear.seat.room    : num  26.5 30 28 31 27 28 30.5 30.5 26.5 35 ...
 $ Luggage.room      : int  11 15 14 17 13 16 17 21 14 18 ...
 $ Weight            : int  2705 3560 3375 3405 3640 2880 3470 4105 3495 3620 ...
 $ Origin            : Factor w/ 2 levels "USA","non-USA": 2 2 2 2 2 1 1 1 1 1 ...
 $ Make              : Factor w/ 93 levels "Acura Integra",..: 1 2 4 3 5 6 7 9 8 10 ...
> 
> Car_Type <- table(Cars93$Type)
> Car_Type

Compact   Large Midsize   Small  Sporty     Van 
     16      11      22      21      14       9 
> Car_Type_Prob <- c(0.2, 0.1, 0.2, 0.2, 0.2, 0.1)
> 
> chisq.test(x=Car_Type, p=Car_Type_Prob)

	Chi-squared test for given probabilities

data:  Car_Type
X-squared = 2.7527, df = 5, p-value = 0.738

 

 

 

이상으로 다양한 형태의 데이터셋을 활용해서 적합도 검정(goodness of fit test)을 카이제곱 검정 기법을 사용해서 하는 방법을 소개하였습니다.

 

=> 카이제곱 적합도 검정으로 관측치가 포아송분포로 부터의 데이터인지 여부를 검정하는 예시는 http://rfriend.tistory.com/362 를 참고하세요.

 

다음번 포스팅에서는 (2) 독립성 검정(test of independence) , (3) 동질성 검정 (test of homogeneity)대해서 알아보도록 하겠습니다.

 

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

 

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

 

 

728x90
반응형
Posted by Rfriend
,