분류 전체보기

  1. 구글 등에서 Mendeley를 검색하여 공식 홈페이지에서 데스크탑 버전을 먼저 다운 받는다.

  2. 설치를 완료하여 데스크탑 버전을 킨 다음, [ Tool 메뉴 > install MS Word Plugin ]을 다운받는다.

    • 기존에 Mendeley에 원하는 참고논문들을 추가해 넣은 상태
    • 플러그인 설치전에는 msword가 종료된 상태여야한다.

  3. MSWORD를 켠 뒤, [참조]탭에서 확인할 수 있다.

각종 주의점들

  1. 필드는 오로지 1줄에 다 입력시키고, 많아지면 데이터 > 그룹을 활용하자
    • 칼럼1, 칼럼2, ..., 칼럼n, 빈 칼럼에서, 빈칼럼 앞까지 그룹을 해주면 빈 칼럼이 그룹명이 된다.





  1. 필드명은 영어로 하자. 통계 프로그램에서 인식 못할 수 도 있다.
    • 대쉬(-)보다는 언더바(_)가 좋다
  1. 범주 2개형의 경우, 성별이나 기저질환 유무 등 하나를 명시하고 다른하나는 0으로 나타내자. 그리고 데이터 유효성 검사를 활용하여, 각 셀에 입력될 값을 제한시키자.

    • 나이의 경우, 정수여야하며, 연령대 제한은 범위로 할 수 있다.







  • 성별(MALE)의 경우, 목록의 형태로 0,1만 오도록 제한할 수 있으며, 반드시 콤마로 나누어주어야한다. 그러면 드랍다운이 생기면서 선택하거나, 직접 0or1을 입력할 수 있다.




최근에 나는 데이터분석 준전문가(ASdP), SQL개발자(Sqld) 자격증 시험을 치룬 상태이다. 나에게 있어서 한의학이외의 관심은 이러한 분야들이다. 특히, 빅데이터에 관심이 있고, 통계가 서툰 상태에서 간단한 교양수준의 통계책을 읽고 싶었는데 마침 독후감 과제도 있었기 때문에 이 책을 감상하는 시간을 가졌다. 저자 니시우치 히로무는 최대한 실무에 활용할 수 있는 간단한 예제와 수식을 제외한 교양 이론으로 읽는 사람으로 하여금 통찰력을 가질 수 있게 해준 것 같다. 책에서 서술한 모든 통계정보를 A4 1장으로 요약하기엔 무리가 있으니, 기존 지식을 새롭게 이해한 내용 위주로 감상을 적으려고 한다.
먼저, 통찰의 통계학에 필요한 세가지 지식을 소개한다.

  1. 평균과 비율의 본질적 의미

    • 숫자로 표현되는 양적변수는 평균으로, 문자로 표현되는 질적변수는 비율로 정리하자.
  2. 데이터를 구간으로 이해하기

    • 평균과 비율로 정의되는 하나의 점만 생각하지말고, 그것들을 떠받히는 구간을 파악해서 어느정도까지 범위에 속해있는지를 항상 생각하자.
  3. 값을 어떻게 정리해야하는지에 대한 지혜

    • 평균과 비율을 이용하여 최종적으로 조절하고 싶은 결과인 아웃컴(성과지표)를 결정하고, 그 아웃컴에 영향을 미치거나 차이를 설명하는 요인인 설명변수를 설정해서 인과관계를 이용해 값을 정리하자.
      비지니스에서 가치있는 데이터분석은 최대화하거나 최소화해야하는 항목이 무엇인지 알아내는 것이라고 한다. 그것이 바로 아웃컴이고 이 아웃컴을 좌웅하는 원인 제공의 대상자인 설명변수는 큰 의미를 가지게 된다. 이러한 설명변수를 설정하는 데 있어서 우선순위도 필요하다.
    1. 인과관계가 너무 당연한 것은 배제할 것 -> 당연해서 시간 들일 필요가 없다
    2. 아웃컴에 영향을 명백히 미치더라도 조절이 가능한 변수여야한다
    3. 분석된 적이 별로 없는 변수여야한다.

이렇게 전반적인 통계에 대한 통찰을 설명해주고, 대표값들(평균, 표준편차) 및 가설검정, 회귀분석의 종류, 결과값이 없는 상태에서 특징을 가지고 분류할 수 있는 인자분석과 군집분석, 마지막에는 총정리를 해주었다. 또한 수학적으로 궁금한 부분은 따로 부록으로 싣어놓아서 필요하면 찾아볼 수 있게 해주었다.
특히나 기억에 남는 서술을 골라서 기술해보면, 제 1종 오류(귀무가설이 옳은데, 대립가설을 채택한 경우)를 덜렁이라 표현하며 제 2종 오류(귀무가설이 틀린데, 귀무가설을 기각하지 않은 경우)를 멍청이라고 표현하였다. 사실 덜렁이는 유의미한 차이가 없는데도 있다고 우기는 경우라고 할 수 있는데, 가설검정의 p-value값이 0.05이하가 나왔다면 덜렁이가 되지 않는다. 이러한 경계선을 유의수준이라고 하며, 제 1종 오류를 범하더라도 인정해주는 최대 한계선이기 때문이다. 즉, p-value 0.05는 덜렁이가 되지 확률 95%를 의미하는 것이기도 하다. 반대로 멍청이는 유의미한 차이가 있는데도 놓치는 경우인데, 제 1종 오류(덜렁이)보다 덜 중요한 문제로서, 멍청이가 되는 경우의 수는 적다고 한다. 다음으로, 회귀분석에 대한 서술이다. 그 중에서도 단일회귀분석에 대한 설명이 참 재미있었다. 만약 어떤 아웃컴을 설명하는 변수가 1개만 있다고 가정하였을 때, 중학수학으로도 회귀식을 표현할 수 있다. 그러나 그 이전에 설명변수를 x축에, 아웃컴을 y축에 놓고 점을 찍는 산포도를 관찰할 수 있는데, 이것은 경향성(상관관계)를 볼 수 있지만, 인과관계는 성립될 수 없다고 하였다. 왜냐하면, 설명변수는 1개라고 설정해놓았을지라도 미처 발견하지 못한 설명변수들이 여러개가 있을 수 있고, 그것들이 통제되지 않는 이상 인과관계는 성립하지 않기 때문이다.
이러한 산포도 위에는 회귀식을 그을 수 있는데 Y= b0 + aX로서, 중학교 직선의 식에 해당한다. 그리고 기울기 a를 해석하는 방법인 "a가 한단위 증가함에 따라 Y가 얼마나 증가하는지"를 관찰하면 된다. 이렇게 단일회귀분석은 간단한 것이었다. 이것 이외에도 다양한 분포와 가정들을 설명해놓았지만, 정보의 양이 너무 많아서 수식과 함께 공부하는 시간은 따로 가지려고 한다.
만약, 내가 한의사가 되고 논문을 쓰거나 사업을 하려면 통계는 필수덕목일 것이다. 요즘같은 빅데이터 시대에서 의학통계에 대한 전공수업이나 교양수업이 부실하다는 것은 너무 안타깝다. 개인적이라도 이러한 분야에 관심을 가지고, 교양으로서 통계를 친숙하게 만들 필요성을 절실하게 느끼는 본과 4학년이다.

로지스틱 회귀분석

지금까지 학습한 선형 회귀분석 단순/다중은 모두 종속변수Y가 연속형이었다.
로지스틱회귀분석종속변수가 범주형이면서 0 or 1인 경우 사용하는 회귀분석이다.

로지스틱 회귀분석을 설명하기 위해서는 먼저 로짓 변환오즈에 대해서 알아야한다.

왼쪽그림의 경우, Y가 0 or 1(사망/생존, 실패/성공, 불합격/합격)이라면 선형회귀로는 fitting하기가 힘들다. 그래서 곡선으로 fitting하기 위해 사용하는 것이 logistic함수 = 로짓변환이다.

Odds이고, Odds ratio비의 비율이다.
Odds의 해석은 확률에서 시작되어 실패에 비해 성공할 확률의 비를 의미하며 , odds = p / 1-p로 계산한다.
예를 들어, 게임에서 이길 확률이 1/5, 질 확률이 4/5이면, 게임에서 이길 odds 1/4이 되며, 계산된 값을 바탕으로 5번 중에, 4번 질 동안 1번 이긴다라고 해석한다.

Odds에 Log를 취한 것이 바로 로짓(Log p/1-p)이다.
예를 들어, 아래와 같은 교차표가 있다고 가정해보자.
Odds라는 것은 각 독립변수(drugA, drugB)에 대해 실패/성공에 대한 확률을 구한 뒤, 각각 구하는 것이며, 각 약에 대한 <생존Odds>가 구해지면 -> Odds ratio도 계산할 수 있다.
Odds ratio는 위험요인과 질병발생간의 연관성을 나타낼 때 사용 + 논문기재시 신뢰구간도 같이 제시해야한다.
예를 들어, 대조군(DrugB)와 실험군(DrugA)를 이용해서, 위험요인(DrugA)와 생존/사망의 연관성을 나타낼때, 각각의 Odds -> Odds ratio를 구한 뒤 제시한다.
오즈비 = 교차비 = 승산비 = 대응위험도 라는 표현도 쓴다.

  1. Odds(A)구하기

    • P(A, A에 대해 생존확률) = 20 / 52 = 0.38
    • 1-P(A, 사망확률) = 0.62
    • Odds(A) = 실패(사망)에 비해 생존(성공)할 확률의 비 = 0.38 / 0.62 = 0.61
    • A약 먹으면, 100명 사망할 동안, 61명 생존
  2. Odds(B)구하기

    • P(B, 생존/성공확률) = 24/66 = 0.63
    • 1-P(B, 사망/실패확률) = 0.37
    • Odds(B) = 0.63 / 0.37 = 1.7
    • B약 먹으면, 100명 사망할 동안, 170명 생존
  3. Odds ratio구한 뒤 해석하기

    • B에 대한 A의 Odds ratio = 0.61 / 1.7 = 0.36
    • 해석 : B에 비해 A일 때, 생존(성공)이 0.36배 = 64%가 생존율(성공률)이 떨어진다

로짓 변환

일반 선형회귀모형 glm : f(x) = b0 + b1X1 + b2X2 + .. bpXp

  • f(x)는 링크함수라고 한다. 이 f(x)자리에 특별한 함수들이 들어가서 구체적인 회귀분석을 만든다.
  • 다중선형 회귀분석 : Y = b0 + b1X1 + b2X2 + .. bpXp
  • 로지스틱 회귀분석 : ln( p/1-p ) = = b0 + b1X1 + b2X2 + .. bpXp

즉, 로지스틱 회귀분석은 일반회귀모형의 링크함수를 로짓으로 변형한 분석이다.

선형회귀분석은 Y값이 -무한대 < Y < +무한대의 연속형이다.
그러나 로지스틱 회귀분석의 경우

  • 먼저, p는 0 < p < 1
  • 0< 1-p < 1
  • p가 0에 가까울 경우 : p/1-p = 0/1 = 0
  • p가 1에 가까울 경우 : p/1-p = 1/0 = 무한대
  • p/1-p(오즈)의 범위 : 0 < Odds < 무한대
  • ln( p / 1-p )의 범위 : -무한대 < 로짓 < +무한대

역시.. 링크함수의 범위가 -무한대 < <+무한대 범위를 가지게 된다.

독립변수가 1개인 상황의 로지스틱 회귀분석의 식을 p(독립변수X가 성공할 확률)에 대해서 정리해보자.

그리고 이 식을 X에 대한 P의 그래프로 나타내면, P는 0과 1의 값을 가진다.

관측치x를 대입하면 P는 0 or 1로 fitting시킬 수 있게 된다.

로지스틱 회귀분석에서 통계량의 이해

ln(p/1-p) = b0 + b1X1 + b2X2 + ... + bpXp의 로지스틱회귀분석에서
X1을 위험인자라 가정하고, 아래에 대한 해석을 해보자.
오즈비(Odd ratio)에 대한 해석

  1. 위험인자 있을 때의 로짓1 = 회귀식1(X1=1대입)을 만든다.
  2. 위험인자 없을 때의 로짓2 = 회귀식2(X1=0대입)을 만든다.
    • 회귀식에서 b0, b2X2, ... bpXp는 모두 같고, X1만 대입한 식이 다를 것이다.
  3. 로짓1과 로짓2를 뺀다.
    • 우항은 b1*1 - b1*0 으로 인해 b1만 남는다.
  4. ln(p/1-p) - ln(p'/1-p') = b1
  5. ln( 위험인자있을때의 오즈 / 위험인자없을때의 오즈) = b1
    • ln( 오즈비 ) = b1 ==> b1을 ln(오즈비)로 생각하자.
  6. b1 = ln(오즈비)이므로, b1에 exp취한 값이 오즈비
    • R에서 실제 구해지는 결과는 b1이므로, exp를 취한 것(exp(b1))을 오즈비로 해석하면 된다.
    • exp(b1) = Odds ratio는 기준이 1이며 효과가 없다로 해석한다. 즉, 위험요인 0에 비해 위험요인1인 경우, 사망율 1배 -> 0% 사망율이 떨어진다. = 효과없다.
    • 1보다 클 경우risky로 해석한다. 즉, 위험요인 0에 비해 위험요인 1인 경우 사망률이 1.5배 -> 사망율이 50% 증가한다. risky는 분모오즈의 상황(base line, 대조군, X1=0)에 비해 분자오즈의 상황(위험인자, 실험군, X1=1)이 더 큰 상황
    • 1보다 작을 경우protective로 해석한다. 즉, 위험요인 0에 비해 위험요인 1인 경우 사망률이 0.5배 -> 사망률이 50% 감소한다.
  7. 유의사항1 : 단, 다른효과가 일정(동일)한 상태라는 가정 이 있어야 한다.
  8. 유의사항2 : 신뢰구간(Confidence Interval)가 같이 제시해야하며, 적은 값을 기준으로 소극적(보수적)으로 해석한다.
    • 예를 들어, OR(odds ratio)의 95% CI = 5.4 ( 3.3 ~ 7.8) => 3.3기준으로 최소 3.3배는 위험하다로 소극적 해석


예를 들어보기

종속변수survived(Y)는 0과 1로 구성되어 있고, 0 = 사망, 1 = 생존 이며, 이것을 예측한다고 가정해보자.
R에서 선형회귀분석을 한다면 lm()이며,
R에서 로지스틱 회귀분석을 한다면, glm() + family = 'binomial'인자를 사용해야한다.

  • model <- glm( survived ~ . , family = 'binomial'(link = 'logit'), data= train)

model을 summary()했을 때, Estimated로 나온 것은 각 변수들의 회귀계수이다. 여기에 exp()를 취해야지 Odds ratio가 나온다.

Coefficiets(회귀계수)에는 절편(Intercept) + 각 변수들의 기울기와 표준편차 + wald's z-statistics + p-value 가 나온다.

만약 범주형 변수 sex의 회귀계수를 보면, 범주형이므로 male은 기준이 되어 male에 비해 female 한 단위 증가시 Y의 변화량이 나타날 것이다.
그 추정치는 -2.6609987이다. 이것은 bp(회귀계수, 기울기)이므로, exp()를 취해야지 Odds ratio가 나온다.
exp(-2.66) = 0.069 => *남자(ref, Xp=0)에 비해 female(Xp=1, female)의 생존율 이 0.069배 감소한다 *. 왜냐하면 회귀계수가 음수이므로 반대로 해석한다. = female의 생존율은 93.1% 증가한다.

직접적으로 바로 exp( model$coefficients)오즈비(Odds ratio)만 관찰할 수 있다.
Age의 오즈비가 0.959가 나왔다? => age가 한단위 증가할수록 생존율이 0.959배 증가한다. -> 0.041 감소 -> 생존률이 4.1% 감소한다.

직접적으로 바로 신뢰구간을 포함한 오즈비를 보는 방법은 exp( confint(model) )를 통해 확인할 수 있다. 2.5% cutoff와 97.5% cut off를 제시한다.

회귀분석 이전에 이상치 + 잔차를 통한 가정사항을 확인해야한다.

회귀모델의 이상치 확인

  1. 이상치 확인1 - outlierTest()
    개별적으로 검사하는 것이 아니라 회귀모델 자체를 넣어주면 outlierTest( 모델 )함수로 확인할 수 있다. 잔차가 2배이상으로 크거나 2배이하로 작은 값을 이상치로 detect한다
    결과에서, 해당변수의 P-value가 0.05이하로 나오면 이상치라 판단하고, 제거한 모델을 만들고 다시 회귀분석한다. ( H0는 이상치가 아니다 인가보다)

    x축 = 지레점(Leverage) => 높으면 지렛점
    y축 = 표준 잔차 => 높으면 이상치
    점선 => i번째 변수가 빠진상태에서, 회귀계수에 얼마나 큰 영향을 주는지에 대한 영향점영향관측치)

  2. 이상치확인2 - influencePlot()
    x축이 크면 지레점
    y축이 크면 이상치
    원의크기가 크면 Cook's Distance에 대한 비율로서 정해지는 영향점



잔차진단을 통한 가정사항 확인

잔차 = 관측치(Y, observed) - 추정치(estimated)

제대로 된 모델에서 나온 잔차

  1. 정규분포를 따르고
  2. 분산이 일정하고
  3. 특별판 추세(패턴)이 없어야한다.

=> 잔차에 추세(패턴)을 보인다면, 회귀모형에 포함되어야할 정보가 빠졌다는 것이다.



회귀분석의 이론적 배경

회귀모형( Y = b0 + b1X + e)의 가정

  1. X는 비확률변수이며, 주어진 어떤 값
  2. e는 서로 독립이며, 정규분포 N(0,o^2)을 따른다.
    • 특히 2번의 가정이 맞지 않는 경우, 단순 회귀분석은 의미가 없다.

e(오차, 엡실론)이 서로 독립이 아닌 경우

  • 예를들어, 시계열 데이터는 각 데이터마다 관계성이 있으므로 오차가 독립이 아니다. 이럴때는 Durbin Watson 의 D통계량을 체크해야한다고 한다.

이분산(Heteroscedasticity)의 문제

  • 잔차가 등분산이 아닌 경우는 OLS(최소제곱법, 오차의제곱의 합이 최소가 되도록하는 회귀식 구하기)가 아닌 WLS(Weighted Least Square), GLS(Generalized Least Square)를 사용해야한다고 한다.

이산적(Discrete) 종속변수인 경우

  • 종속변수Y가 연속형(숫자형)이 아닌 경우에는, 정규분포를 따르지 않으므로 로짓/프로빗 모형을 적용해야한다고 한다.


결과적으로 잔차( resid( 모델 )) 에서 체크해야할 사항 3가지

  1. 정규성 검정 : 만약 만족하지 않는 경우, Log, Root를 취해서 정규분포를 취하도록 만든다.
    1) resid( lm() )을 shapiro.test()에 넣으면 된다. H0 : 정규분포를 따른다.
    2) normal Q-Q plot에서 45도 line에 있다면 정규분포를 따른다.
    3) car패키지의 qqPlot()으로 정규성을 평가하면 신뢰구간(CI)까지 표시해준다.

  2. 등분산성 검정 : 잔차가 동일하 분산을 가지지 않는 경우, 가중치를 고려한 WLS, 혹은 GLS를 이용해서 회귀식(선)을 만든다.

    • Fitted Value와 Residual값(or 표준화된 Residual값)plot으로 그려서 평행하게 분포(특별한 패턴X)되면 등분산성을 만족한다. 만약 cone모양으로 점들이 퍼진다면, 등분산성을 만족하지 않는 것이다.
  3. 독립성 검정 : 잔차가 독립적이지 않다면, 시계열분석밖이다.

    • 잔차의 독립성은 durbinWatsonTest( fit )를 통해한다.
    • dw test의 결과인 D-W statistic0 ~ 4의 값을 가지며, 0은 양의 상관관계를 가져 독립x, 4는 음의 상관관계를 가져 독립x, 2가 독립적인 것을 의미한다. 검정결과가 2에 가깝다면 독립인 것이다.
    • p-values는 rho(자기상관관계)에 대한 것이다. 만약 p값이 0.05보다 작다면, 자기상관관계에 있다.
  4. 한번에 검정하기 by gvlma패키지의 gvlma( 모델 )

    • gvlma(모델) 이후 summary()를 하면 Global stat값이 나온다. 이때는 p-value가 0.05보다 커서 H0가 기각 안되어야 H0 : 3+1(비선형회귀가능성)가정이 다 만족하는 것
    • 만약 summary( gvlma(모델) )의 결과 Global stat의 p-value가 0.05보다 작다면 H0가 기각되어 적어도 1가지 이상이 잔차의 가정을 만족하지 x => 1,2,3 다 하나씩 봐야한다. 하지만 잔차의 진단에서는 p-value를 너무 strict하게 보지 않아도 된다.

다중 공선성의 해결방안

어느정도의 X들(설명변수들)끼리의 영향은 존재할 수 밖에 없으나, 그 교집합이 심하게 겹칠 경우, 다중 공선성이 발생한다. 이 때는 적절한 변수를 선택(제거)를 하는 변수 선택법이 필요하다.

사실 X의 개수(컬럼수, 변수수)가 증가할수록 예측error가 감소하는 구간이 있다.
하지만, 적절한 복잡성을 넘어선 모델 복잡성(다중공선성으로 인해)은 과적합 = 오버피팅을 일으켜 오류가 증가하게 된다. 비록, 현재의 데이터(훈련데이터) 오류가 계속 줄어든다고 하더라도, 모르는 집단(실제데이터, 테스트 데이터)의 오류가 줄어들다가 다시 증가하는 순간이 생긴다. 이 때, 적절한 변수만 선택하여 과적합, 오버피팅을 발생시키지 않게 해야한다. 의미있는 변수들만 선택해야한다.

변수 선택법 3가지

  1. 전진 선택법(Forward selection) : X 0개부터 시작하여, 가장 유의한 변수들부터 하나씩 추가하면서, 매번 모형의 유의성을 판단한다

  2. 후진 제거법(Backward elimination) : 모든 X를 넣어놓고, 하나씩 제거해간다.

  3. 단계적 방법(stepwise selection) : 0개부터 추가했다가 뺐다가 왔다갔다함.

    • r에서 step( 모델, direction = "backward")식으로 step()함수를 이용한다. 그러면 결과창에 AIC가 찍히므로 확인한다.
  4. 의학분야 : 임상적인 의미를 가진 변수, 기존 연구에서 검증된 변수를 유의성과 관계없이 넣을 수 있다. 실제 임상에서는 단변수 -> 다변수로 체크하면서 변수를 추가하는 경우가 많다.
    예를들어, Y에 X1, X2, X3, X4, X5가 있다고 치면,
    단변수 체크 : Y ~ X1, Y ~ X2... 다 체크해서 유의성 높은 것 뽑느다. 0.05는 너무 STRICT한 기준으로 보고 0.1 ~ 0.2까지 허용해준다.
    다변수 체크 : 단변수에서 유의성이 높은 변수들을 가지고 유의성을 체크한다 Y ~ X2 + X4+ X5

### 모형 선택 기준 1. `결정계수(R^2) / 수정된 결정계수(adj R^2)` : best모델의 설명력(나머지는 설명안되는 error) / 여러모델 비교시 사용(높을수록 좋다) 2. `AIC(Akaike Information Criterion)` : -2logL + 2K (낮을수록 좋다) 3. `BIC(Baysian Information Criterion)` : -2logL + Klogn(낮을수록 좋다)

### AIC + ANOVA를 이용한 회귀모델 비교 1. `AIC(모델1, 모델1)` : 2모델 중 AIC가 낮은 모델을 선택한다. 2. `summary(모델1) (모델2)$r.squared` : 각 모델마다 설명력을 확인한다(선택에는 영향x) 3. `summary(모델1) (모델2)$adj.r.squared` : 각각의 모델에서 나온 수정된 결정계수를 비교하여 높은 모델을 선택한다. 3. `anova( 모델1, 모델2)` : `H0 : 모형1 = 모형2`가 기각되는지 `F pr(>F)`를 확인하여, `H0`가 기각되지 않고 같은 모형으로 판단된다면 **`변수가 작은 모델을 선택한다`**

### 가능한 모든 경우의수의 모델 조사 by leaps패키지 많이 쓰진 않는다고 한다. regsubsets()함수를 이용하여 nbest = size당 n개의 최적 모형을 저장하여, plot으로 그리면 된다. y축에는 adj R^2가 나오는데, 높은 모델을 사용하면 될 것이다.

회귀모형(lm())으로 만든 결과(summary())를 해석

lm( Y ~ X , data = )를 summary()로 요약하면

  1. Residual :
    Residual이라는 것이 나온다. 회귀모형으로 예측한 Y 와 실제 데이터의 Y의 차이를 의미하는 것이며, 찌꺼기로 취급해야한다.
    이 찌꺼기에는 패턴이나 추세가 있으면 안된다. 만약, plot으로 그린 결과에 어떠한 패턴이 있다면, 중요한 변수를 빠트려서 잔차에 남아있는 것이므로 찾아서 반영해줘야한다.

  2. coefficient :
    절편(b0)와 각 X의 기울기가 나온다. (Intercept)의 Estimate절편을 추정 한 값이며, 각 변수들의 Estimate기울기를 추정한 값이다.
    또한, 추정한 해당 회귀계수(절편, 기울기)에 대한 유의성을 단일표본t-test, mu = 0 (H0 : 해당 회귀계수 = 0)에 대한 p-value가 Pr(>|t|)로 표시된다.

    • 0 ~ 0.001 ***, 0.001 ~ 0.01 **, 0.01 ~ 0.05 *, 0.05 ~ 1 아무표시없음
    • 0.05이하로 H0가 기각되어야지, 해당 X의 회귀계수 0이 아니다 = 의미있는 회귀계수
    • 0.05보다 커서 H0가 기각안되면, 해당 X의 회귀계수는 0이다! -> 제외시켜야한다.
  3. multiple / Adjusted R-squared :
    먼저, ANOVA에서 나왔던 개념인 변동에 대한 이해가 필요하다. ANOVA에서는 총변동 = 그룹간변동(effect) + 그룹내변동(error)였다.
    회귀분석에서의 총 변동 = 회귀식이 설명할 수 있는 변동(effect) + 회귀식으로 설명할 수 없는 변동(error)로 구성되어있다.
    1) SST(결과변수의 총변동) = SSR + SSE
    2) SSR(Sum of Square of Regression) : 회귀식으로 설명되는 변동
    3) SSE(Sum of Square of Error) : 회귀식으로 설명하고 남은, 설명되지 않은 변동
    4) R-square(R^2) : SSR / SST (전체변동 중 회귀식으로 설명되는 변동) = 0 ~ 1

    그림으로 나타내면, 1 - SSE/SST정도인가보다.

3-1. Multiple R-squared :
결정계수 R^2(혹은 설명력 R^2)으로서, 앞에 말한 SSR/SST를 의미한다. 추정된 회귀식이 변동을 얼마나 잘 설명하는가로 해석된다.
만약, 값이 1이 나왔다면, 실제 관측값들이 회귀선상에 정확히 일치함을 의미한다.
만약, 값이 0.65 => 35%는 회귀식으로 설명할 수 없는 error임을 의미한다.
하지만, 치명적인 단점으로서 X의 개수가 증가할수록(아무상관없는 X를 추가하더라도) Multiple R-square(회귀식의 모형설명력 결정계수)도 무조건 증가하는 단점

3-2. Adjusted R-squared :
Multiple R-squared의 단점을 보완하기 위해, 보정된 R-squared(회귀식의 모형설명력, 결정계수)
변수(X)가 많아진 것을 분모에 반영하여, 변수의 개수가 고려된 R-squared(모형설명력, 결정계수)이다.

분모에 k는 변수의 개수를 의미한다.

3-3. 정리
1. Multiple R-square : BEST모델(Best회귀식)의 설명력(결정계수)
2. Adjusted R-square : 여러개의 모델을 만들어 놓은 상태에서, 좋은 모델을 찾기위해 비교시 사용되는 설명력

  1. F-statistic :

H0 : 모든 회귀계수에 대해 b1(=b2=b3=..bk)= 0을 가정한 뒤, H0를 기각시켜 회귀계수 0이 아닌 것이 적어도 하나 존재한다회귀모형의 유의성F-test한 것이다. F-검정통계량 계산은 ANOVA에서 처럼, F= SSR의평균 / SSE의평균으로 검정한다.

  • cf) 만약, 단순회귀(X1개 -> H0 : b1 = 0)일 시, coefficient의 t-test(단일표본t-test)나 모형의 F-test나 동일하다.
  • cf) F = t^2 (??) : X가 1개일때만?? 아니면 원래 분포자체가??

범주형의 X에 대한 회귀계수의 해석

앞서 숫자형 변수 X에 대해서 회귀계수(기울기)는 X가 한단위 증가시 Y의 증가량으로 해석했다.
범주형 변수 X의 경우, 범주안의 집단들 중에 1집단은 reference가 되고, 나머지 집단들이 나눠서 다른 변수처럼 표시되며, ref집단에 비해 나머지집단들의 X 한단위 증가시 Y의 증가량으로 해석된다.
예를 들어, mtcars데이터의 cyl(범주형)변수에는 4, 6, 8의 집단이 있다. 이 변수를 X(설명변수, 독립변수)로 넣어서 회귀모델을 만들면, 4집단은 표시가 안되면서 (cyl)6, (cyl)8에 대한 회귀계수가 생긴다. 해석은 (cyl)4에 비해, (cyl)6, (cyl)8집단이 Y에 어떻게 영향을 미치는지로 해석해야한다.

  • R에서 각 범주의 집단확인은 table(칼럼인덱싱)으로 하면된다. 마치 교차표 만들듯이

내부에서 작동하는 방식은 범주형을 -> 숫자형처럼 만드는 가변수Dummy 변수를 만들어서 처리한다. 집단이 3개인 경우, refence집단0,0으로 만들어놓고 나머지 2집단은 1,00,1이 될 것이다. 한 집단은 0으로 구성되어 기준이 되니, 집단수-1개의 자리수를 가지는 0과 1로 구성된다.

  • reference(기준)변수는 보통 빈도가 많거나 일반적인 집단을 기준집단으로 선택한다.

예를 들어, 약물의 종류가 P(lacebo), A(treatment A), B(treatment B)가 있다고 치자. 이때는 일반적인 집단인 P가 기준집단으로서 0으로만 구성될 것이다.

  • R에서 범주칼럼인덱싱 = relevel( factor(범주칼럼인덱싱) , ref = )로 기준집단을 직접 선택하여 더미변수화 시킬 수 있다.

표준화 회귀계수 for 회귀계수의 상대적 중요성

X-m/ sd(X)표준화시키는 목적은 단위를 고려한 비교를 위해서이다.
회귀분석전 모든 숫자형 변수들을 표준화시킨 후 회귀분석lm()을 돌리면, 표준화된 회귀계수가 얻어지며, 이것들을 plotting해서 단위를 배제시킨 회귀계수의 상대적 중요성을 확인할 수 있다. 이것은 가장 쉬운 변수 선택 방법이다.

  • SPSS로 회귀분석을 돌리면 표준화(B)+비표준화 회귀계수(b,베타)가 모두 출력된다. 그때는 표준화 회귀계수를 보고 상대적 중요성으로 해석하자!

신뢰구간은 보통 95%로 제시한다. 5%의 유의확률은 양측검정시 왼쪽 2.5%이하 <--> 오른쪽 97.5%이상의 임계치를 제시해준다. 추정한 회귀계수(기울기)가 2.5% ~ 97.5%는 확신을 가지고 있다는 것을 의미한다.
해석은 해당변수X가 1% 변화시, 종속변수Y는 추정회귀계수(2.5% ~ 97.5%)만큼 변화한다고 <<확신을 가지고 이야기할 수 있다>>

  • R에서는 confint()에 lm()결과를 넣어주면, 각 회귀계수별 신뢰구간을 제시해준다.

그러나, 이러한 해석을 위해서는 회귀분석의 가정 몇가지를 만족시켜야한다.

5-2. 회귀분석(모델링)

2019. 2. 25. 23:34

회귀분석과 로지스틱 회귀분석

  1. 회귀분석(Regression) : 관심있는 Y에 대해, X가 어떻게? 얼마나? Y에 영향을 주는지 알아보는 분석방법
  2. 로지스틱 회귀분석(Logistic Regression) : 관심있는 Y가 0 or 1일때, X가 어떻게? 얼마나? Y에 영향을 주는지 분석
    • 보통 ML과 DL의 성능에 대해서 reference로서 로지스틱 회귀분석을 먼저 제시하기도 한다.

회귀분석의 Workflow

  1. 모델링
    1) 단순 선형회귀 : Y에 영향을 주는 X가 1개
    2) 다중or중 선형회귀 : Y에 영향을 주는 X가 2개 이상

    • 다중공선성의 문제발생 : 2개 이상의 X가 서로 영향을 끼침.
  2. 모형 해석
    1) 계수 해석
    2) 표준화 계수b(베타)
    3) 계수의 유의성

  3. 모형 설명력
    1) 모형의 유의성
    2) 모델 설명력

  4. 모형 선택
    1) 변수선택기법 : 모형비교, 측도로서 AIC와 BIC

  5. 예측
    1) 잔차(residual) : 실제값 - 예측값 -> 찌꺼기인데 추세/패턴이 있는지 가정체크후 있다면 모델에 변수로 추가

1. 모델링

수학적 모형(오차없이 y=ax+b)과 다르게 통계적 모형으로서, 같은 input(x들)이라도 다른 y를 가지는 경우를 설명하기 위해 회귀방정식에는 항상 오차(엡실론)를 수반하는 형태를 가진다.

이 때, b0 = 절편, b1 = 기울기, e=오차(통계적 모형에서만 가짐)를 의미한다.
기울기 b1x가 1단위 증가할 때, y의 증가량으로 해석해야한다.
오차 e는 확률변수이며, 통계적 모형에만 적용되는 변수이다.

선형회귀는 독립변수와 종속변수의 1차 방정식을 전제로 한 선형성의 관계를 전제하여 분석한다. 그리고 데이터를 이용하여, 오차(e)를 최소화하는 방식으로 회귀방정식의 회귀계수와 회귀선을 구해야한다.
회귀선은 오차를 최소로하는 회귀방정식의 회귀계수들(b0, b1, b2 ..)을 다 구하여, 그래프상 가장 잘 설명하는 직선 하나를 찾은 것이다.

사실 회귀분석은 영국의 화학자 프랜스시 골턴이 발견한,
키 큰 부모 -> 자식들의 키는 점점 더 커지는 것이 아니라 평균으로 회귀하는 경향을 보고 만들어낸 개념이다.
독립변수 : 아버지의 키
종속변수 : 아들의 키
의 단순 선형회귀 문제로서, 키작은부모 -> 키큰자식(평균으로 가게끔), 키큰부모-> 키작은자식(평균으로가게끔)평균을 유지하려는 현상을 발견한 개념이 회귀인 것이다.

그리고 회귀분석에서 아버지의 키의 X가 1개만 있는 것이 아니라 여러개의 X들이 존재하는 것을 알 수 있다.
즉, 회귀분석이란 변수들의 관련성 규명을 위해, 어떤 수리적 모형을 가정, 측정된 데이터로 모형의 변수들을 추정하는 통계적인 방법이다.
Y에 영향을 주는 독립변수들 X에 어떤 것이 있고, 얼마나 영향을 주는지를 알아낸 뒤, 다음 Y값을 예측하는데 사용된다. 사실 예측력은 떨어지나 설명이 가능해서 자주 사용하는 분석방법이라고 한다.

예를 들어,
광고비용(X)가 매출액(Y)에 얼마나 영향을 줄까? -> 오차를 최소로하는 회귀선의 기울기로 해석한다. 광고비용 1단위 증가시 매출액 얼마나 증가했나?
공부(X)가 성적향상(Y)에 얼마나 영향을 줄까? -> 회귀선의 기울기로 해석
즉, X가 어떤 것들이 있는지 안 상태로서 선형회귀 = 회귀방정식+오차를 구한 뒤, 오차를 최소로하는 회귀선에서 기울기로 영향을 해석하면 된다.


위의 데이터에서는, X가 1개인 단순 선형회귀이며,
같은 speed(X)라도 다른 dist(Y)를 가지는데, 이는 오차(e)가 존재하기 때문이다.
dist = b0 + b1*speed + e 로 회귀식을 놓고, e를 최소화하는 b0, b1를 구해야할 것이다.

회귀식(Y = b0 + b1X + e)를 구하는 방법을 앞에서 계쏙 오차를 최소화한다고 했다.
사실, 구체적으로는 각 데이터의 실제값y와 회귀식의 값 yi의 차이인 오차(e)에 대해, 모든 오차들의 제곱의 합이 최소가 되는 회귀계수들을 구하는 방법인 최소제곱법(OLS, Ordinary Least Square)를 이용한다.
가상의 단순선형 회귀식 Y = b0 + b1X +e 에서 오차e를 좌변에 나머지를 우변에 이항한다.
e = Y - b0 - b1X
그리고, 각 데이터들을 입력하여 나오는 모든 e값들을 제곱해서 더한다.
오차의 제곱의 합을 S라 두면, 아래와 같은식이 나온다.

  1. 오차제곱의합(S, SSE)를 절편(a, b0)에 대해 편미분한다.
  2. 오차제곱의합(S, SSE)를 기울기(b, b1)에 대해서 편미분한다.
  3. 두 연립방정식을 이용하여, 절편(a)와 기울기(b)를 구할 수 있다.

이러한 방식으로 회귀계수(b0, b1-절편과 기울기)를 구하면, 회귀식이 완성된다.
이 때, 더이상 오차는 없어지고 하나의 회귀식이 만들어지는 것이다.

  • R에서는 lm()함수를 사용하면, 안에서 자동으로 b0, b1을 구해준다.
  • coefficients가 회귀계수 = Intercept(절편) + df$칼럼명(해당X앞에 곱해진 기울기)
  • lm()의 결과를 그대로 산점도 위에 abline()으로 그려주면, 그래프 상으로 회귀식을 그려줄 수 있다.
  • ggplot으로 그릴 때는, 산점도(geom_point()) + geom_smooth(method="lm")으로 자동으로 그려준다. 이 때, 신뢰구간까지 같이 그려지는데, n수가 적은 곳은 conf.interval이 넓고, n수가 많은 곳은 conf.interval이 좁다고 해석하면 된다. 좁을 수록 명사수.

사실 lm()에는 회귀계수( 절편, 해당X의 기울기들)만 있는 것이 아니라, 결과물로서 모델이 생성된다.

  • lm()결과를 reg_model변수에 저장하고, 새로운 데이터를 만들어서 predict( 모델변수명, 새로운데이터) 형식으로 predict()함수를 사용한다.


다중 공선성

이제 다중 회귀분석에 대해 다시 알아보자. Y에 영향을 미치는 X가 여러개 이므로 X to Y뿐만 아니라 X끼리 영향을 미치는 문제인 다중 공선성(Multicollinearity)의 문제가 발생한다. 즉, 독립변수(X)들간에 강한 상관관계가 나타나는 문제이다. 이것들을 어떻게 파악하고 어떻게 해결해야할까?

X들끼리의 강한 상관관계로 정보가 많아져서 Y를 예측하는데 더 좋다고 생각할 수 도 있으나, 최소제곱 추정지(OLS)의 공식상, 설명변수(X)들간의 상관관계가 높으면 OLS추정치의 분산이 불안정해진다고 한다.
다른 관점으로서 각 X들 앞에 곱해져있는 회귀계수들의 해석으로서 생각해 볼 수 있다.
X앞에 곱해져있던 b1은 가정이 있다.
다른 변수(X)들이 일정한 상황에서, 해당 X가 한단위 증가할 때 Y의 증가량이다. 즉, 다른 X들은 일정해서, 해당 X에 영향을 주지 않는 상태에서의 계수를 의미한다. 회귀계수의 정의자체가 서로 영향을 주지 않는 상태가 기본 전제인 것이다.
예를 들어, Yi = b0 + b1*X1 + b2*X2 인 상황에서, b1(기울기, X1의 회귀계수)의 해석은, X2가 일정한 상황에서 정해진 계수로서 해석하는 것이다.

만약, X1 : 태어난 날짜, X2 : 나이 처럼 X1과 X2가 완전한 상관관계(하나가 일정할 때, 를 가정이 불가능)에 있는 경우, 다중공선성 문제가 크게 발생하므로 -> 회귀계수를 정확하게 추정하기 어렵다.
또한, X1 : 상품의 내구도, X2 : 불량률 처럼 X1과 X2가 negative 상관관계가 너무 높은 경우도 다중공선성이 발생한다(하나를 일정하게 둘 수가 없다)

이러한 X들간의 상관관계를 확인하는 방법은 2가지 정도로 압축된다.

  1. X들간의 산점도그리기 : 직선성으로서 강한 상관관계를 보이면, 동시에 넣으면 안된다
  2. X들간의 VIF(Variation Inflation Factor) 혹은 공차한계(Tolerance) 계산

참고사이트 : https://datascienceschool.net/view-notebook/36176e580d124612a376cf29872cd2f0/

VIF의 역수가 Tolerance = 1 - R^2이다.

이 때, Ri^2i번째 X를 종속변수(Y)로 취급하여 회귀모형을 만든 후 Xi에 대한 나머지 X들의 모형 설명력이다.
예를 들어, Y ~ X1+ X2+ X3 모델에서, X1에 대한 VIF를 구하기 위해서 R1^2을 구해야한다. 구하는 방식은 X1 ~ X2 + X3의 모델을 만든 뒤 R1^2이 계산된다. 그 값이 0.9가 나왔다면, R1^2 0.9 = X2와 X3로 X1을 90% 설명할 수 있다.
그러면 VIF는 1 / ( 1- 0.90) 이 되고, Tolerance는 1-0.90이 된다.
이 때, VIF가 4이상 혹은 6이상 넘어가면, 해당 X는 다중공선성을 의심할 수 있다. VIF가 10(R^2 - 0.9)이상이면 보통 다중공선성 문제가 있다고 한다.

같은 이치로서, X1이 만약 다중공선성이 의심되면, 그것과 짝을 이루는 Xi가 반드시 존재할 것이다. 만약, X1과 강한 상관관계의 다중공선성을 가진 X2가 있다면, 둘중에 하나는 제거해야한다.

  • R에서는 car패키지안의 vif()함수를 이용해서 각 X에 대한 VIF를 나열해준다.

  • vif( model ) > 4라는 마스크를 만들고, True인 것을 찾아내면 된다.

  • 만약, 그러한 X가 존재한다면 최소 2개가 있을 것이다.(서로 연관)-> 1개만 선택해야한다.



다중회귀분석

Y에 영향을 끼치는 X가 여러개 있을 때는,linear combination이라는 선형결합 or 가중합으로 연결된다. 이것은 앞서, 1차 방정식으로 전제가 된다고 했다. 또한 다른 변수들이 일정할 때를 가정한 뒤, 해당 변수의 가중치(회귀계수)를 곱하는 보정이 들어간 것이다.
여기까지 만들어지는 것이 다중회귀 모형이다. 여기에는 e(엡실론)이 포함되어있다.

다음으로, 측정된 데이터들을 이용하여 OLS를 통해 오차의 제곱의 합을 최소로 하는 다중 회귀식을 구한다. e가 사라지고, 하나의 회귀식이 만들어진다.

하지만, 측정된 데이터 = 표본을 통해 만들어진 것이니, 결국에는 모집단을 추정하는 추정 다중회귀식이 만들어진다.

R에서는 lm( Y ~ X, data=)로 만들어진 단순선형회귀 모형(모델)에 대해, X를 추가해나가면서 다중선형회귀 모델이 되는 과정에서 자동으로 내부에서 회귀계수(기울기)를 보정하여 각각의 보정된 회귀계수가 summary()의 결과로 나타난다. 이 때 $coef로서, 회귀계수만 인덱싱해서 살펴보면 단순회귀와 다중회귀일때, 같은X에 대해서 다른 회귀계수가 나오는 것을 확인할 수 있다.

  • summary( lm( BMI ~ TC, data = acs) )$coef
  • summary( lm( BMI ~ TC + age + sex , data = acs) )$coef
  • TC의 Estimate=추정된 회귀계수값을 확인해보면, 변수를 추가할 때마다 회귀계수가 달라지며, 해석은 TC 한단위가 증가시 BMI의 증가량이 어떻게 변하는지 보면 된다.
  • 같이 등장하는 Pr( >|t|)는회귀계수에 대한 p-value이며,
    맨 마지막줄의 F-statistic의 p-value모형 자체의 p-value이다.
  • 즉, lm( Y ~ X의 선형모델 )은 OLS로 회귀모형을 생성
  • summary( lm())은 회귀계수 및 통계수치를 확인할 수 있다.
  • 각 회귀계수는 해당변수에 대해 다른 변수들이 일정할 때~의 가정이 있는 상태에서, 해당변수가 한 단위 증가시 종속변수가 얼마나 증가.감소하는지 로 해석하면 된다.

Tidyvers 소개

ggplot의 개발자인 Hadley Wickam의 Tidyverse패키지안에
dplyr, tidyr, tidytext, ggplot2, readr들이 다 포함되어있다.

  1. tidyverse : 포함된 패키지들 모두에서 pipe operator를 사용할 수 있다.

    • %>%를 사용하여 가독성 및 naming이 필요없어진다. 마치 ggplot의 + 연산자와 똑같음
  2. readr : read_csv()을 사용할 수 있어
    1) read.csv()보다 더 빠르다( string을 factor로 자동변환이 없음 )
    2) data.frame대신 tibble을 사용하여 최종적으로는 4~5배 빠르다고 한다.
    3) cf) data.table 패키지는 100배 이상 빠르게 사용할 수 있다고 한다.


Dplyr

가장 큰 특징 2가지는 pipe linegroup별 계산인데, 세부적으로 아래와 같다.

  1. select() : 쉽게 필요한 칼럼만 인덱싱 or 새로운칼럼명 = 으로 칼럼명 변경

    • select( data, 칼럼명)
    • select( data, c()없이 여러개 칼럼명 연결, 1, 2, 3)
    • select( data, starts_wwith("칼럼명 시작단어 선택")
    • select( data, ends_wwith("칼럼명 끝단어 선택")
    • select( data, contains("포함 단어 선택")
    • select( data, 새로운칼럼명 = 뽑아낸 칼럼 )
  2. filter() : 쉽게 필요한 행(특정범주)만 인덱싱

    • filter( data, 칼럼명 조건 )
    • filter( data, 칼럼명 조건 & 칼럼명2 조건2 )
  3. arrange() : 쉽게 오름/내림차순 정렬

    • arrange( data, 오름차순 정렬기준칼럼명 )
    • arrange( data, desc(오름차순 정렬기준칼럼명) )
  4. mutate() : 쉽게 새로운 칼럼 생성 by 기존칼럼를 변형

    • mutate( data, 새로운칼럼명 = 기존칼럼명 + 수식 )
  5. group_by() & Summarise() : 그룹별 통계량

    • filter()로 필요행(특정범주)인덱싱 후 summarise( filtedData, mean(칼럼명, na.rm=T)
    • filter() -> group_by( filterdData, 기준칼럼명 ) -> summarise( groupedData, 통계함수 ) 그룹별 통계
  6. %>% : 파이프라인

    • data %>% filter( 칼럼명과 조건) %>% group_by( 기준칼럼명 ) %>% summarise( 집계함수))
    • 각 함수에서 data부분만 생략됨
  7. 그외 count, sample


Tidyr으로 데이터 reshape

  1. spread() : to wide
  2. gather() : to long

R markdown

Untitled

Dplyr

  • dplyr has 5 verbs which make up the majority of the data manipulation tasks we perform.
  • 1.Select: used to select one or more columns;
  • 2.Filter: used to select some rows based on specific criteria;
  • 3.Arrange: used to sort data based on one or more columns in ascending or descending order;
  • 4.Mutate: used to add new columns to our data;
  • 5.Summarise: used to create chunks(묶음) from our data.

파이프라인 및 전체 연습

library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(moonBook)

#### 파이프라인 없이 해보기
# 1. filter()로 행인덱싱(특정범주 가져오기)
acs2 <- filter( acs, sex == "Female" & age > 40)
# 2. arrange()로 내림차순 정렬
ordered_acs2 <- arrange( acs2, desc(height))
# 3. group_by()로 그룹화 -> tibble로 변형 & 그룹화를 직접적으로 하진 않음
obesity_to_acs2 <- group_by(ordered_acs2, obesity)
head(obesity_to_acs2) # tibble상에서는 그룹화 안되어있음. 그룹별 집계가 가능한 상태만 Groups로 표시
## # A tibble: 6 x 17
## # Groups:   obesity [2]
##     age sex   cardiogenicShock entry Dx       EF height weight   BMI
##   <int> <chr> <chr>            <chr> <chr> <dbl>  <dbl>  <dbl> <dbl>
## 1    67 Fema~ No               Radi~ Unst~  57.2    168     52  18.4
## 2    89 Fema~ No               Femo~ STEMI  21.8    165     50  18.4
## 3    64 Fema~ No               Radi~ Unst~  57.7    165     74  27.2
## 4    53 Fema~ No               Radi~ Unst~  58.4    165     75  27.5
## 5    64 Fema~ No               Radi~ STEMI  59      165     71  26.1
## 6    48 Fema~ No               Radi~ STEMI  61.5    165     58  21.3
## # ... with 8 more variables: obesity <chr>, TC <dbl>, LDLC <int>,
## #   HDLC <int>, TG <int>, DM <chr>, HBP <chr>, smoking <chr>
# 4. summarise()로 그룹별 집계
obesity_to_BMI_mean <- summarise(obesity_to_acs2, mean( BMI, na.rm = T))


#### 파이프라인을 통해 naming걱정없이 바로 해보기
acs %>%
  filter( sex == "Female" & age > 40) %>%
  arrange( desc(height))  %>%
  group_by( obesity ) %>%
  summarise( mean( BMI, na.rm = T))
## # A tibble: 2 x 2
##   obesity `mean(BMI, na.rm = T)`
##   <chr>                    <dbl>
## 1 No                        22.2
## 2 Yes                       27.7

select() : 필료한 칼럼인덱싱 & 칼럼명 변경도 가느

library(dplyr)
head(iris) # 내장 데이터
##   Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 1          5.1         3.5          1.4         0.2  setosa
## 2          4.9         3.0          1.4         0.2  setosa
## 3          4.7         3.2          1.3         0.2  setosa
## 4          4.6         3.1          1.5         0.2  setosa
## 5          5.0         3.6          1.4         0.2  setosa
## 6          5.4         3.9          1.7         0.4  setosa
# select() - 필요한 칼럼만 인덱싱
head( select( iris, Sepal.Length) )
##   Sepal.Length
## 1          5.1
## 2          4.9
## 3          4.7
## 4          4.6
## 5          5.0
## 6          5.4
head( select( iris, Sepal.Length, Sepal.Width) )
##   Sepal.Length Sepal.Width
## 1          5.1         3.5
## 2          4.9         3.0
## 3          4.7         3.2
## 4          4.6         3.1
## 5          5.0         3.6
## 6          5.4         3.9
# select() - 특정조건의 칼럼명 뽑아내기
head( select( iris, starts_with( "Sepal")) ) # Sepal 로 시작하는 칼럼명만인덱싱
##   Sepal.Length Sepal.Width
## 1          5.1         3.5
## 2          4.9         3.0
## 3          4.7         3.2
## 4          4.6         3.1
## 5          5.0         3.6
## 6          5.4         3.9
head( select( iris, ends_with( "Width")) )
##   Sepal.Width Petal.Width
## 1         3.5         0.2
## 2         3.0         0.2
## 3         3.2         0.2
## 4         3.1         0.2
## 5         3.6         0.2
## 6         3.9         0.4
head( select( iris, contains( "etal")) ) # contain아님 contains
##   Petal.Length Petal.Width
## 1          1.4         0.2
## 2          1.4         0.2
## 3          1.3         0.2
## 4          1.5         0.2
## 5          1.4         0.2
## 6          1.7         0.4
# select() - 칼럼명 변경해서 인덱싱
head( select( iris, petal = Petal.Length) )
##   petal
## 1   1.4
## 2   1.4
## 3   1.3
## 4   1.5
## 5   1.4
## 6   1.7
head( select( iris, petal = starts_with("Petal")) ) # 인덱싱 칼럼 2개를 하나의 칼럼명에 넣으면 뒤에 숫자 1, 2 붙어서 칼럼명 생성 ex> petal1, petal2
##   petal1 petal2
## 1    1.4    0.2
## 2    1.4    0.2
## 3    1.3    0.2
## 4    1.5    0.2
## 5    1.4    0.2
## 6    1.7    0.4

filter() - 행인덱싱, 특정범주 인덱싱

head(airquality)
##   Ozone Solar.R Wind Temp Month Day
## 1    41     190  7.4   67     5   1
## 2    36     118  8.0   72     5   2
## 3    12     149 12.6   74     5   3
## 4    18     313 11.5   62     5   4
## 5    NA      NA 14.3   56     5   5
## 6    28      NA 14.9   66     5   6
head( filter( airquality, Temp > 70) )
##   Ozone Solar.R Wind Temp Month Day
## 1    36     118  8.0   72     5   2
## 2    12     149 12.6   74     5   3
## 3     7      NA  6.9   74     5  11
## 4    11     320 16.6   73     5  22
## 5    45     252 14.9   81     5  29
## 6   115     223  5.7   79     5  30
head( filter( airquality, Temp > 70 & Month > 5) )
##   Ozone Solar.R Wind Temp Month Day
## 1    NA     286  8.6   78     6   1
## 2    NA     287  9.7   74     6   2
## 3    NA     186  9.2   84     6   4
## 4    NA     220  8.6   85     6   5
## 5    NA     264 14.3   79     6   6
## 6    29     127  9.7   82     6   7

arrange()

# default는 오름차순
head( arrange( airquality, Ozone) )
##   Ozone Solar.R Wind Temp Month Day
## 1     1       8  9.7   59     5  21
## 2     4      25  9.7   61     5  23
## 3     6      78 18.4   57     5  18
## 4     7      NA  6.9   74     5  11
## 5     7      48 14.3   80     7  15
## 6     7      49 10.3   69     9  24
# desc(기준칼럼명) 내림차순
head( arrange( airquality, desc(Ozone)) )
##   Ozone Solar.R Wind Temp Month Day
## 1   168     238  3.4   81     8  25
## 2   135     269  4.1   84     7   1
## 3   122     255  4.0   89     8   7
## 4   118     225  2.3   94     8  29
## 5   115     223  5.7   79     5  30
## 6   110     207  8.0   90     8   9
# 기준을 2개로
head( arrange( airquality, desc(Ozone), Day) )
##   Ozone Solar.R Wind Temp Month Day
## 1   168     238  3.4   81     8  25
## 2   135     269  4.1   84     7   1
## 3   122     255  4.0   89     8   7
## 4   118     225  2.3   94     8  29
## 5   115     223  5.7   79     5  30
## 6   110     207  8.0   90     8   9

mutate() - 기존칼럼 + 수식으로 새로운 칼럼 생성(한 df return)

head( mutate( airquality, TempInc = (Temp - 32) * 5 / 9) )
##   Ozone Solar.R Wind Temp Month Day  TempInc
## 1    41     190  7.4   67     5   1 19.44444
## 2    36     118  8.0   72     5   2 22.22222
## 3    12     149 12.6   74     5   3 23.33333
## 4    18     313 11.5   62     5   4 16.66667
## 5    NA      NA 14.3   56     5   5 13.33333
## 6    28      NA 14.9   66     5   6 18.88889

summarise() with group_by()

# group_by()안하더라도, 집계한다.
summarise(airquality, mean(Temp, na.rm = T))
##   mean(Temp, na.rm = T)
## 1              77.88235
# group_by()와 함께
summarise(group_by(airquality, Month), mean(Temp, na.rm = T))
## # A tibble: 5 x 2
##   Month `mean(Temp, na.rm = T)`
##   <int>                   <dbl>
## 1     5                    65.5
## 2     6                    79.1
## 3     7                    83.9
## 4     8                    84.0
## 5     9                    76.9

연습해보기

  • COPD데이터에서 성별 남자만 추출 -> 흡연이력으로 group_by()후 BMI의 평균을 summarise
library(moonBook)
data(acs)
head(acs)
##   age    sex cardiogenicShock   entry              Dx   EF height weight
## 1  62   Male               No Femoral           STEMI 18.0    168     72
## 2  78 Female               No Femoral           STEMI 18.4    148     48
## 3  76 Female              Yes Femoral           STEMI 20.0     NA     NA
## 4  89 Female               No Femoral           STEMI 21.8    165     50
## 5  56   Male               No  Radial          NSTEMI 21.8    162     64
## 6  73 Female               No  Radial Unstable Angina 22.0    153     59
##        BMI obesity  TC LDLC HDLC  TG  DM HBP smoking
## 1 25.51020     Yes 215  154   35 155 Yes  No  Smoker
## 2 21.91381      No  NA   NA   NA 166  No Yes   Never
## 3       NA      No  NA   NA   NA  NA  No Yes   Never
## 4 18.36547      No 121   73   20  89  No  No   Never
## 5 24.38653      No 195  151   36  63 Yes Yes  Smoker
## 6 25.20398     Yes 184  112   38 137 Yes Yes   Never
acs %>%
  filter( sex == "Male") %>%
  group_by( smoking ) %>%
  summarise( mean(BMI, na.rm = T))
## # A tibble: 3 x 2
##   smoking   `mean(BMI, na.rm = T)`
##   <chr>                      <dbl>
## 1 Ex-smoker                   24.5
## 2 Never                       24.4
## 3 Smoker                      24.2

다운로드 : https://public.tableau.com/en-us/s/download
예시데이터 : https://github.com/is2js/jupyter_python/blob/master/Global%20Superstore.xls

단축키 : F7 - 전체화면(프리젠테이션 모드)
기본개념 : 범주칼럼은 열 / 숫자형칼럼은 행 / 위도,경도,지역명 등은 더블클릭으로 자동배정

범주별 숫자를 tableu로 정렬 / 색상(범례) / 필터링

  1. tablue public 버전을 설치한다.


  2. 예시데이터(엑셀)을 아래와 같이 열거나 or 파일을 tableu에 드래그 한 뒤,
    기존 엑셀파일 속 sheet가 왼쪽 메뉴에 나열되어있다.
    하나를 선택해 준다.



  1. 왼쪽메뉴에서 시트를 선택해주면, 아래 [시트1]이 생긴다.

    • 차원 : 범주형에 속하는 칼럼들
    • 측정값 : 숫자형에 속하는 칼럼들이 자동 배정된다.
  2. 차원에 있는 범주 중 category칼럼을 -> 열에 ,
    측정값에 있는 숫자형 sale칼럼을 -> 행으로 드래그해보자.
    열로 주면, 그래프의 x축의 범주로 / 행으로 주면 y축의 값으로 간다.

  3. 열(범주), 행(숫자)가 나열된 상태에서 우측상단의 [표현방식]을 바꾸면 그래프의 모양을 선택할 수 있다.

  4. 메뉴아이콘 중에 [마크 레이블 표시]를 누르면, 그래프상의 수치를 다 표현해준다.

  1. 메뉴아이콘 중 [오름/내림차순]을 선택하면, 원본 순서이외에 크기순으로 정렬해준다.

  2. 범주로서 열에 배정한 category칼럼을, 다시 왼쪽에서 가져와 [색상]에 드래그 해주면 색을 나타낸다.
    우측상단에 카테고리칼럼이 색상을 나타내고 있다는 범례도 표시된다.


  3. 범주로서 열에 배정한 category칼럼을, 다시 왼쪽에서 가져와 [필터]에 드래그 한뒤, category칼럼의 화살표를 내려 [필터표시]를 한 뒤, 우측에서 [클릭을 통해 필터링]해보자.









위도,경도,지역명을 가진 데이터에서 숫자형을 크기+색상gradient로 표시

  1. 위도/경도 데이터는 더블클릭만으로 자동 배정된다.

  1. country 칼럼도 더블클릭만으로 자동으로 마크의 [원]에 배정되는 것 같다.

  1. 현재까지는 위도+경도+county명으로 지도에 표시가 된 상태다.
    여기에 county별(원) 숫자형인 sale를 <원의 크기>로 나타내고 싶어서 [크기]에 드래그했다.
    sale(여러 행에 같은 것이 있을 수 있으니) 합계로 표시가 된다.
    county별(원별) sale(원의 크기)별로 마크의 크기가 결정되었다.

  2. Discount를 [색상]으로 추가하면, county(원)별 sale(원의크기)별 discount(원의 색상)으로 나타날 것이다.
    색상의 기본은 discount의 합계별로 파란색gradient로 나타난다.

큰 discount에는 빨간색 <-> 낮은 discount에는 파란색이 나타나도록 [색상] - [색상편집] - [빨간색 - 파란색 다중] - [반전](빨간색이 높은수치로 이동)을 선택해주자.






범주별 숫자데이터를, 연도별(페이지별)로 보기 + 시계열로 그리기

  1. [시트3]을 생성한 뒤, category -> 열 / sale -> 행에 배정하자

  2. ship date(년월일, 날짜)칼럼을 -> [페이지]에 배정하면 [플레이]버튼이 생긴다.
    날짜칼럼은 기본으로 [년]을 기준으로 그려진다.


  3. 플레이 버튼을 누르면, 년도별로 1년씩 증가하면서 해당 데이터가 변한다
    직접 년도를 선택할 수 있다.


  4. [표현방식]을 lineplot으로 바꾸면, 페이지로 년도별로 볼 순 없지만, 시게열 데이터가 된다.



대쉬보드 = 지금껏 만든 시트들 통합

  1. 하단이 메뉴에서 시트가 아닌 [새 대쉬보드]를 만든다.

  1. 기본적으로 크기가 작게 배정되어있으므로 [크기] - [범위] - [자동]을 선택해 최대화 시켜준다.

  2. 지금까지 만들어둔 시트1, 2, 3을 원하는 위치에 드래그해놓고 크기를 조절한다.





  1. 시트1에서 만든 category의 범주 [필터]는 시트1에만 작동한다.
    대쉬보드안의 모든 시트(1,2,3)에 작동할 수 있도록, [필터 우클릭] - [워크시트에 적용]
    • [이 데이터원본을 사용하는 모든 항목]을 적용해주면, 대쉬보드 안의 모든 시트에 다 적용된다.
    • 년도를 조절하는 [페이지]도 시트2에만 적용되어있다. 이것은 모든 페이지에 옵션이 안보이는 것 같다.



  1. [F7]로 프리젠테이션모드(전체화면)에서 사용해본다.


Untitled

카이제곱 (독립성) 검정 3+1가지 방법

  • 범주-범주 교차표를 이용한 연관성 확인
  • COPD데이터로 sex별 DM(당뇨)의 연관성 확인
# 데이터 준비
df1 <- read.csv("week4/4주차_코드/COPD_Lung_Cancer.csv", header=T)
head(df1)
##   PERSON_ID SEX DTH CTRB_PT_TYPE_CD   BMI BP_LWST BLDS TOT_CHOLE  HMG
## 1  10000065   1   0               8 25.73      90  112       177 12.2
## 2  10000183   1   0               9 22.77     100  105       177 15.7
## 3  10000269   1   0               8 24.54      70   76       156 13.1
## 4  10000471   1   0               1 21.63     100   70       228 14.0
## 5  10000788   2   0               7 25.65      54  113       199 12.1
## 6  10001096   2   0               3 19.84      60  137       182 12.0
##   FMLY_CANCER_PATIEN_YN SMK_STAT_TYPE DRNK_HABIT EXERCI_FREQ lungC copd
## 1                     1             1          1           1     0    0
## 2                     1             1          5           1     0    1
## 3                     1             1          1           1     1    1
## 4                    NA             1          1           1     0    1
## 5                     1             1          1           1     0    0
## 6                     1             1          1           1     0    0
##   Asthma DM Tub before_op_score after_op_score
## 1      0  1   0              44             52
## 2      1  1   0              44             53
## 3      0  1   1              37             60
## 4      0  1   0              50             44
## 5      1  1   0              31             43
## 6      0  0   0              34             49
# 방법1 : table( 칼럼인덱싱, 칼럼인덱싱 ) 을 chisq.test()에
chisq.test( table( df1$SEX, df1$DM )) # 0.8943
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  table(df1$SEX, df1$DM)
## X-squared = 0.017658, df = 1, p-value = 0.8943
# 방법2 : table없이 바로
chisq.test( df1$SEX, df1$DM ) # 0.8943
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  df1$SEX and df1$DM
## X-squared = 0.017658, df = 1, p-value = 0.8943
# 방법3 : xtabs( ~ 칼럼1+칼럼2, data= ) 를 summary()
x <- xtabs( ~ SEX + DM , data = df1 )
summary(x) #0.8412
## Call: xtabs(formula = ~SEX + DM, data = df1)
## Number of cases in table: 944 
## Number of factors: 2 
## Test for independence of all factors:
##  Chisq = 0.04013, df = 1, p-value = 0.8412
# 방법4 : mytable() 칼럼순서 상관x
library(moonBook) 
mytable( DM ~ SEX , data= df1) # 0.894
## 
##    Descriptive Statistics by 'DM'  
## ____________________________________ 
##             0           1        p  
##          (N=375)     (N=569)  
## ------------------------------------ 
##  SEX                           0.894
##    - 1 158 (42.1%) 236 (41.5%)      
##    - 2 217 (57.9%) 333 (58.5%)      
## ------------------------------------

비모수검정 2집단(mann-whitney, wilcoxon rank-sum) 3집단(kruskal wallis)

  • COPD 중 100개 행만 뽑아 정규분포를 안이루도록 만든다음
  • sex(2집단)별 BMI 평균차이 검정
  • SMK_STAT_TYPE(3집단)별 BMI 평균차이 검정
# 데이터 준비
df2 <- df1[1:100,]
head(df2)
##   PERSON_ID SEX DTH CTRB_PT_TYPE_CD   BMI BP_LWST BLDS TOT_CHOLE  HMG
## 1  10000065   1   0               8 25.73      90  112       177 12.2
## 2  10000183   1   0               9 22.77     100  105       177 15.7
## 3  10000269   1   0               8 24.54      70   76       156 13.1
## 4  10000471   1   0               1 21.63     100   70       228 14.0
## 5  10000788   2   0               7 25.65      54  113       199 12.1
## 6  10001096   2   0               3 19.84      60  137       182 12.0
##   FMLY_CANCER_PATIEN_YN SMK_STAT_TYPE DRNK_HABIT EXERCI_FREQ lungC copd
## 1                     1             1          1           1     0    0
## 2                     1             1          5           1     0    1
## 3                     1             1          1           1     1    1
## 4                    NA             1          1           1     0    1
## 5                     1             1          1           1     0    0
## 6                     1             1          1           1     0    0
##   Asthma DM Tub before_op_score after_op_score
## 1      0  1   0              44             52
## 2      1  1   0              44             53
## 3      0  1   1              37             60
## 4      0  1   0              50             44
## 5      1  1   0              31             43
## 6      0  0   0              34             49
# mann-whitney = wilcoxon rank sum test
# cf) wilcoxon signed rank -> signed = 절대값 = 차이를 이용 = 대응표본t-test의 비모수
wilcox.test( BMI ~ SEX, data = df2) # 0.6791
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  BMI by SEX
## W = 1189.5, p-value = 0.6791
## alternative hypothesis: true location shift is not equal to 0
wilcox.test( BMI ~ SEX, data = df2, exact = F) # 만약 동률이 생길 경우, exact = F 옵션줄 것
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  BMI by SEX
## W = 1189.5, p-value = 0.6791
## alternative hypothesis: true location shift is not equal to 0
# kruskal - wallis test 
kruskal.test( BMI ~ SMK_STAT_TYPE, data = df2) #0.1803
## 
##  Kruskal-Wallis rank sum test
## 
## data:  BMI by SMK_STAT_TYPE
## Kruskal-Wallis chi-squared = 3.4263, df = 2, p-value = 0.1803
Untitled

비모수 검정(Non-parametic)

  • 정규성을 만족하지 않거나 데이터가 작을 때 사용한다.

1. Mann-Whitney test (or Wilcoxon rank-sum test)

  • 독립표본 t-test (1범주-2집단에 대한 숫자형의 평균차이)에 대한 비모수적 검정
  • 2가지 방법이 있다.
# 데이터 준비
library(moonBook)
data(acs)

acs2 <- acs[1:100, ]
head(acs2)
##   age    sex cardiogenicShock   entry              Dx   EF height weight
## 1  62   Male               No Femoral           STEMI 18.0    168     72
## 2  78 Female               No Femoral           STEMI 18.4    148     48
## 3  76 Female              Yes Femoral           STEMI 20.0     NA     NA
## 4  89 Female               No Femoral           STEMI 21.8    165     50
## 5  56   Male               No  Radial          NSTEMI 21.8    162     64
## 6  73 Female               No  Radial Unstable Angina 22.0    153     59
##        BMI obesity  TC LDLC HDLC  TG  DM HBP smoking
## 1 25.51020     Yes 215  154   35 155 Yes  No  Smoker
## 2 21.91381      No  NA   NA   NA 166  No Yes   Never
## 3       NA      No  NA   NA   NA  NA  No Yes   Never
## 4 18.36547      No 121   73   20  89  No  No   Never
## 5 24.38653      No 195  151   36  63 Yes Yes  Smoker
## 6 25.20398     Yes 184  112   38 137 Yes Yes   Never
# 방법1 : mann-whitney test 1
#  - 물결로 칼럼명만 입력 & exact옵션주기 (a logical indicating whether an exact p-value should be computed.)
# **** exact 옵션 - 순서대로 나열 시, 같은 값이 존재하면 순서를 정하는데 문제가 생겨 디폴트인 exact test로 p-value를 구하지 못하고, 대신 정규 분포에 근사 시켜 p-value를 구한다. Exact=FALSE를 옵션으로 추가하여 정규분포에 근사 시키는 방법을 선택하면 디폴트인 exact test를 시도하지 않으므로 Warning을 없앨 수 있다.
wilcox.test( BMI ~ sex, data = acs2, exact = FALSE ) 
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  BMI by sex
## W = 946, p-value = 0.6327
## alternative hypothesis: true location shift is not equal to 0
# 방법2 : mann-whitney test 2
# - 직접 2개 집단(class, 수준)을 인덱싱하여 따로따로 데이터를 만든 뒤, 넣어주기
# - 이 경우는 exact = False옵션을 안주어도 된다.
male_bmi <- acs2[ acs2$sex == "Male", c("BMI")]
female_bmi <- acs2[ acs2$sex == "Female", c("BMI")]
wilcox.test(male_bmi, female_bmi)
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  male_bmi and female_bmi
## W = 1069, p-value = 0.6327
## alternative hypothesis: true location shift is not equal to 0

2. wilcoxon signed rank test

  • 대응표본 t-test(실험 전/후 데이터의 평균차이)에 대한 비모수적 검정
# 데이터 준비( 실험 전/후 -> 데이터 개수 같아야함 )
x1 <- c(51.4, 52.0, 45.5, 54.5, 52.3, 50.9, 52.7, 50.3, 53.8, 53.1)
x2 <- c(50.1, 51.5, 45.9, 53.1, 51.8, 50.3, 52.0, 49.9, 52.5, 53.0)

# wilcoxon signed rank test : paired t-test의 비모수 검정
wilcox.test(x1, x2, 
            alternative = c("greater"), paired = TRUE, conf.level = 0.95, # 95% 신뢰구간 -> 5%유의수준
            exact = F) # 오류나서 내가 추가해준 것
## 
##  Wilcoxon signed rank test with continuity correction
## 
## data:  x1 and x2
## V = 52.5, p-value = 0.006172
## alternative hypothesis: true location shift is greater than 0

3. kruskal wallis test 과 사후분ㅅ

  • 1way ANOVA의 비모수 검정
# 데이터 준비
head(acs2)
##   age    sex cardiogenicShock   entry              Dx   EF height weight
## 1  62   Male               No Femoral           STEMI 18.0    168     72
## 2  78 Female               No Femoral           STEMI 18.4    148     48
## 3  76 Female              Yes Femoral           STEMI 20.0     NA     NA
## 4  89 Female               No Femoral           STEMI 21.8    165     50
## 5  56   Male               No  Radial          NSTEMI 21.8    162     64
## 6  73 Female               No  Radial Unstable Angina 22.0    153     59
##        BMI obesity  TC LDLC HDLC  TG  DM HBP smoking
## 1 25.51020     Yes 215  154   35 155 Yes  No  Smoker
## 2 21.91381      No  NA   NA   NA 166  No Yes   Never
## 3       NA      No  NA   NA   NA  NA  No Yes   Never
## 4 18.36547      No 121   73   20  89  No  No   Never
## 5 24.38653      No 195  151   36  63 Yes Yes  Smoker
## 6 25.20398     Yes 184  112   38 137 Yes Yes   Never
# 1. kruskal wallis test
kruskal.test( weight ~ factor( smoking ), data=acs2)
## 
##  Kruskal-Wallis rank sum test
## 
## data:  weight by factor(smoking)
## Kruskal-Wallis chi-squared = 11.699, df = 2, p-value = 0.002882
# 2. 사후분석
# -  kruskal wallis test에서 3집단 중 적어도 하나의 집단이 차이가 있다고 할 경우
# 2_1. PMCMR 패키지 : kruskal willis test -> Tukey
# install.packages("PMCMR")
library(PMCMR)
## PMCMR is superseded by PMCMRplus and will be no longer maintained. You may wish to install PMCMRplus instead.
posthoc.kruskal.nemenyi.test(x=acs2$weight, g=as.factor(acs2$smoking), method="Tukey")
## Warning in posthoc.kruskal.nemenyi.test.default(x = acs2$weight, g =
## as.factor(acs2$smoking), : Ties are present, p-values are not corrected.
## 
##  Pairwise comparisons using Tukey and Kramer (Nemenyi) test  
##                    with Tukey-Dist approximation for independent samples 
## 
## data:  acs2$weight and as.factor(acs2$smoking) 
## 
##        Ex-smoker Never 
## Never  0.2784    -     
## Smoker 0.3899    0.0019
## 
## P value adjustment method: none
# install.packages("PMCMRplus") # 사용할려니 이것으로 깔아라고함( PMCMR은 더이상 지원 안할 예정)
# library(PMCMRplus)



# 2_2. n(on)par(ametic) comp(arison) 패키지 - mctp()함수
# 마찬가지로 Tukey로 하는 것 같다.
# - Estimation Method: Global Pseudo ranks 
# - Type of Contrast : Tukey 

# install.packages("nparcomp")
require(nparcomp)
## Loading required package: nparcomp
## Loading required package: multcomp
## Loading required package: mvtnorm
## Loading required package: survival
## Loading required package: TH.data
## Loading required package: MASS
## 
## Attaching package: 'TH.data'
## The following object is masked from 'package:MASS':
## 
##     geyser
result=mctp(weight~smoking, data=acs2)
## 
##  #----------------Nonparametric Multiple Comparisons for relative effects---------------# 
##  
##  - Alternative Hypothesis:  True differences of relative effects are less or equal than 0 
##  - Estimation Method:  Global Pseudo Ranks 
##  - Type of Contrast : Tukey 
##  - Confidence Level: 95 % 
##  - Method = Fisher with 33 DF 
##  
##  #--------------------------------------------------------------------------------------# 
## 
summary(result) # Analysis 탭으로 확인
## 
##  #----------------Nonparametric Multiple Comparisons for relative effects---------------# 
##  
##  - Alternative Hypothesis:  True differences of relative effects are less or equal than 0 
##  - Estimation Method: Global Pseudo ranks 
##  - Type of Contrast : Tukey 
##  - Confidence Level: 95 % 
##  - Method = Fisher with 33 DF 
##  
##  #--------------------------------------------------------------------------------------# 
##  
##  #----Data Info-------------------------------------------------------------------------# 
##              Sample Size    Effect     Lower     Upper
## Ex-smoker Ex-smoker   20 0.5054632 0.4243261 0.5863135
## Never         Never   39 0.3836336 0.3237192 0.4473026
## Smoker       Smoker   37 0.6109032 0.5455423 0.6725069
## 
##  #----Contrast--------------------------------------------------------------------------# 
##                    Ex-smoker Never Smoker
## Never - Ex-smoker         -1     1      0
## Smoker - Ex-smoker        -1     0      1
## Smoker - Never             0    -1      1
## 
##  #----Analysis--------------------------------------------------------------------------# 
##                    Estimator  Lower Upper Statistic     p.Value
## Never - Ex-smoker     -0.122 -0.309 0.074    -1.522 0.289976447
## Smoker - Ex-smoker     0.105 -0.093 0.296     1.296 0.403133946
## Smoker - Never         0.227  0.081 0.364     3.767 0.001777051
## 
##  #----Overall---------------------------------------------------------------------------# 
##   Quantile     p.Value
## 1 2.443167 0.001777051
## 
##  #--------------------------------------------------------------------------------------#

mytable()로 한다면? (사후분석은x)

  • method = 1 : 정규분포 가정 모수적 -> 2개 평균 +- 표준편차 제시하는 t-test
  • method = 2 : 비모수적인 통계방법 가정 -> 중앙값 [ 1, 3] 제시하는 비모수 2개 -> wilcox(대응표본), mann(독립표본)
  • method = 3 : 실제로 잔차의 정규성 검정(샤피로) 등등 해놓고 나서 알아서 판단하게 한다
head(acs2)
##   age    sex cardiogenicShock   entry              Dx   EF height weight
## 1  62   Male               No Femoral           STEMI 18.0    168     72
## 2  78 Female               No Femoral           STEMI 18.4    148     48
## 3  76 Female              Yes Femoral           STEMI 20.0     NA     NA
## 4  89 Female               No Femoral           STEMI 21.8    165     50
## 5  56   Male               No  Radial          NSTEMI 21.8    162     64
## 6  73 Female               No  Radial Unstable Angina 22.0    153     59
##        BMI obesity  TC LDLC HDLC  TG  DM HBP smoking
## 1 25.51020     Yes 215  154   35 155 Yes  No  Smoker
## 2 21.91381      No  NA   NA   NA 166  No Yes   Never
## 3       NA      No  NA   NA   NA  NA  No Yes   Never
## 4 18.36547      No 121   73   20  89  No  No   Never
## 5 24.38653      No 195  151   36  63 Yes Yes  Smoker
## 6 25.20398     Yes 184  112   38 137 Yes Yes   Never
library(moonBook)
# method = 2 : 강제로 비모수적 검정
# - 2집단  mann-whitney(독립표본), wilcox(대응표본), 3집단 kruskal - willis test
mytable(smoking ~ weight, data = acs2, method = 2) # 0.003 - 직접수행한것에 반올림만 하면 똑같다!
## 
##               Descriptive Statistics by 'smoking'              
## ________________________________________________________________ 
##            Ex-smoker          Never            Smoker        p  
##              (N=20)           (N=42)           (N=38)     
## ---------------------------------------------------------------- 
##  weight 60.0 [54.0;69.5] 59.0 [50.0;64.0] 65.0 [60.0;71.0] 0.003
## ----------------------------------------------------------------
# - 직접 kruksal
kruskal.test( weight ~ factor( smoking ), data=acs2) #  0.002882
## 
##  Kruskal-Wallis rank sum test
## 
## data:  weight by factor(smoking)
## Kruskal-Wallis chi-squared = 11.699, df = 2, p-value = 0.002882
# method = 3 : 잔차의 정규성 검정을 통해 판단후 모수적or비모수적 검정 ( 2집단 : t-test or mann/wilcoxon ) (3집단 : anova or kruskal )
mytable(smoking ~ weight, data = acs2, method = 3) # 직접 모수/비모수인지 잔차로 정규성을 검정한 뒤 -> 그에 따라 검정 -> 평균+-표준편차면 정규성O 모수적 /  중앙값[1Q, 3Q]면 정규성X 비모수적
## 
##        Descriptive Statistics by 'smoking'      
## _________________________________________________ 
##          Ex-smoker     Never      Smoker      p  
##           (N=20)      (N=42)      (N=38)   
## ------------------------------------------------- 
##  weight 62.6 ± 10.2 57.7 ±  9.0 66.0 ±  9.7 0.001
## -------------------------------------------------
Untitled

범주-범주의 교차표(contingency table) 만드는 2가지 방법

  • table()
  • descr패키지의 CrossTable()
# 데이터 준비
library(moonBook)
data(acs)


# 1. table()
table(acs$sex)
## 
## Female   Male 
##    287    570
table(acs$sex, acs$obesity)
##         
##           No Yes
##   Female 194  93
##   Male   373 197
# 2. descr패키지-CrossTable
# install.packages("descr")
library(descr)
a <- CrossTable(acs$sex, acs$obesity) # 빈도 뿐만 아니라 비율도 보여준다.
a
##    Cell Contents 
## |-------------------------|
## |                       N | 
## | Chi-square contribution | 
## |           N / Row Total | 
## |           N / Col Total | 
## |         N / Table Total | 
## |-------------------------|
## 
## ================================
##            acs$obesity
## acs$sex       No     Yes   Total
## --------------------------------
## Female       194      93     287
##            0.089   0.175        
##            0.676   0.324   0.335
##            0.342   0.321        
##            0.226   0.109        
## --------------------------------
## Male         373     197     570
##            0.045   0.088        
##            0.654   0.346   0.665
##            0.658   0.679        
##            0.435   0.230        
## --------------------------------
## Total        567     290     857
##            0.662   0.338        
## ================================

카이제곱 검정 - 관측값과 특정한 확률(기존의 기대값)와의 적합도 검정

  • H0 : 관측값과 기대값이 동일하다 ( chisq.test ( 범주별 빈도, p = 기대값 ))
  • 동전의 앞면/뒷면의 확률은 1/2인데, 관측결과 90/200, 110/200 이 나왔다면
  • 그게 적합하게 나온 건지 확인하는 것
# 1. 관측값 
obs <- c(19, 41, 40) # 19 / 전체 수 로 확률계산이 내부에서 될 것이고, 기대값과 비교할 듯..

# 2. 특정한 확률(기대값)
null.probs <- c(2/10, 3/10, 5/10)

# 3. 적합도 검정
chisq.test(obs, p = null.probs)
## 
##  Chi-squared test for given probabilities
## 
## data:  obs
## X-squared = 6.0833, df = 2, p-value = 0.04776
#### acs데이터로 해보기
# 1. table()로 범주의 빈도를 관측값으로 가져오기
smk_type <- table(acs$smoking) # 범주1개를 table()을 이용하여, 범주별 빈도를 변수에 담음
smk_type
## 
## Ex-smoker     Never    Smoker 
##       204       332       321
# 2. 기대값(특정확률)은 직접 변수에 입력해주기 (집단의 수만큼)
smk_type_prob <- c(0.3, 0.35, 0.35)

# 3. 관측값(범주별 빈도) + 기대값으로 카이제곱 적합도 검정
chisq.test( smk_type, p = smk_type_prob)
## 
##  Chi-squared test for given probabilities
## 
## data:  smk_type
## X-squared = 15.869, df = 2, p-value = 0.0003582

카이제곱 검정(Chisq test) - 2개 범주 독립성 검정의 3가지 방법대

  • chisq.test( 칼럼인덱싱, 칼럼인덱싱 )
  • H0 : 두 범주 X, Y는 독립이다.
  • 범주-범주의 교차표를 바탕으로 2개의 범주의 연관성(독립성)검정
  • X-squared : 카이제곱 검정통계량 -> 교차표상 각 cell의 ( (x-E)^2 / E)를 구해서 다 더한 것
# 1. 교차표(table) 생성 후 검정
table( acs$sex, acs$obesity )
##         
##           No Yes
##   Female 194  93
##   Male   373 197
chisq.test(table( acs$sex, acs$obesity ))
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  table(acs$sex, acs$obesity)
## X-squared = 0.30627, df = 1, p-value = 0.58
# 2. 바로 검정
chisq.test(acs$sex, acs$obesity)
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  acs$sex and acs$obesity
## X-squared = 0.30627, df = 1, p-value = 0.58
# 3. xtabs(교차표) & summary로 검정
xtabs( ~ sex + obesity, data = acs )
##         obesity
## sex       No Yes
##   Female 194  93
##   Male   373 197
summary( xtabs( ~ sex + obesity, data = acs ) )
## Call: xtabs(formula = ~sex + obesity, data = acs)
## Number of cases in table: 857 
## Number of factors: 2 
## Test for independence of all factors:
##  Chisq = 0.3968, df = 1, p-value = 0.5288

카이제곱 검정 다중비교 - 범주-범주이나, 독립변수의 집단이 3개 이상이어서 nC2로 -> 종속변수(2집단)에 독립인지 아닌지 본다.

  • fifer 패키지의 chisq.post.hoc( table() )
# 깃허브로 설치
# library(devtools)
# install_github("cran/fifer")

# 다중검정 in 3x3
library(fifer)
## Loading required package: MASS
table( acs$smoking, acs$Dx )
##            
##             NSTEMI STEMI Unstable Angina
##   Ex-smoker     42    66              96
##   Never         50    97             185
##   Smoker        61   141             119
chisq.post.hoc( table( acs$smoking, acs$Dx ) )
## Adjusted p-values used the fdr method.
##             comparison  raw.p  adj.p
## 1  Ex-smoker vs. Never 0.1112 0.1112
## 2 Ex-smoker vs. Smoker 0.0233 0.0350
## 3     Never vs. Smoker 0.0000 0.0000
chisq.post.hoc( table( acs$Dx, acs$smoking ) )
## Adjusted p-values used the fdr method.
##                   comparison  raw.p  adj.p
## 1           NSTEMI vs. STEMI 0.2964 0.2964
## 2 NSTEMI vs. Unstable Angina 0.0114 0.0172
## 3  STEMI vs. Unstable Angina 0.0000 0.0000

Fisher’s exact Test

  • 범주-범주의 교차표에서 각 cell의 Expected(기대값)이 5가 안되는 것이 25%이상(1/4이상)시 사용
# 교차표
table(mtcars$carb, mtcars$cyl)
##    
##     4 6 8
##   1 5 2 0
##   2 6 0 4
##   3 0 0 3
##   4 0 4 6
##   6 0 1 0
##   8 0 0 1
# 카이검정제곱 -> 카이제곱 approximation은 정확하지 않을수도 있습니다 -> Fisher로!
chisq.test(table(mtcars$carb, mtcars$cyl))
## Warning in chisq.test(table(mtcars$carb, mtcars$cyl)): Chi-squared
## approximation may be incorrect
## 
##  Pearson's Chi-squared test
## 
## data:  table(mtcars$carb, mtcars$cyl)
## X-squared = 24.389, df = 10, p-value = 0.006632
# fisher's exact test
fisher.test(mtcars$carb, mtcars$cyl)
## 
##  Fisher's Exact Test for Count Data
## 
## data:  mtcars$carb and mtcars$cyl
## p-value = 0.0003345
## alternative hypothesis: two.sided

Trend test

  • 순서를 가진 3집단이상의 범주(독립변수)에 대해 종속변수(2집단의 범주)가
  • H0 : 비율이 동일하다(일정하다)
  • H1 : 증가/감소 추세가 있다
# 1. trend-test 기본 
# 독립변수가 열 / 종속변수의 2개 집단이 행으로 나열된 교차표를 생각하자.
# 독립변수의 집단이 3개라고 가정 -> 3열
# 종속변수의 집단 2개 -> 2행 
prop.trend.test( c(15, 10, 30), # 종속변수 집단1 
                 c(35, 22, 44)) # 종속변수 집단2
## 
##  Chi-squared Test for Trend in Proportions
## 
## data:  c(15, 10, 30) out of c(35, 22, 44) ,
##  using scores: 1 2 3
## X-squared = 5.2588, df = 1, p-value = 0.02184
# 2. DescTools 패키지 - CochranArmitageTest()
# install.packages("DescTools")
library("DescTools")

# 데이터(matrix) 준비  2행 3열(독립변수 집단개수 = 열)
test <- matrix(c(2372,1859,2373,                 3065,2837,3065), byrow=TRUE, nrow=2)

# Desc( test )
# trend test 
CochranArmitageTest(test)
## 
##  Cochran-Armitage test for trend
## 
## data:  test
## Z = 0.011174, dim = 3, p-value = 0.9911
## alternative hypothesis: two.sided


5명의 발표자가 있었고, 나는 해당사항이 없었다.

나는 거의 BioInformatics에 필요한 기초 자료 학습에 집중하였고, 나만의 분석라이브러리도 만들었다. 또한, 민감한 의료데이터의 샘플을 받아 처리과정도 살펴보았다.

어느 인턴들보다 많은 학습을 했다고 자부하고, 교수님 & 연구실분들에게 크게 감사하고 있다.특히 논문에 사용되는 통계기법들을 학습하고 직접 r로 구현해본 것.
내 논문실험에 직접 대입해본 것
데이터를 보는 눈이 생긴 것
sql을 학습해본 것
등등 이다.

발표에 당선된 연구들은 거의다 해당 연구실의 대학원생들이 진행하던 연구를 그대로 받아온 것 같다. Figure나 연구결과들은 단기간(인턴기간)동안 나올 수 있는 결과물이 아니었다..
나는 직접 계획하고 실험하여 머신러닝 딥러닝 연구를 한의학에 적용하였지만, 의과학과 연구발표와는 거리가 있으므로 크게 발표나 수상에 염려를 두지 않았다.
하지만, 연구실마다 다양한 발표가 있었고, 그분들은 직접 실험하는데 보조를 했을 것이다.

나는 수상은 못했지만, 경험하지 못했고, 앞으로도 경험할 수 없는 것들을 얻고 간다!
부상으로 서울대 마크의 독서대를 받았다!

'한의대 생활 > 본4 - SNUBI 인턴생활' 카테고리의 다른 글

02. 19 분석 피드백  (0) 2019.02.20
02. 07 새로운 의료데이터 받기  (0) 2019.02.09
02. 01 첫 cdm 스터디  (0) 2019.02.01
01. 29 식권 대량구매  (0) 2019.01.30
01. 24 Pandas세미나 듣고옴  (0) 2019.01.24

비모수 통계 분석 ( Non parametic )

  1. 모집단의 분포가 정규분포가 아님 -> 표본은 도저히 정규분포를 가정 불가능
  2. 표본의 수가 n < 10으로 도저히 정규분포가 될 수 없는 경우
  3. 10 <= n <= 30 이나 but 정규성 검정(shapiro test)를 통과 못한 경우
  4. 순서를 가진 경우

모수 : 모평균, 모표준편차, 모분산 -> 현실적으로는 추출된 표본의 평균,표준편차,분산으로 -> 모집단의 모수들을 추정

모수 검정 : 정규성=모수적 특성=정규분포를 따르니, 서로 다른집단이라도 정규성을 이용하여 평균차이로 집단의 차이를 밝힌다. 표본이 정규성을 만족한다면, 모집단도 정규분포 -> 표본 집단들의 평균차이 = 모집단들의 평균차이 = 모집단들이 차이난다

비모수 검정 : 실제 표본들로 만들어진 평균이 아니라 표본들을 나열한 뒤, 부호,절대값(signed)이나 순위(rank)를 이용해서 검정한다. 대신 검졍력은 모수적 방법들에 비해 떨어진다.

제시 통계량 : 모수 검정시에는 평균 +- 표준편차을 제시, 비모수 검정시에는 중앙값 + [1QR, 3QR]을 제시

cf)

비모수 통계 분석의 종류

1~2집단에 대한 숫자형의 평균차이

  1. 대응표본 t-test -> Wilcoxon signed rank test
  2. 독립표본 t-test -> Mann-Whitney test or Wilcoxon rank-sum test

3집단의 숫자형의 평균차이(요인1개)

  1. one way ANOVA -> Kruskal-Wallis test

1. Wilcoxon signed rank test ( 대응표본 t-test의 비모수 검정 )

H0 : before와 after의 크기가 같다.

  1. 짝지어전 실험 전/후 자료를 빼서 diff를 구한다.
  2. diff에 절대값(signed)를 취한 뒤, 크기를 내림차순(rank)나열한다
  3. 동률 처리로서 절대값 diff값이 같은 애들은 순서를 더해서 평균을 준다.
    • 예를들어, 1 -1 -> rank 1, 2 -> rank 1.5, 1.5
  4. 절대값을 취하기전의 부호끼리 그 순위들을 다 더한다(diff순위의 부호별 합).
    • diff가 양이었던 애들의 순위들 합 vs 음이었던 애들의 순위들 합
  5. 양의 순위합과 음의 순위합을 비교하여 크기 차이가 있는지 검정(자세한 내용 확인필요)

2. Mann-Whitney test ( 독립표본 t-test의 비모수 검정 )

Wilcoxon rank-sum test라고도 하는 것 같다.

H0 : 두군의 크기가 같다

  1. 독립된 2집단의 자료들을 집단 관계없이 모두 섞는다(n수 늘리는 작업)
  2. 내림차순으로 정렬하여 순위(rank)
  3. 동률 처리
  4. 섞기전의 집단별로 순위들의 합(rank-sum)을 구한다
  5. 집단별로 순위합을 비교하여 크기 차이가 있는지 검정

3. Kruskal-wallis test( one way ANOVA의 비모수 검정 )

H0 : 세군의 평균이 모두 같다

  1. 3집단의 자료들을 집단과 관계없이 모두 섞는다.

  2. 내림차순으로 정렬한다.

  3. 동률처리는 더해서 평균

  4. 다시 집단별 순위합의 평균

  5. 비모수라도 사후검정이 가능하다.

+ Recent posts