SAS 코드

data corr;
input method1 method2;
cards;
132 130
138 134
144 132
146 140
148 150
152 144
158 150
130 122
162 160
168 150
172 160
174 178
180 168
180 174
188 186
194 172
194 182
200 178
200 196
204 188
210 180
210 196
216 210
220 190
220 202
;run;

상관계수 출력 및 모상관계수=0을 검정

proc corr;
var method1 method2 ;
run;

상관계수 출력 및 모상관계수=0.98을 피셔의 방법으로 검정

proc corr fisher(rho0=0.98);
var method1 method2 ;
run;

크루스칼-왈리스 검정은 모수적 방법에서의 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이다. 따라서 예측구간이 신뢰구간에 비해 넓은 것이다.