생물학, 통계학, 데이터 과학

이제서야 올리는 오타와 튤립축제 포스팅

오타와 튤립축제는 5월 초에 열립니다. 매해 날짜가 살짝 달라지니 여기서 참고하시면 되겠어요. http://tulipfestival.ca/ 에버랜드 튤립축제에 비해 훨씬 떨어진다라는 얘길 듣고 갔는데 네, 그 말이 맞았어요. 한국과 비교하면 안됩니다. 페스티벌은 2주간에 걸쳐 열리는데 일찍 가면 꽃이 다 안 필 수 있어요.

오타와 관광은 1박 2일 정도면 어느정도 다 돌 수 있을만큼 작은 도시고요. 호텔을 다운타운 근처에 잡으면 걸어다니면서 구경하는것도 가능합니다. 다만 튤립축제는 다운타운에서 열리는게 아니라 근처 공원에서 열리는거라 차 타고 가야해요. 셔틀버스도 운행한다고 들었지만 전 그냥 차 끌고 갔습니다.

주말이어서 과연 주차할 수 있을지 걱정했는데, 아니나 다를까 공원 입구부터 차가 많더라고요. 공원 입구 맞은편은 주택가인데 주택가 돌아다니다보면 street parking 혹은 parking lot이 있는데 제가 주차했을땐 3불 이었어요. 이것도 어디에 주차하느냐에 따라 가격이 달라지니 잔돈은 여유롭게 늘 준비하셔야 합니다. 

주말엔 불꽃놀이도 한다고 들었지만 전 그냥 호텔에 들어가서 쉬었습니다. 가볍게 한바퀴 돌다 오자란 마음으로 간거여서 여행 일정에 대한 기록보단 사진만 올려봅니다. 단, 공원 근처에 먹을게 풍성하지 않았어서, 돗자리랑 간식거리 들고왔으면 더 좋았을텐데 아쉬움은 있었어요.  

http://statnmath.blogspot.ca/2015/08/case-study-2x2-contingency-table-in-sas.html (개념설명) 참고하세요 :)  

이번 포스팅은 2x2 contingency table(분할표)에 대한 내용입니다. 뭔가 어려워보이지만 간단하게 설명해보자면 우리가 기대하는것과 실제로 얻은 값을 비교해 정말 차이가 있는지 없는지를 보는건데요. 예를들어 1부터 6까지 나올 확률이 동일한 fair 주사위를 60번 던져봅시다. 우리가 예상하는건 각각 1부터 6번까지 나올 확률이 동일하니까 숫자 1은 10번 나올테고 숫자 2는 10번 나올테고 등등...숫자 6은 10번 나오겠지~라고 예상할 수 있겠죠. 하지만 실제로 해보니까 1번은 8번 나오고 2번은 14번 나오고... 이렇게 들쭉날쭉합니다. 그래서 우리가 예상한것(expected)과 관찰한것(observed)이 얼마다 다른지 한번 비교해본게 바로 Chi-Squared test라고 보면 되겠어요. 

예제로 바로 들어가보겠습니다. Framingham Heart Study입니다.   

Reference: https://www.framinghamheartstudy.org/, https://en.wikipedia.org/wiki/Framingham_Heart_Study  

저희가 알아볼 데이타는 총 샘플이 1329 남자이고요. 1948년 콜레스테롤을 재서 High 와 Low  두가지로 나눴습니다. 그리고 10년 후, cardiovascular disease (CVD, 심혈관계질환)이 발생했는지 아닌지 나눴습니다. 이때 중요한점은 1329명의 샘플 숫자는 고정되어있고요. 그리고 콜레스테롤과 심혈관계질환은 카테고리 데이타 입니다. 이게 왜 중요하냐면, 1392명이 고정되지 않으면 Poisson distribution을 띄기 때문이예요. 그리고 count data가 아니라 proportions(비율이나 분수)이면 계산 자체가 안됩니다.

그래서 우리가 알고싶은건 과연 콜레스테롤 레벨이 심혈관계 질환과 연관되어있는지 아닌지 살펴보는 겁니다. 즉 두 개의 변수가 독립적인지 아닌지 살펴보려고 해요.

데이터 불러오는 SAS 코드는 다음과 같습니다. 이렇게 입력하면 아래 이미지처럼 데이타가 입력되어요.

data fram;
  input chol $ cvd $ count;
  datalines;
low present 51
low absent 992
high present 41
high absent 245; 

C와 D의 변수가 서로 독립적인지 아닌지 살펴보려면 다음을 살펴보면 됩니다. P(CD)=P(C) x P(D) 이때 P(CD)는 joint distribution(결합분포)이고 P(C)와 P(D)는 marginal distribution(주변분포)입니다. 두 값을 구해서 서로 같은지 아닌지 살펴보면 되겠지요.

The Joint distribution:

The probability that an observation falls into row i, column j, for i & j=1, 2 $= P(C=i, D=j)= \pi_{ij}$

The Marginal distribution:

The probability an observation falls into row I $= P(C=i)= \pi_{i \cdot}$

The probability an observation falls into column j $= P(D=i)= \pi_{ \cdot j}$

$H_0$ : $\pi_{ij}=\pi_{i\cdot}\pi_{\cdot j}$ There is no relationship between Cholesterol and CVD.

$H_1$ : $\pi_{ij} \neq \pi_{i\cdot}\pi_{\cdot j}$  

SAS에서 proc freq로 Chi-squared value를 구할 수 있습니다. 이걸로 서로 두 변수가 독립적인지 아닌지 알아볼 수 있어요. 

proc freq;

 weight count;

 table chol*cvd / chisq;

run;

SAS 결론은 아래 링크 참고하세요.

http://statnmath.blogspot.ca/2015/08/intro-to-log-linear-model-ixj.html (예제) 

이번 포스팅은 Estimator(추정량)과 Estimates에 대한 용어설명 입니다. 

먼저, Population (모집단) & Parameter (모수) // Sample (표본) & Statistic (통계치) // Model (모형,모델) 용어는 http://statnmath.tistory.com/95 여기 참고하세요.

Population(모집단)에 대해 알고싶다면, 이를 잘 설명할 수 있는 모델을 찾아야겠죠. 이 모델은 parameter(모수)로 이뤄져 있고요. 모집단으로 모델을 찾기에 불가능하니까, sample(표본)을 구해 statistic(통계치)를 구합니다. 그렇게 해서 모집단을 설명할 수 있는 모델을 만드는데요. 이 모델에서 모수가 어떤 값이라고 추정하는 값을 바로 estimate, 추정값이라고 합니다. 

이때 estimator는 추정값을 추정하는 방법(method)을 말합니다. 간단하게 생각해서 estimates는 수치(값)를 말하는거고, estimator는 이 수치를 구하게 된 방법이라고 생각하면 되겠습니다.

그래서 estimator는 보통 $\hat{\mu}$ 이렇게 표시하는데요. 이걸 영어로 hat이라고 부릅니다. 단순히 $\mu$ 만 있으면 이건 parameter를 말하는거고요. hat이 씌워져있으면 ($\hat{\mu}$) 추정량을 말합니다.  

Population, Sample, Model에 대한 용어 정리입니다.

Population은 모집단이고, Sample은 표본 그리고 Model은 모델(모형)이라고 합니다. 예를들어 대한민국 대학 졸업생의 평균나이를 구하려고 합니다. 그렇다면 모든 대학교의 졸업생 전체를 다 알아봐야겠죠. 이게 바로 Population, 즉 모집단을 말합니다. 하지만 이렇게 전체 알아보는게 불가능하기때문에 전국 10개 대학의 졸업생을 따로 조사해봅니다. 이게 바로 Sample, 표본이예요.

자, 여기서 생각해야할게 Population에서 구하고자 하는 값(평균나이)을 parameter, 모수라 합니다. 이 모수는 알 수가 없어요. 알 수 없다는게 이 값을 구한다는게 불가능하죠. 만약 가능했으면 sample(표본)을 뽑지도 않았겠죠. 그럼 sample(표본)을 전국 10개 대학교를 뽑을 수도 있고, 예산이 부족해서 3개 대학을 뽑을 수도 있잖아요. 이처럼 10개 대학을 뽑을때와, 3개 대학을 뽑을때의 수치(평균나이)가 달라질 수 있습니다. 이 수치를 바로 statistic(통계치)라고 합니다.  

Model(모델, 모형)은 population이 어떤지 리서처들이 생각하는 바를 식으로 나타낸다고 생각하면 될것 같네요. 물론 model이 population과 정확히 일치하진 않지만 그래도 표본을 통해 population과 비슷하도록 모형화합니다.

다시 정리해보자면, $\mu$는 모집단의 평균나이를 말하는 거고요. 알 수 없는 수치입니다. 하지만 sample, 표본을 통해 얻은 평균나이는 알 수 있죠. 이걸 $\bar{X}$ 로 표현합니다. 모집단(population)의 모수(parameter)를 모델로 만든다면 이때 $\mu$와 다르게 표현하기 위해 $\hat{\mu}$ 이라고 적습니다.

다음 포스팅에서 Estimator와 Estimates등 이어 설명하도록 할게요.

Reference : Ramsey, F.L. and Schafer, D.W. (2013). The Statistical Sleuth: A Course in Methods of Data Analysis (3rd ed), Cengage Learning. (Data from Vaisanen and Jarvinen, “Dynamics of Protected Bird Communities in a Finnish Archipelago,” Journal of Animal Ecology 46 (1977): 891-908.)

[1] Data and Model

아래 이미지에서 원래 데이타는 island, area (섬 사이즈), nspecies (1949년 각 섬마다 발견된 새 종의 수, X), nextinct (1959년 멸종된 새 종, m)까지만 있고요.  그 이후는 제가 binomial logistic regression돌린 값입니다. Response variable, Y는 nextinct 입니다.

$\pi_i$ is probability of extinction of each island, assuming species survival is independent.

Then, $Y_i\sim Binomial (m_i,\pi_i)$

Observed response proportion is $\hat{\pi_i}=\frac{y_i}{m_i}$ = EXTINCT / ATRISK

Observed logit : $\log(\frac{\pi_{S,i}}{1-\pi_{S,i}})$, where S means "Saturated" 

Model : $\log(\frac{\pi_{S,i}}{1-\pi_{S,i}})= \beta_0+\beta_1 \cdot AREA_i$

[2] Initial Assessment

먼저 Binomial Logistic Regression에 들어가기 전, logit과 X변수간의 linear relationship이 있는지 살펴봐야합니다. 입니다!! Scatter Plot으로 logit과 area를 살펴보니 선형관계가 아니군요. 

Area에 log를 취하고 다시 scatter plot을 해보니 선형관계가 보이네요. 그래서 최종 모델은 아래와 같습니다.

Our model :  $\log(\frac{\pi_{S,i}}{1-\pi_{S,i}})= \beta_0+\beta_1 \cdot \log(AREA_i)$

[3] SAS Result

Number of observations are 18, and sum of frequencies is 632 (=$\sum_{i=1}^{18}m_i$ )

Fitted model: $logit(\hat{\pi})= - 1.196-0.297 \log(AREA)$

AREA 변수가 과연 significant 한지 알아보려면 WALD Test를 이용하면 됩니다!!  

Testing Global Null Hypothesis : BETA=0 (Wald Test) (Log Area is significant?) 

- Likelihood Ratio statistic is 33.2765 which is computed by  578.013 - 544.736, and the degree of freedom is 1 as we are testing only $\beta_1$.

- 결론은 P-value가 0.0001이라서 log(AREA)의 계수가 0이 아니라는거죠. 그래서 log(AREA) 변수가 logit 값을 위해 필요하다는걸 알 수 있네요.

$\beta_1$의 해석하기, 멸종확률 구하기, Deviance Goodness of Fit Test 에 대한 SAS 결과값 해석방법은 위 링크 참고하시면 정리되어 있습니다. 궁금하신 사항은 댓글 주세요 :)

참고로 Deviance GOF - reduced 모델과 full 모델 중에 어느 모델이 더 적합한지 알아보는 테스트입니다. 결과값은 링크된 곳에 정리되어 있습니다.

$H_0$: Our fitted model fits data as well as saturated model. 

$H_1$: Saturated model $logit(\pi)= \beta_0+ \beta_1 I_1+...+ \beta_{n-1}I_{n-1}$ is better. 

Test Statistics: Deviance = $-2 \log \frac{L_R}{L_F}=12.06$ with n-(p+1) = 18 -1 + 1 =16 d/f, p: # testing $\beta$