크루스칼-왈리스 검정은 모수적 방법에서의 one-way anova와 같은 목적으로 쓰입니다.

대신 그룹별로 평균이 아닌 중위수가 같은지 아닌지를 검정합니다.

A (98,95,76,95,83,99,82,75,88)
B (95,71,80,81,77,70,80,72,81)
C (77,79,91,84,80,93,87,81,83)

예를 들어, 아래와 같은 세 그룹이 있는 이 세 그룹간의 중위수의 차이가 있는가를 알아보기 위하여 쓰입니다.

즉, H0 : A의 중위수 = B의 중위수 = C의 중위수 vs H1 : 최소 한 그룹의 중위수가 다르다.

이를 검정하기 위한 검정통계량은 아래와 같습니다.

이 때 N은 전체 표본의 수이고, ni는 각각의 그룹의 표본 수 입니다. ri는 표본들을 모두 합친 후 순위를 냈을 때의 그룹별 데이터의 순위값입니다. 만약 중위수가 정말 같다면(H0=True) ri의 평균은 비슷할 것입니다. 다르다면 이 값이 커지게 되고, 유의하게 컸을 때, H0를 기각하게 됩니다.

H의 정확한 분포를 구하는 것은 매우 복잡하기 때문에 그룹당 표본수가 충분하면 H가 자유도 (k-1)인 카이제곱 분포를 따른다는 것을 이용하여 검정할 수 있습니다. (이 때 k=그룹의 수)

R 을 통하여 크루스칼 왈리스 검정을 하는 법은 아래와 같습니다.

value1 <- c(98,95,76,95,83,99,82,75,88)
value2 <- c(95,71,80,81,77,70,80,72,81)
value3 <- c(77,79,91,84,80,93,87,81,83)
category <- c(rep("A",9),rep("B",9),rep("C",9))
data <- data.frame(value=c(value1,value2,value3), category)

kruskal.test(value~category, data=data)

결과

Kruskal-Wallis rank sum test


data:  value by category


Kruskal-Wallis chi-squared = 5.6972, df = 2, p-value = 0.05793

모수검정과 비모수검정 비교표

참고

http://rfriend.tistory.com/130

윌콕슨 순위합 검정은 두 표본의 중위수를 비교하는데 쓰인다. 모수적 방법에서 이표본 t-test와 같다.

A, B 회사로부터 각각 생산된 계피가루 한 통의 무게를 X,Y라고 하자 n1=8, n2=8로 표본을 추출하였다.

X = {117.1, 121.3, 127.8, 121.9, 117.4, 124.5, 119.5, 115.1}

Y = {123.5, 125.3, 126.5, 127.9, 122.1, 125.6, 129.8, 117.2}

H0 : X의 중위수 = Y의 중위수 vs H1 : X의 중위수 != Y의 중위수

이 때 X, Y가 결합된 표본을 만들고 이 표본을 오름차순으로 정렬해보자.

XY = {115.1, 117.1, 117.2, 117.4, 119.5, 121.3, 121.9, 122.1, 123.5, 124.5, 125.3, 125.6, 126.5, 127.8, 127.9, 129.8}

그리고 Y 값들에 대하여 밑줄을 쳐놓고, 이 Y 값들의 순위를 모두 합쳐보자.

W = 3+8+9+11+12+13+15+16 = 87

이 때 W를 윌콕슨 통계량이라고 한다. 만약 이 W 값이 너무 작거나 너무 크면 H0를 기각할 유의한 증거가 있는 것이다. 그러면 우선 H0하에서 W의 분포를 알아야하는데, W의 정확한 분포를 알기는 매우 힘드므로 중심극한 정리를 이용한다.

우리가 알 수 있는 것은 E[W] = n2(n1+n2+1)/2, V[W] = n1n2(n1+n2+1)/12 이다. 그러므로 중심극한정리를 이용하여

W ~ N(n2(n1+n2+1)/2, n1n2(n1+n2+1)/12)

n1=8, n2=8을 이용하여 평균과 분산을 각각 구하면, 68, 90.66 이다.

위에서 구한 W 값을 표준화 시키면,

z = (87-68)/sqrt(90.66) = 1.99

양측검정에서 기각역은 z > 1.96 or z < -1.96이므로 H0를 기각한다.

Wilcoxon Signed Rank Test는 일표본에 대하여 중위수를 검정하는데에 쓰인다. 또는 paired 된 데이터의 중위수의 차이가 있는지를 검정하는데 쓰인다. 모수적 방법에서 일표본 t-test와 같은 목적으로 쓰인다.

표본이 176.9, 158.3, 152.1, 158.8, 172.4, 169.8, 159.7, 162.7, 156.6, 174.5, 184,4, 165.2, 147,8, 177.8, 160.1, 160.5 라고 하자.

H0 : 중위수 = 160 vs H1 : 중위수 != 160 를 검정하도록 해보자.

우선 | Xi - 160 | 을 오름차순으로 차례대로 줄 세운다.

0.1, 0.3, 0.5, 1.2, 1.7, 2.7, 3.4, 5.2, 7.9, 9.8, 12.2, 12.4, 14.5, 16.9, 17.8, 24.4

그러고 나서, Xi-160 < 0 인 데이터에 대해서는 표시를 해놓는다. 그리고 부호순위를 모두 더한다. 밑줄친 부분에 대한 순위는 음수로 계산한다.

W = 1-2+3-4-5+6-7+8-9+10-11+12+13+14+15+16 = 60

만약 W의 값이 너무 크거나 너무 작으면 중위수가 160이라는 귀무가설을 기각하게 된다. 그러면 그 값이 유의하게 큰지 작은지는 어떻게 알 수 있을까? 이를 위해선 우선 H0하에서 W의 분포를 알아야한다. 하지만 실제로 W의 분포를 알기는 매우 어렵고 복잡하다. 이 때 우리가 알 수 있는 것은 W통계량의 평균은 0, 분산은 n(n+1)(2n+1)/6 이라는 것이다. (H0하에서 중위수의 성질을 이용해 평균과 분산을 구할 수 있다.) 이 값들로 중심극한정리를 사용하여 W의 분포를 정규분포로 근사시킨다. 이렇게 할 수 있는 이유는 W가 어떠한 분포를 갖는 값들의 합이기 때문이다. 예를 들어, 이 경우에 H0하에서 위 식의 첫번째 term은 1/2 확률로 1, 1/2 확률로 -1/2 이다. 두 번째 term은 1/2확률로 2, 1/2 확률로 -2이다. 이 값들을 모두 합친 W 통계량은 종모양의 분포가 나올 것이다. 그러므로 중심극한 정리를 적용시킬 수 있다.

n=16이므로, W를 표준화 시키면 Z = 60-0/sqrt((16)(17)(33/6)) = 1.525 이다.

양쪽 검정이므로 유의수준 0.05에서 기각역은 Z>1.96 or Z<-1.96이다. 계산한 z값이 기각역에 포함되지 않으므로 중위수가 160이 아니라는 유의한 증거가 없으므로 H0를 기각하지 못한다. 또한 p-value = 2*P[z > 1.525] 이다.

비모수 통계학이란 모집단의 분포를 가정하지 않고, 오직 데이터를 통해서만 데이터를 검정하는 방법을 의미한다. 비모수 통계학은 영어로 non-parametric statistics라 하고, distribution-free 방법이라고도 말한다.

모수적 방법에서 t-test의 예를 들어보면, 기본 가정중의 하나가 바로 모집단의 "정규성" 이다. 모집단이 정규분포가 아니라면 t-test의 신뢰도는 낮아지게 된다. 그러므로 보통 비모수 방법을 사용하는 이유가 바로 모수적 방법에서의 가정이 만족하지 않을 때이다. t-test의 경우에는 모집단에 대하여 정규성 검정을 수행한 후 정규분포가 아니라는 근거가 있을 때는 비모수적 방법을 사용하면 된다.

비모수 통계학의 장점은 아래와 같은 것들이 있다.

1. 모집단의 분포가 어떻든지 사용할 수 있다.
2. 모수적 방법에 비해 계산이 간단하며, 적용하기 쉽다.
3. 관측값에 신뢰성이 적어 단순히 순위를 통해 검정하고 싶을 때 이용하면 좋다.

반면 비모수적 방법은 아래와 같은 단점이 있다.

1. 모수적 방법으로 검정할 수 있는 데이터에 비모수적 방법을 이용하면, 효율성이 떨어진다.
2. 표본의 크기가 커질 수록 비모수적 방법의 계산량은 빠르게 늘어난다.

<출처 - 위키피디아>

위 그림은 분자생물학에서 non-coding RNA를 포함한 central dogma를 나타낸 것이다. 

DNA에는 mRNA를 만드는 gene만 있는 것이 아니라 다른 non-coding RNA를 만드는 gene도 있다. tRNA의 예를 들면 tDNA라고 표시된 tRNA를 만드는 DNA 상의 gene이 RNase P라는 효소에 의해 분해되어서 tRNA를 만든다.

단순한 central dogma (DNA=>mRNA=>Protein)만 알고 있었는데, 이렇게 non-coding RNA까지 포함하여 한 번에 나타내주는 그림이 있어서 공부하는데 도움이 되었다. 

tRNA란 무엇인가?

tRNA란 mRNA의 코돈에 맞는 아미노산을 리보솜으로 운반하는 역할을 합니다. tRNA는 리보솜과 함께 단백질을 생산합니다. 

tRNA에는 안티코돈과 아미노산이 결합하는 자리가 있습니다.

tRNA의 구조

• tRNA의 구조는 3'-ACC-5' 로 되어있으며 3' 쪽에 아미노산이 공유결합으로 붙어서 리보솜으로 전달됩니다.

• ACC 서열은 CCA서열은(5'->3'쪽으로 읽어야하므로) 모든 tRNA에 공통적으로 존재합니다.

• 안티코돈은 코돈에 상보적인 염기서열입니다.

• 안티코돈은 32종입니다. 안티코돈은 45종입니다. 코돈은 64종인데 안티코돈은 45종인이유는 똑같은 기능을 하는 안티코돈이 있기 때문입니다. 안티코돈중 세번째 염기가 U인 안티코돈은 A,G 모두와 쌍을 이룹니다. (wobble pair) (출처 - 캠벨 생명과학)

• 위의 두 그림을 보면 보다 자세히 tRNA의 구조를 볼 수 있습니다.

• 몇몇 시퀀스가 상보적인 결합(수소결합)을 이루고 있기 때문에 위처럼 입체적인 구조를 갖고 있습니다.

• 루프가 존재하며 이 루프는 상보적인 결합을 이루고 있지 않습니다.

• 3'과 5'를 손으로 잡고 쭉편다고 생각하면 이것이 tRNA의 시퀀스가 됩니다.

tRNA는 리보솜과 함께 단백질을 생산한다.

• 가장 위 그림처럼 개시코돈을 보면 이와 상보적인 안티코돈을 갖는 tRNA가 와서 붙습니다.

• 위에서 두 번째 그림에서 보면 두 번째 tRNA가 그 다음 코돈에 가서 붙습니다. 그러면 첫번째 tRNA는 자신의 아미노산을 두번째 tRNA로 옮깁니다. 

• 그 다음 리보솜이 한 칸 이동합니다. 위에서 두 번째에서 세 번째 그림

• 이와 같은 과정을 계속 반복하다가 종결코돈을 보면 Release factor가 붙게되고 마지막으로 단백질 리보솜 tRNA가 모두 분리됩니다.

참고

http://snustudy.com/220404538909

http://blog.naver.com/dyner/100174698698

R - 질적변수의 교호작용을 고려한 회귀모델

/* 2017.5.21  */

데이터의 형식은 아래와 같습니다.

  case gram weeks smoke

1    1 3147    40     0

2    2 2977    41     0

3    3 3119    38     0

4    4 3487    38     0

5    5 4111    39     0

6    6 3572    41     0

gram : 태아의 몸무게 

weeks : 임신주차

smoke : 임신기간 중 흡연 여부 

임신주차와 흡연여부 별로 태아의 몸무게가 어떻게 변하는지 설명하는 회귀식을 찾는 것입니다. 임신주차가 길 수록 태아의 몸무게는 무겁고, 흡연을 안할 수록 태아의 몸무게가 무겁다는 것을 예상해볼 수 있습니다.

이 때, smoke는 질적변수이며, 회귀모델을 구축하는 것에는 두 가지 종류가 있을 수 있습니다.

첫 번째는 smoke라는 질적변수에의해 회귀모델의 절편만 차이가 나는 경우

두 번째는 smoke라는 질적변수에 의해 절편과 기울기가 둘 다 차이가 나는 경우 (교호작용을 고려한 경우)

즉, 흡연여부에 의해 임신주차에 따른 태아의 몸무게의 변화가 영향을 주는 것까지 고려한 경우입니다.  이처럼 한 변수가 다른변수의 종속변수에 대한 효과에 영향을 끼치는 것을 다른말로 교호작용이라고 합니다. 예를 들어 흡연을 하는 경우, 임신주차에 따른 태아 몸무게의 증가 폭이 완만할 수 있겠죠.

두 번째 경우에 대해서만 R을 이용하여 모델을 구축하여보겠습니다.

data <- read.csv("C:/태아몸무게.csv")
head(data)

# 교호작용도 고려한 질적변수 처리
model2 <- lm(gram~weeks*smoke, data=data)
coef <- coefficients(model2)

inter <- coef[1]
weeks <- coef[2]
smoke <- coef[3] # 절편의 차이
weeks_smoke <- coef[4] # 기울기의 차

y1 <- inter+weeks*0:50 # non-smoker
y2 <- inter+(weeks+weeks_smoke)*0:50+smoke # smoker

plot(data$weeks, data$gram)
lines(y1, col="red") 
lines(y2, col="blue")

• model2 <- lm(gram~weeks*smoke, data=data) 라고하면 weeks와 smoke의 모든 경우에 수를 고려하여 회귀식을 만들겠다는 것입니다.
• 그러면 weeks의 회귀계수(beta1)와 smoke의 회귀계수(beta2), weeks*smoke(beta3)의 회귀계수가 생기게 됩니다.
• beta2는 smoke효과로 인한 절편의 차이이며, beta3는 smoke효과로 인한 기울기의 차이입니다.
• 일반적인 모델에서는 beta3가 없고 절편의 차이만 있는데 교호작용까지 고려하면 기울기의 차이도 고려하는 것입니다.

붉은색은 smoke=0 즉 비흡연자이고, 파란색은 smoke=1 즉 흡연자입니다. 흡연자는 임신주차가 증가하더라도 태아의 몸무게가 완만하게 증가하는 경향을 볼 수 있습니다. 절편값은 흡연자가 더 높은 것을 볼 수 있는데, 이것은 의미가 없습니다. (0주차에 출산하는 경우는 없기 때문) 절편의 효과를 보고 싶다면 centering 기법을 사용하면 됩니다.

R 단순회귀분석

단순회귀모형 만들기 코드

weight <- c(64,75.3,73,82.1,76.2,95.7,59.4,93.4,82.1,78.9,76.7,82.1,83.9,73,64.4,77.6)
blood <- c(108,109,104,102,105,121,79,107,101,85,99,100,108,104,102,87)
model <- lm(blood~weight)
summary(model) # 회귀분석 결과 요약
anova(model) # anova table
confint(model) # 회귀계수의 confidence interval

summary(model)의 결과

Residuals: ## 잔차정보 요약

     Min       1Q   Median       3Q      Max 

-17.0962  -2.9774   0.8138   5.5072  13.4990 

Coefficients:  ## 회귀계수의 추정과 회귀계수의 유의성 검정

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

(Intercept)  61.8769    19.1890   3.225  0.00611 **

weight        0.5098     0.2462   2.070  0.05740 . 

F-statistic: 4.286 on 1 and 14 DF,  p-value: 0.0574      ## 모델이 유의한지에 대한 F-Test (분산분석)

anova(model)의 결과

Analysis of Variance Table

Response: blood

          Df Sum Sq Mean Sq F value Pr(>F)  

weight     1  368.8  368.80  4.2861 0.0574 .

Residuals 14 1204.6   86.05                 

• ANOVA F-test의 p-value와 회귀계수에 대한 t-test의 p-value가 같다. (0.0574)
• 또한 단순회귀분석에서는 모델이 유의한지에 대한 F-test의 p-value도 위의 p-value와 같다.

confint(model)의 결과

2.5 %     97.5 %

(Intercept) 20.7205466 103.033319

weight      -0.0183444   1.037845

• 회귀계수의 confidence interval을 구할 수 있다.

Prediction

newdata <- data.frame(weight=95)
predict(model, newdata, interval="predict")
predict(model, newdata, interval="confidence")

• weight=95일 때, blood의 값을 예측한다. 만약 여러개의 값을 예측하고 싶은경우는 weight=c(95,90,80, ...) 처럼 vector 형식으로 값을 주면 된다.
• interval=prediction으로 하면 예측구간이고 interval=confidence로하면 신뢰구간이다.

Prediction 결과

 fit      lwr      upr

1 110.3032 87.77959 132.8269

  fit      lwr      upr

1 110.3032 99.74413 120.8623

• fit이 예측값이고, lwr, upr은 각각 신뢰구간의 lower bound와 upper bound이다.

• 위와 같이 예측구간과 신뢰구간을 구할 수 있다. 예측구간이 신뢰구간에 비해 넓다. 예측구간은 y에 대한 95% interval이고 신뢰구간은 E(y|x) 에 대한 95% interval이다. 따라서 예측구간이 신뢰구간에 비해 넓은 것이다.

회귀분석의 가정

회귀분석에는 데이터에 대하여 어떠한 가정을 합니다. 그 가정하에서만 회귀분석은 유효하다고 할 수 있습니다. 가정이 틀리면 아무리 분석을 하더라도 그것이 타당하지 않게 되는 것입니다.

회귀분석의 가정에는 아래와 같은 것들이 있습니다.

1. 독립변수 X는 사전에 주어진 값으로 간주한다. 즉 독립변수는 미리 결정된 값이거나 고정된 값이다. 이를 비확률변수라고 한다.

2. 독립변수 X는 측정오차가 없다고 가정한다.

3. 독립변수 X 별로 Y의 부분모집단이 존재한다고 가정한다. 예를 들어 X=1일 때 Y의 부분모집단이 존재하며 이 부분모집단은 정규분포를 따른다.

4. Y의 부분모집단의 분산은 같고 이를 라고 한다. (등분산성)

5. Y의 부분모집단의 모평균은 일직선상에 있다. 이를 선형성의 가정이라고 한다.

6. Y는 서로 독립이다. 즉 e ~ N(0,) 라고 정리할 수 있다. (독립성)