R - Multinomial Logistic Regression

Multinomial Logistic Regression이란 y의 범주가 3개 이상(multi)이며 명목형(nomial)일 때 사용하는 로지스틱 회귀분석이다. 본 포스팅은 이 사이트의 튜토리얼을 보고 실습하고 해석한 자료이다.

데이터 설명

해당 데이터셋은 hsbdemo 데이터셋으로 총 n은 200이다.

ml <- read.dta("https://stats.idre.ucla.edu/stat/data/hsbdemo.dta")

> head(ml)
   id female    ses schtyp     prog read write math science socst       honors awards cid
1  45 female    low public vocation   34    35   41      29    26 not enrolled      0   1
2 108   male middle public  general   34    33   41      36    36 not enrolled      0   1
3  15   male   high public vocation   39    39   44      26    42 not enrolled      0   1
4  67   male    low public vocation   37    37   42      33    32 not enrolled      0   1
5 153   male middle public vocation   39    31   40      39    51 not enrolled      0   1
6  51 female   high public  general   42    36   42      31    39 not enrolled      0   1

Y : prog (프로그램 타입)

X : ses(social economic status), write (writing score)

다중회귀분석을 통해 X가 Y에 어떤 영향을 미치는지 알아본다.

빈도 테이블과 그룹 평균을 통해 X와 Y의 관계 알아보기

with(ml, table(ses, prog)) 
# 또는 table(ml$ses, ml$prog) 
# 또는 xtabs(~ses+prog, data=ml)

> with(ml, table(ses, prog))
        prog
ses      general academic vocation
  low         16       19       12
  middle      20       44       31
  high         9       42        7

with(ml, do.call(rbind, tapply(write, prog, function(x) c(M = mean(x), SD = sd(x)))))

> with(ml, do.call(rbind, tapply(write, prog, function(x) c(M = mean(x), SD = sd(x)))))
                M       SD
general  51.33333 9.397775
academic 56.25714 7.943343
vocation 46.76000 9.318754

nnet 패키지의 multinom 함수를 통해 다중 로지스틱 회귀분석 실시

ml$prog2 <- relevel(ml$prog, ref = "academic")
test <- multinom(prog2 ~ ses + write, data = ml)

multinom을 실행하기 전에 relevel이라는 함수를 통해 reference를 지정해준다. 이 reference를 기준으로 결과를 해석할 수 있다. 따라서 이 reference는 어떤 baseline이 되어야한다. 이 경우에 academic이 baseline이며 이를 기준으로 하여 분석한다.

multinom 함수를 실행하면 모델이 fitting되어 계수가 결정되고 이 정보는 test에 담긴다.

summary(test)

> summary(test)
Call:
multinom(formula = prog2 ~ ses + write, data = ml)

Coefficients:
         (Intercept)  sesmiddle    seshigh      write
general     2.852198 -0.5332810 -1.1628226 -0.0579287
vocation    5.218260  0.2913859 -0.9826649 -0.1136037

Std. Errors:
         (Intercept) sesmiddle   seshigh      write
general     1.166441 0.4437323 0.5142196 0.02141097
vocation    1.163552 0.4763739 0.5955665 0.02221996

Residual Deviance: 359.9635 
AIC: 375.9635 

이 예제의 경우 다중 로지스틱 회귀분석의 계수는 2세트가 나오게 된다.

1. academic vs general 의 log odds에 관한 계수

2. academic vs vocation의 log odds에 관한 계수

결과 해석

• write가 1단위 증가할 때, academic이 general이 될 log odds가 -0.058이 된다.
• write가 1단위 증가할 때, academic이 vocation이 될 log odds가 -0.1136037이 된다.
• ses:row에서 ses:high로 변할 때, academic이 general이 될 log odds가 -1.1628이 된다.
• ses:row에서 ses:middle로 변할 때, academic이 general이 될 log odds가 -0.5332가 된다.
• ses:row에서 ses:high로 변할 때, academic이 vocation이 될 log odds가 -0.9826이 된다.
• ses:row에서 ses:middle로 변할 때, academic이 vocation이 될 log odds가 +0.29가 된다.
  

변수의 유의성 판단

coefficient가 0라는 귀무가설 하에서 coef ~ N(0, sd^2)라는걸 이용하여 z-test 한다.

z <- summary(test)$coefficients/summary(test)$standard.errors
z

> z <- summary(test)$coefficients/summary(test)$standard.errors
> z
         (Intercept)  sesmiddle   seshigh     write
general     2.445214 -1.2018081 -2.261334 -2.705562
vocation    4.484769  0.6116747 -1.649967 -5.112689

# 2-tailed z test
p <- (1 - pnorm(abs(z), 0, 1)) * 2
p

> p <- (1 - pnorm(abs(z), 0, 1)) * 2
> p
          (Intercept) sesmiddle    seshigh        write
general  0.0144766100 0.2294379 0.02373856 6.818902e-03
vocation 0.0000072993 0.5407530 0.09894976 3.176045e-07

OR 추정값

exp(coef(test))

> exp(coef(test))
         (Intercept) sesmiddle   seshigh     write
general     17.32582 0.5866769 0.3126026 0.9437172
vocation   184.61262 1.3382809 0.3743123 0.8926116

회귀 계수들에 exponential을 취하면, 이는 OR에 대한 추정값이 된다.

Write가 평균일 때, ses에 따라 예측값이 어떻게 달라지는가?

dses <- data.frame(ses = c("low", "middle", "high"), write = mean(ml$write))
dses
predict(test, newdata = dses, "probs")

> dses <- data.frame(ses = c("low", "middle", "high"), write = mean(ml$write))
> dses
     ses  write
1    low 52.775
2 middle 52.775
3   high 52.775

> predict(test, newdata = dses, "probs")
   academic   general  vocation
1 0.4396845 0.3581917 0.2021238
2 0.4777488 0.2283353 0.2939159
3 0.7009007 0.1784939 0.1206054

Write를 변화시켜가면서 예측값의 변화 관찰하기

dwrite <- data.frame(ses = rep(c("low", "middle", "high"), each = 41), write = rep(c(30:70),
                                                                                   3))

## store the predicted probabilities for each value of ses and write
pp.write <- cbind(dwrite, predict(test, newdata = dwrite, type = "probs", se = TRUE))

## calculate the mean probabilities within each level of ses
by(pp.write[, 3:5], pp.write$ses, colMeans)

## melt data set to long for ggplot2
lpp <- melt(pp.write, id.vars = c("ses", "write"), value.name = "probability")
head(lpp)  # view first few rows

> head(lpp)  # view first few rows
  ses write variable probability
1 low    30 academic  0.09843588
2 low    31 academic  0.10716868
3 low    32 academic  0.11650390
4 low    33 academic  0.12645834
5 low    34 academic  0.13704576
6 low    35 academic  0.14827643

## plot predicted probabilities across write values for each level of ses
## facetted by program type
ggplot(lpp, aes(x = write, y = probability, colour = ses)) + geom_line() + facet_grid(variable ~
                                                                                        ., scales = "free")

• 이 그래프는 식을 전개하여 y를 sigmoid 함수 형태로 만든 것이다. 2개의 적합식에서 3개의 함수를 만들어낼 수 있다.

 Y~write의 결과랑 무엇이 다른가?

예측의 관점에서 보면 ses를 빼고 write만을 통해 Y를 예측하면 예측을 세분화하지 못하여 예측력이 떨어진다고 볼 수 있다. 이는 ses를 고려하지 않고, write와 Y의 관계만을 보기 때문이다. 따라서 ses가 추가됨으로써 write의 계수가 바뀐다는 것을 알 수 있다. 이 과정을 ses를 보정한다고도 말한다.

출처 - https://stats.idre.ucla.edu/r/dae/multinomial-logistic-regression/

R - 로지스틱 회귀분석

데이터 탐색

GRE, GPA, RANK이 입학(admission)에 어떤 영향을 주는지 로지스틱 회귀분석을 통해 분석한다.

library(aod)
library(ggplot2)

mydata <- read.csv("https://stats.idre.ucla.edu/stat/data/binary.csv")

# view the first few rows of the data
head(mydata)

# 데이터의 대략적인 분포 확인
summary(mydata)

# 데이터 구조 확인
str(mydata)

# 변수별 표준편차 확인
sapply(mydata, sd)

# contingency table : xtabs(~ x + y, data)
xtabs(~admit+rank, data=mydata)

> head(mydata)
  admit gre  gpa rank
1     0 380 3.61    3
2     1 660 3.67    3
3     1 800 4.00    1
4     1 640 3.19    4
5     0 520 2.93    4
6     1 760 3.00    2

> summary(mydata)
     admit             gre             gpa        rank   
 Min.   :0.0000   Min.   :220.0   Min.   :2.260   1: 61  
 1st Qu.:0.0000   1st Qu.:520.0   1st Qu.:3.130   2:151  
 Median :0.0000   Median :580.0   Median :3.395   3:121  
 Mean   :0.3175   Mean   :587.7   Mean   :3.390   4: 67  
 3rd Qu.:1.0000   3rd Qu.:660.0   3rd Qu.:3.670          
 Max.   :1.0000   Max.   :800.0   Max.   :4.000  

> str(mydata)
'data.frame':	400 obs. of  4 variables:
 $ admit: int  0 1 1 1 0 1 1 0 1 0 ...
 $ gre  : int  380 660 800 640 520 760 560 400 540 700 ...
 $ gpa  : num  3.61 3.67 4 3.19 2.93 3 2.98 3.08 3.39 3.92 ...
 $ rank : Factor w/ 4 levels "1","2","3","4": 3 3 1 4 4 2 1 2 3 2 ...

> sapply(mydata, sd)
      admit         gre         gpa        rank 
  0.4660867 115.5165364   0.3805668   0.9444602 

> xtabs(~admit+rank, data=mydata)
     rank
admit  1  2  3  4
    0 28 97 93 55
    1 33 54 28 12

GLM을 통한 로지스틱 회귀모형 구축

rank 변수를 factor 타입으로 변경시킨다.

glm 을통해 로지스틱회귀모형을 fitting시킨다. family='binomial' 인자를 통해, glm으로 로지스틱 회귀모형을 쓸 수 있다. (link function이 logit function이라는 의미)

mydata$rank <- factor(mydata$rank)
mylogit <- glm(admit ~ gre + gpa + rank, data = mydata, family = "binomial")

summary(mylogit)

> summary(mylogit)

Call:
glm(formula = admit ~ gre + gpa + rank, family = "binomial", 
    data = mydata)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.6268  -0.8662  -0.6388   1.1490   2.0790  

Coefficients:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept) -3.989979   1.139951  -3.500 0.000465 ***
gre          0.002264   0.001094   2.070 0.038465 *  
gpa          0.804038   0.331819   2.423 0.015388 *  
rank2       -0.675443   0.316490  -2.134 0.032829 *  
rank3       -1.340204   0.345306  -3.881 0.000104 ***
rank4       -1.551464   0.417832  -3.713 0.000205 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 499.98  on 399  degrees of freedom
Residual deviance: 458.52  on 394  degrees of freedom
AIC: 470.52

Number of Fisher Scoring iterations: 4

로지스틱 회귀모형 결과 해석

Call : 구축한 모형에 대해 다시 상기시켜준다.

Deviance Residuals : Deviance residual에 대한 정보를 알려주는데, model fitting이 잘 되었는지에 대한 measure이다. 이를 통해 모델이 잘 적합됐는지를 평가할 수 있다.

Coefficient : 회귀계수와 그것들의 표준편차, z-statistics(wals's z-statistics), p-value를 나타낸다. 위 결과에서는 모든 변수가 유의한 것을 알 수 있다. 로지스틱 회귀모형에서는 회귀계수가 변수가 한 단위 증가했을 때 log(odds)의 증가량으로 해석할 수 있다.

example

gre가 1증가할 때, admission의 log odds(non-admission에 대한)가 0.0022 증가한다.

gpu의 경우도 마찬가지로 해석한다.

더미변수인 rank의 경우 해석이 약간 다른데, 예를 들어 rank2의 회귀계수 -0.67은 rank1에서 rank2로 바뀌었을 때, log(odds)의 변화량이다. 즉, rank1에 비해 rank2가 admission에 안좋은 영향을 준다는 것을 알 수 있다.

회귀계수들의 신뢰구간 얻기

Log-likelihood를 통해 구하는 법

confint(mylogit)

> confint(mylogit)
Waiting for profiling to be done...
                    2.5 %       97.5 %
(Intercept) -6.2716202334 -1.792547080
gre          0.0001375921  0.004435874
gpa          0.1602959439  1.464142727
rank2       -1.3008888002 -0.056745722
rank3       -2.0276713127 -0.670372346
rank4       -2.4000265384 -0.753542605

주어진 회귀계수의 표준편차를 이용해 신뢰구간을 구하는 법

## CIs using standard errors
confint.default(mylogit)

> confint.default(mylogit)
                    2.5 %       97.5 %
(Intercept) -6.2242418514 -1.755716295
gre          0.0001202298  0.004408622
gpa          0.1536836760  1.454391423
rank2       -1.2957512650 -0.055134591
rank3       -2.0169920597 -0.663415773
rank4       -2.3703986294 -0.732528724

이 신뢰구간은 이렇게도 계산 가능하다.예를 들어, gre의 95% 신뢰구간은 아래와 같이 구한다.

> 0.002264-1.96*0.001094
[1] 0.00011976
> 0.002264+1.96*0.001094
[1] 0.00440824

Wald Test를 통해 범주형 변수의 overall effect 파악하기

rank 변수를 dummy 변수로 만들어, reference인 rank1과 비교해 어떤 effect가 있는지 분석하였다. 하지만 aod 패키지의 walt.test를 이용하면, rank의 overall effect를 파악할 수 있다. [b:회귀계수, Sigma:error term의 공분산행렬, Terms:rank변수가 있는 열] wald.test는 순서를 통해 해당 범주형 변수의 위치를 파악하기 때문에 순서를 잘 신경 써야 한다.

wald.test(b = coef(mylogit), Sigma = vcov(mylogit), Terms = 4:6)

> wald.test(b = coef(mylogit), Sigma = vcov(mylogit), Terms = 4:6)
Wald test:
----------

Chi-squared test:
X2 = 20.9, df = 3, P(> X2) = 0.00011

카이스퀘어 테스트의 p-value가 0.00011이므로 rank 변수는 유의하다.

변수간의 회귀계수가 같은지 검정하기

wald test를 통해 rank가 admission에 유의한 효과가 있다는 것을 파악하였고, rank1이라는 reference에 비해 rank2와 rank2가 유의한 영향을 준다는 것도 확인하였다. 한가지 의문은 rank2, rank3, rank4 의 회귀계수가 동일한지에 대한 것이다. 회귀 계수가 동일하다면, rank1일 때에 비해 rank2, rank3, rank4의 효과가 동일하다는 것을 의미한다.

이 구문은 rank2의 효과와 rank3의 효과가 동일한지를 검정한다.

l <- cbind(0, 0, 0, 1, -1, 0)
wald.test(b = coef(mylogit), Sigma = vcov(mylogit), L = l)

> l <- cbind(0, 0, 0, 1, -1, 0)
> wald.test(b = coef(mylogit), Sigma = vcov(mylogit), L = l)
Wald test:
----------

Chi-squared test:
X2 = 5.5, df = 1, P(> X2) = 0.019

검정 결과, rank2, rank3의 효과는 유의하다 다른 것으라 나타났다.

OR과 회귀계수의 관계

앞서 로지스틱 회귀분석에서의 회귀계수는 log odds의 증가량이라고 언급하였다. 그러면 회귀계수에 exponential을 취하면, OR의 증분이 된다.

exp(coef(mylogit))

> exp(coef(mylogit))
(Intercept)         gre         gpa       rank2       rank3       rank4 
  0.0185001   1.0022670   2.2345448   0.5089310   0.2617923   0.2119375 

해석 ex) gpa가 1증가하면 admission의 non-admission에 대한 OR이 2.23배 증가한다.

예측하기

1. 다른 변수들을 고정하고 rank가 변할 때 예측값의 변화를 보기

newdata1 <- with(mydata, data.frame(gre = mean(gre), gpa = mean(gpa), rank = factor(1:4)))

> newdata1
    gre    gpa rank
1 587.7 3.3899    1
2 587.7 3.3899    2
3 587.7 3.3899    3
4 587.7 3.3899    4

newdata1$rankP <- predict(mylogit, newdata = newdata1, type = "response") # probability 점추정값
newdata1

> newdata1
    gre    gpa rank     rankP
1 587.7 3.3899    1 0.5166016
2 587.7 3.3899    2 0.3522846
3 587.7 3.3899    3 0.2186120
4 587.7 3.3899    4 0.1846684

2. rank, gpa를 고정한 후, gre의 효과 보기

newdata2 <- with(mydata, data.frame(gre = rep(seq(from = 200, to = 800, length.out = 100),
                                              4), gpa = mean(gpa), rank = factor(rep(1:4, each = 100))))
newdata2

> newdata2
         gre    gpa rank
1   200.0000 3.3899    1
2   206.0606 3.3899    1
3   212.1212 3.3899    1
4   218.1818 3.3899    1
5   224.2424 3.3899    1
6   230.3030 3.3899    1
7   236.3636 3.3899    1
8   242.4242 3.3899    1
9   248.4848 3.3899    1
10  254.5455 3.3899    1

newdata3 <- cbind(newdata2, predict(mylogit, newdata = newdata2, type = "link", 
                                    se = TRUE))
newdata3 <- within(newdata3, {
  PredictedProb <- plogis(fit)
  LL <- plogis(fit - (1.96 * se.fit))
  UL <- plogis(fit + (1.96 * se.fit))
})

> head(newdata3)
       gre    gpa rank        fit    se.fit residual.scale        UL        LL PredictedProb
1 200.0000 3.3899    1 -0.8114870 0.5147714              1 0.5492064 0.1393812     0.3075737
2 206.0606 3.3899    1 -0.7977632 0.5090986              1 0.5498513 0.1423880     0.3105042
3 212.1212 3.3899    1 -0.7840394 0.5034491              1 0.5505074 0.1454429     0.3134499
4 218.1818 3.3899    1 -0.7703156 0.4978239              1 0.5511750 0.1485460     0.3164108
5 224.2424 3.3899    1 -0.7565919 0.4922237              1 0.5518545 0.1516973     0.3193867
6 230.3030 3.3899    1 -0.7428681 0.4866494              1 0.5525464 0.1548966     0.3223773

이 때 type="link" 옵션을 주면, link scale로 예측값을 내준다. 이는 logit 변환을 하기전의 예측값이다. 이를 통해 예측한 확률의 신뢰구간을 구할 수 있다. 왜냐하면 이 link scale의 예측값 (fit 변수)의 표준편차를 쉽게 알 수 있기 때문이다. 따라서 link scale의 예측값의 신뢰구간을 구한후 이를 다시 logit 변환하여 예측 확률의 신뢰구간을 얻을 수 있다. 이를 구현한 것이 바로 위의 코드이다.

plogis 함수를 link scale의 예측값에서 확률 추정값 p를 구해준다. 즉 log(p/1-p) = Xb 의 방정식을 풀어서 p를 구한다. 예를 들어 첫번째 라인의 경우, log(p/(1-p)) = -0.811 의 방정식을 양변에 exponential을 취해 풀면, 대략 p=0.308이 나온다.

newdata3의 예측값 시각화하기

rank를 고정시켜놓고 gre가 admission 확률에 어떤 영향을 주는지를 시각화

ggplot(newdata3, aes(x = gre, y = PredictedProb)) + geom_ribbon(aes(ymin = LL,
                                                                    ymax = UL, fill = rank), alpha = 0.2) + geom_line(aes(colour = rank),
                                                                                                                      size = 1)

출처

https://stats.idre.ucla.edu/r/dae/logit-regression/

Variant Calling

정의

variant call 은 어떤 개인의 genome혹은 transcriptome에서 reference와 뉴클레오타이드 서열의 차이가 있는지를 결정하는 것이다. 이 과정은 variant frequency와 confidence measure가 고려된다.

유즈케이스

DNA-seq : variants

• 질병과의 유전적 연관성
• 암에서의 돌연변이 알아내기
• 이형접한 인구 빈도 알아내기

RNA-seq : allele-speicifc expression

• Allelic imbalance
• Association with isoform usage (splicing QTLs)
• RNA editing (allele absent from genome)

Chip-seq : allele-specific binding

• protein binding이 allele-specific하게 일어나는 현상

핵심은 variant call을 통해 어떤 locus에 'variant가 있다' 라는 정보를 알 수 있기 때문에 이를 통해 위와 같은 연구를 할 수 있음

Variant Call은 Genotype 보다 더 일반적인 접근이다.

genotype은 어떤 locus의 allele의 집합을 나타낸다. 근데 배수성을 보통 2n으로 가정하고 한다.  또한 genotype은 주로 SNP에 대해서만 고려된다. 하지만 RNA-seq 같은 경우 2n이 아니며 genotype이 없다. 또 암유전체학 등에서는 copy number variation이 중요하기 때문에 variantl call이 genotype보다 일반적인 접근법이라 할 수 있다.

참고

https://www.bioconductor.org/help/course-materials/2014/CSAMA2014/3_Wednesday/lectures/VariantCallingLecture.pdf

파일첨부

VariantCallingLecture.pdf

참고할 사이트들

Best Practices workflow

https://software.broadinstitute.org/gatk/best-practices/bp_3step.php?case=GermShortWGS

범주형 자료분석 - 맥니마 검정

맥니마 검정은 1947년 Psychometrika 이라는 학술지에 처음 발표된 것으로 미국의 심리학자이자 통계학자인 Quinn McNemar에 의해 만들어졌습니다. 맥니마 검정이 어떤 상황에서 쓰이고 무엇을 검정하기 위한 것인지 알아보겠습니다.

맥니마 검정은 짝지은 명목형 데이터에서 Column과 Row의 marginal probability가 같은지를 검정하는데에 쓰입니다. 2x2 테이블일 때 적용가능한 방법으로, 즉 변수가 2개일 때 쓸 수 있는 방법입니다. '짝지은' 이라는 말은 두 변수의 값이 서로 연관되어있으며, (동일한 사람에 대해 두 번 측정하거나, 부모-자식 처럼 관련있는 사람들에 대해 측정함) 두 변수의 값을 측정할 때 총 n수가 변하지 않았다는 의미로 볼 수 있습니다. 예를 들어 다음의 데이터를 봅시다.

표1

위 데이터의 161명의 환자를 대상으로 약을 처방하기 전의 질병의 상태(present, absent)와 처방한 후의 질병의 상태를 나타내는 집계 데이터입니다. 이 집계 데이터의 row 데이터는 아래와 같이 생겼을 것입니다.

이러한 모양의 행이 161개인 데이터를 집계한 결과가 위의 2x2테이블인 셈입니다. 짝지었다(paired)는 말은 이 경우에는 동일한 사람에 대해 Before와 After을 측정하였으니 이를 짝지은 데이터라고 할 수 있습니다. 또 Before와 After이 범주형이기 때문에 맥니마 검정은 범주형 자료분석의 범위에 속해있는 것입니다. 만약 Before와 After의 변수의 값이 연속형이라면 paired t-test를 통해 before와 after의 모평균의 차가 0인지를 검정할 수 있습니다.

맥니마 검정에서의 귀무가설과 대립가설은 무엇인가?

맥니마 검정에서 검정하고자하는 것은 marginal probability가 같은지를 검정하는 것입니다. 즉, 식으로 p_a + p_b = p_a + p_c , p_c + p_d = p_b + p_{d 으로 표현할 수 있습니다. (이 때 p(a) = a/(a+b), p(b) = b/(b+c) ... 입니다. ) 식을 정리하면 결국, 귀무가설은} p_{b = pc가 됩니다. 이 귀무가설을 말로 표현해보면 "전체에서 Before:present의 비율과, After:present의 비율이 같은가?" 또한 "Before:absent와 After:present의 비율이 같은가?" 로 표현할 수 있습니다. 공교롭게도 이것은 하나의 식으로 표현가능합니다. 바로 pb = pc 이죠. 이 케이스에서는 약 처방 후 환자의 상태가 유의하게 변했는가? 라고 표현할 수도 있습니다.}

_{귀무가설과 대립가설}

검정통계량

(자유도 = 1)

맥니마 검정의서의 귀무가설과 대립가설, 검정통계량은 위와 같습니다. 이해를 돕기 위해 극단적인 케이스를 살펴보겠습니다.

표2

데이터가 이렇게 생겼다고 해봅시다. 이 경우에 환자는 약의 처방과 관계 없이 이전과 이후의 상태가 완벽하게 같습니다. 따라서 Pb와 Pc의 추정량이 같으며, 검정통계량은 0이 됩니다. p-value를 계산하면 1이되며, 귀무가설을 기각할 수 없게됩니다.

표1

다시 원래의 표로 돌아가 검정 통계량을 계산해보겠습니다. b=6, c=16이므로 검정통계량은 100/22 = 4.55입니다. p-value는 카이제곱분포에서 4.55이상의 값을 관측할 확률이므로 0.033으로 계산할 수 있습니다. 따라서 marginal probability가 같다는 귀무가설을 기각할 수 있습니다.

참고

http://www.statisticssolutions.com/non-parametric-analysis-mcnemars-test/

http://www.statisticshowto.com/mcnemar-test/

https://en.wikipedia.org/wiki/McNemar%27s_test

범주형 자료에서 RR, OR을 구할 때, 혼란 변수를 보정하는 한 가지 방법 - 코크란 멘텔 헨젤 방법

앞서 혼란변수가 변수간의 관계를 추론할 때 어떠한 문제를 일으키는 지 알아보았습니다. (http://3months.tistory.com/228) 이러한 혼란 변수는 다양한 방법으로 보정할 수 있습니다. 이번 포스팅에서는 2x2 테이블에서 Relative Risk와 Odds Ratio를 구할 때, 혼란 변수를 보정하여 구하는 방법인 코크란 멘텔 헨젤 방법에 대해 알아보겠습니다.

OR = Odds(Obes)/Odds(Not Obese) = (46/254) / (60/640) = 1.93

나이가 혼란변수로 의심되 나이에 따라 Stratification 합니다.

Age < 50

Age >= 50

Age < 50 : OR = (10/90) / (35/465) = 1.476

Age >= 50 : OR = (36/164) / (25/175) = 1.53

Stratification 후에 각각 계산된 OR 값을 통해 혼란 변수로 인해 OR 값이 과추정 되었음을 알 수 있습니다. 이를 감안하여 두 OR 값을 합치기 위한 방법이 바로 코크란-멘텔-헨젤 방법을 통해 공통 OR을 구하는 것입니다.

RR과 OR에 대한 코크란-멘텔-헨젤 추정량은 위와 같습니다. 위 공식에 넣어 공통 OR을 계산하면 1.52를 얻을 수 있습니다.

위에서는 나이 50살을 기준으로해서 단지 2개의 그룹으로 나누었습니다. 하지만 이 기준은 모호하고, 그룹을 세분화하지 않았기 때문에 아직 confounding effect가 남아있을 수 있습니다. 이를 더욱 제거하기위해 아래와 같이 더 세분화하여 그룹을 나눌 수도 있습니다. 이 경우에도 마찬가지로 코크란-멘텔-헨젤 방법의 공식을 통해 공통 OR, RR을 구할 수 있습니다.

또한 혼란 변수가 보정 된 코크란 멘텔 헨젤 추정량을 얻기 전 OR 또는 RR에 대한 동질성 검정(homogeneity test)을 실행하는 것도 좋습니다. 위 경우에는 두 그룹의 OR 값이 비슷했기에 이를 합쳐 공통 OR을 구하는 것이 합리적인 방법이였지만, 만약 두 그룹의 OR의 차이가 많이 난다면, 따로따로 관계를 해석하는 것이 바람직할 것입니다. 이를 테스트하기위해 사용하는 한가지 방법이 breslow-day test입니다. 이 검정의 p-value가 작으면, 두 그룹의 OR 값이 유의하게 다른것이며, 이 경우에는 코크란 멘텔 헨젤 추정량을 이용하기보다 그룹별로 따로 구한 OR을 각각 해석하는 것이 좋습니다.

참고

http://sphweb.bumc.bu.edu/otlt/mph-modules/bs/bs704-ep713_confounding-em/BS704-EP713_Confounding-EM7.html

심슨의 역설(Simpon's Paradox)

심슨의 역설로 유명한 영국 통계학자 에드워드 심슨(1922~)

통계학에서의 심슨의 역설이란 그룹을 나누어서 봤을 때는 나타나는 경향이 그룹을 합쳤을 때는 사라지거나 오히려 경향이 역전되는 것을 말합니다. 이는 경향이 역전된다니 그것이 무슨 말일까요?

예를 들어, 신장 결석을 치료하는 두 가지 수술 방법이 있다고 해봅시다. 수술방법 A는 가격이 비싼 대신 치료 성공률이 높습니다. 하지만 가격이 비싸 기 때문에 결석이 큰 중증 환자만 수술방법 A를 택합니다. 반면, 수술 방법 B는 가격이 싸서 비교적 증상이 약하고 결석의 크기가 작은 환자들이 많이 하는 수술방법입니다. 하지만 가격이 싼 만큼 A에 비해서는 치료 성공률도 낮습니다. 이 때, 신장 결석의 크기 별로 수술방법 A와 수술방법B의 치료율을 알아보겠습니다.

작은 결석

수술방법 A의 치료율 = 81/87 = 93%

수술방법 B의 치료율 = 234/270 = 87%

큰 결석

수술방법 A의 치료율 = 192/263 = 73%

수술방법 B의 치료율 = 55/80 = 69%

결과를 보면 작은결석과 큰 결석 모두에서 수술방법 A의 치료율이 조금씩 높습니다. 그러므로 결론적으로 수술방법 A의 치료율이 높습니다.

작은 결석 + 큰결석

수술방법 A의 치료율 = 273/350 = 78%

수술방법 B의 치료율 = 289/350 = 83%

하지만 결석의 크기에 따라 치료율을 알아보지 않고 합쳐서 봤을 때는 어떨까요? 이 때, 수술방법 A의 치료 성공률은 273/350 = 0.78, 수술방법 B의 치료 성공률은 289/350 = 0.83입니다. 합치고 보니 오히려 수술방법 B가 더 좋은 것입니다. 이처럼 그룹(결석)을 나누었을 때, 보였던 경향(수술방법A의 치료율이 더 좋다)이 그룹을 합쳤을 때는 오히려 역전(수술방법 B의 치료율이 더 좋다)되는 것을 심슨의 역설이라고 합니다.

그렇다면 왜 이런 현상이 나타날 것일까요? 바로 수술을 받는 환자들의 결석 크기의 분포가 다르기 때문입니다. 서두에도 말했든 수술방법 A는 중증 환자가 받는 수술이기 때문에 큰 결석을 갖고 있는 환자가 많이 받았고, 수술 방법 B는 작은 결석 환자가 많이 받았습니다. 그래서 수술방법 A는 대부분 치료가 어려운 환자를 상대했기 때문에 치료율이 낮을 수 밖에 없는 것입니다.

이를 변수를 도입하여 모델링 해봅시다. 수술방법을 X1, 결석 크기를 X2, 치료 여부를 Y라고 놓아봅시다. 우리가 알고 싶은 것은 X1과 Y의 관계입니다. 하지만 X2는 X1과 Y에 모두 영향을 주면서, (즉 X2<->X1과 X2<->Y가 독립이 아닙니다.) X1와 Y의 인과관계 판단에 영향을 줍니다. 이러한 변수 X2를 우리는 혼란변수(confounding variable)이라고 합니다. 

심슨의 역설은 데이터를 해석할 때는 주의를 기울여야 한다는 교훈을 주었습니다. 왜냐하면 이 경우에는 작은 결석과 큰 결석으로 나눈 데이터를 먼저 보았기 때문에 이러한 결론을 내릴 수 있지만, 실제로 우리가 보는 데이터가 그룹별로 나뉘어진 데이터가 아닐 수도 있기 때문입니다. 그럴 경우에 우리는 '수술방법 B가 더 낫다' 라는 잘못된 결론을 내릴 수도 있습니다. 따라서 데이터를 해석하고 그 안에 숨겨진 인과관계를 추론할 때는 이러한 혼란 변수를 반드시 고려하여 결론을 내려야 합니다.

참고

https://en.wikipedia.org/wiki/Simpson%27s_paradox

원래 데이터

해당 지역의 10000명당 기관의 수를 나타낸 데이터입니다. (출처-통계청)

           sgg_kor facility_type 2003fac 2004fac 2005fac 2006fac 2007fac
 4               계  상급종합병원  0.0086  0.0086 0.00852  0.0087  0.0087
6               계      종합병원  0.0494  0.0494 0.05054  0.0510  0.0530
8               계          병원  0.1645  0.1755 0.18450  0.1937  0.2127
10              계      요양병원  0.0139  0.0231 0.04120  0.0727  0.1200
14              계      치과병원  0.0211  0.0221 0.02517  0.0274  0.0311
18              계      한방병원  0.0311  0.0320 0.03024  0.0292  0.0288
22              계        조산원  0.0143  0.0129 0.01055  0.0103  0.0104
24              계        보건소  0.0473  0.0475 0.04750  0.0478  0.0481
26              계      보건지소  0.2601  0.2605 0.25859  0.2571  0.2602
28              계    보건진료소  0.3865  0.3881 0.38504  0.3843  0.3875
30              계    보건의료원  0.0035  0.0035 0.00345  0.0034  0.0035
36            서울  상급종합병원  0.0195  0.0195 0.01942  0.0193  0.0196
38            서울      종합병원  0.0418  0.0409 0.04079  0.0406  0.0412
40            서울          병원  0.1158  0.1168 0.12431  0.1304  0.1432
42            서울      요양병원  0.0029  0.0058 0.01068  0.0251  0.0461

Melt

Melt는 Column을 Row로 바꾸는 함수입니다. (Key=id에 지정, 나머지값들은 Row로 변환)

facility <- melt(facility, id=c("sgg_kor", 'facility_type'))

sgg_kor facility_type variable  value
 1      계  상급종합병원  2003fac 0.0086
 2      계      종합병원  2003fac 0.0494
 3      계          병원  2003fac 0.1645
 4      계      요양병원  2003fac 0.0139
 5      계      치과병원  2003fac 0.0211
 6      계      한방병원  2003fac 0.0311

reshape

Reshape함수를 통해 Row를 Column으로 바꿀 수 있습니다. (Key를 idvar로 놓고, 컬럼으로 바꾸고 싶은 변수를 timevar에 지정)

facility <- reshape(facility, idvar=c("sgg_kor", "variable"), timevar="facility_type", direction="wide")

  sgg_kor variable value.상급종합병원 value.종합병원 value.병원 value.요양병원 value.치과병원      
계         2003fac             0.0086         0.0494     0.1645         0.0139         0.02      
    서울  2003fac             0.0195         0.0418     0.1158         0.0029         0.0350         
  강남구  2003fac             0.0373         0.0560     0.2425             NA         0.1679         
  강동구  2003fac             0.0209         0.0417     0.1043             NA             NA         
  강서구  2003fac                 NA         0.0185     0.1668         0.0185         0.0185         
  관악구  2003fac                 NA         0.0190     0.0949             NA             NA             

Allele and Genotype Frequency

하디 바인베르크의 법칙

특정 조건 안에서 allele frequency와 genotype frequency는 영속적으로 보존된다.

위 그림과 같이 allele frequency가 freq(A) = p, freq(a) = q일 때, 다음 세대의 genotype frequency의 기댓값은 freq(AA) = p^2, freq(Aa)=2pq, freq(aa)=q^2 이다. 이 genotype frequency를 통해 다시 allele frequency를 계산하면 freq(A)=p, freq(a)=q가 나온다. 즉, allele frequency와 genotyep frequency는 계속해서 유지된다.

allele frequency와 genotype frequency에서 중요한점

1. genotype frequency를 알면 allele frequency를 알 수 있다.

2. allele frequency를 알아도 genotype frequency를 반드시 알 수는 없다.

2번의 경우 예를 들어, freq(A) = 0.5, freq(a) =0.5 라하자. 그러면 하디 바인베르크 법칙을 만족하면 genotype frequency는 freq(AA)=0.25, freq(Aa)=0.5, freq(aa)=0.25일 것이다 하지만 freq(AA)=0.5, freq(Aa)=0, freq(aa)=0.5 여도 주어진 조건에 맞는다. 따라서 allele frequency를 알아도 genotype frequency를 반드시 알 수 있는 것은 아니다.allele frequency는 재료이다. 그것이 어떻게 조합되어 genotype을 구성할지는 확정적이지 않다. 이는 하디 바인베르크 법칙에 위배된다.

그렇다면 하디 바인베르크 법칙은 언제 만족하는가?

1. random mating

: 위에서 언급한 allele frequency를 알아도 genotype frequency를 반드시 알 수 없다. 하지만 random mating이라면 genotype frequency의기댓값은 하디 바인 베르크에서의 값과 같다.

2. no mutation, selection, migration

: 모집단의 임의적 변화가 없어야한다.

3. 무한한 모집단 사이즈

: genetic drift가 없어야 한다.

하디 바인베르크는 위 가정들을 만족할 수 없다. 그렇다면 하디 바인베르크 법칙은 왜 필요할까? 하디 바인베르크 법칙이 중요한 이유는 그것이 귀무가설(Null hypothesis)이 될 수 있기 때문이다. 하디 바인베르크 법칙을 귀무가설로 놓고 진화적 과정이 어떻게 일어나는지를 살펴볼 수 있다.

하디 바인베르크 법칙이 깨지는 상황

하디 바인베르크 법칙은 어떠한 집단 내에서 random mating이 일어나지 않을 때이다. 즉, 인종이 섞인 집단을 하나의 집단으로 놓았을 때, 하디 바인베르크 법칙이 깨지게 된다. 예를 들어, 한국인, 중국인이 각각 하디 바인베르크 법칙을 만족하더라도, 한국인, 중국인을 합쳐서 하나의 집단을 만들고 이 집단에 대해 하디 바인베르크 법칙을 만족하는지 테스트하면 하디 바인베르크 법칙을 만족하지 않는 결과가 나올 수도 있다.

wahlund effect

인종을 합쳤을때 heterozygous genotype이 하디-바인베르크 법칙으로 구한 hetero genotype의 기댓값보다 적게 나타나는 현상.

하디 바인베르크의 법칙이 중요한 이유는 GWAS를 할 때, 연구집단이 하디-바인베르크 법칙을 만족한다는 것을 가정하기 때문이다. 하디-바인베르크 법칙이 만족되지 않으면 GWAS의 결과가 부정확할 수 있다.

예를 들어, 인종1과 인종2의 SNP의 allele frequency가 다르고, 인종1의 유병률이 인종2의 유병률보다 높다고 하자. (유전적 원인이 아닌 환경적 원인에 의해 인종1의 유병률이 높다.) 그러면 인종1이 해당 마커에 많이 갖고있는 allele 근처의 gene이 질병에 영향을 준다는 잘못된 결론을 내릴 수 있다. 그러므로 이를 방지하기 위해 하디-바인베르르 법칙을 만족하는지 테스트를 해야한다. 연구집단에 대해 하디-바인베르크 법칙을 만족한다면 그 연구집단내에서 random mating이 이일어난다는 것을 알 수 있고, 그 연구집단이 하나의 인종을 이룰 수 있다는 것을 알 수 있기 때문에 다른 환경적 요인을 보정하여 bias를 방지할 수 있다. 

Gene vs Environment

어떠한 사람의 특성이 유전자로부터 비롯된 것인지, 환경으로부터 비롯된 것인지를 아는 것은 중요하다. 예를 들어, 만약 그 특성이 질병이라면 환경을 의식적으로 컨트롤함으로써 질병을 예방할 수 있다.

이러한 것을 확인하는 한가지 방법은 쌍둥이를 연구하는 것이다. 유전 vs 환경을 연구할 때, 유전, 환경 둘 중 하나를 고정시키고 다른쪽의 effect만 보면 보다 정확한 유전, 환경 효과를 파악할 수 있다. 쌍둥이 연구는 유전, 환경 중 하나를 고정하고 다른 요인의 효과를 파악하는데 적합하다.

유전적 요소 파악하기

쌍둥이를 통한 연구에서 보통 쌍둥이는 shared environment라고 가정한다. 이것이 쌍둥이 연구의 이점이다. 일란성 쌍둥이에서의 상관도(correlation)이 이란성 쌍둥이에 비해 높다면 그 특성은 유전적 요소가 환경적 요소에 비해 더 많이 개입할 것이라고 추론할 수 있다.

실제 예로 IQ의 경우 일란성 쌍둥이의 상관계수는 0.85, 이란성 쌍둥이의 상관계수는 0.42이다.

환경적 요소 파악하기

쌍둥이를 통한 연구에서 어떠한 특성의 환경적 요소를 파악하는 방법은 함께 살은 일란성 쌍둥이와 떨어져서 살은 일란성 쌍둥이를 비교하는 것이다. 함께살은 쌍둥이는 유전, 환경 모두 동일하지만, 떨어져 살은 쌍둥이는 유전, 환경 중에 유전만 동일하고 환경이 다르다. 따라서 함께 살은 쌍둥이의 상관계수가 더 높을수록 그 특성은 환경적 요소가 크게 작용한다고 생각할 수 있다.

보통 쌍둥이간의 상관계수는 같이 살은 일란성 > 떨어져 살은 일란성, 같이 살은 일란성 > 같이 살은 이란성 이 성립한다.

예를 들어서, 언어적 능력이라는 특성에 대하여 상관계수가 다음과 같이 나온 경우를 보자.

같이 살은 일란성(유전 같음, 환경 같음) : 0.76

떨어져 살은 일란성(유전 같음,환경 다름) : 0.51

같이 살은 이란성(유전 다름, 환경 같음) : 0.43

이 경우 0.76-0.51을 환경적 요소에 의한 특성의 차이로 볼 수 있고, 0.76-0.43을 유전적 요소에 의한 특성의 차이로 볼 수 있다.

유전, 환경 정량화하기

어떠한 정량화된 특성의 분산을 유전에 의한 분산과 환경에 의한 분산으로 나눌 수 있다.

simple formula : V = V(환경) + V(유전)

V는 분산을 구하는 공식으로 계산할 수 있는데, 어떻게 V(환경), V(유전)을 계산할 수 있을까? 방법은 생각보다 간단하다. V(환경)만 구하면, V(유전)도 구할 수 있고, V(환경)은 genotype을 고정시킨 후의 특성의 분산을 통해 구할 수 있다. 어떠한 특성에 관련있는 유전자가 6개라고 하자. 이 6개의 유전자에 대한 genotype이 모두 같은 sample안에서 분산을 구하면 이 분산은 오직 환경으로 인한 분산이다. 따라서 이 때의 분산 V1 = V(환경)이다. 이를 기존에 랜덤하게 뽑은 sample의 V에서 빼면 V(유전)을 구할 수 있다. 이 때 Heritability는 아래와 같이 정의된다.

Heritability = V(유전) / V(유전)+V(환경)

<모든 유전자가 heterozygous인 F1세대를 교배하여 만든 F2 세대를 통해 heritability 계산>

하지만 이런 방법을 사람에게 적용할 수는 없다. 실제로 사람에 대해 유전, 환경적 요소를 정량화하는 방법을 알아보자.

부모-자손 상관계수를 통해 Heritability 계산하기

실제 사람을 대상으로 heritability를 계산하는 방법 중 하나는, 부모-자손 상관계수를 이용하는 방법이다. 하지만 이 경우에 많은 한계점이 있다. 우선 첫 번째로, 부모-자손 상관계수를 이용하여 Heritability를 추정하는 경우 overestimated될 가능성이 있다. 부모-자손은 환경을 많은 부분 공유하기 때문이다. 우리가 보고 싶은건 유전적 요소 뿐인데 이러한 방법을 이용하게 되면 환경이 개입하여 heritability를 실제보다 더 높게 추정하게 된다. 또한 V(유전)을 알고 싶은 경우, 환경을 고정시킨 후 보아야하는데 환경을 고정시킬 수가 없다. 환경은 인구 집단에 따라 인종에 따라 다르다. 즉, V(환경)이 일정하지 않다. 그래서 연구 집단으로 어느 집단을 선택하냐에 따라 heritability 추정값이 달라진다. 

<키가 유전적 요소에 의해서만 좌우되는 경우>

<실제 키의 부모-자손 상관계수>

Breeder's Equation

인위 선택을 통해 Heritability를 계산할 수 있다. 옥수수 키의 평균이 5인치라고 하자. 이 방법에서는 7인치인 옥수수를 교배해서 나온 옥수수의 평균이 7인치라면 heritability = 1이다. 7인치 옥수수를 교배해서 나온 옥수수의 평균이 5인치면, 유전의 영향을 전혀받지 않는 것으로 heritability = 0이다. 만약 평균이 6인치가 나왔으면 1/2 = 0.5이다. 6.5인치면 1.5/2 = 0.75이다. 왜냐하면 7인치 옥수수를 뽑은 것은 모집단에서 뽑은 것이기 때문에 7-5=2 = V(유전)+V(환경)이다. 하지만 7인치를 교배해서 나온 6.5인치에서 평균을 뺀 6.5-5=1.5=V(유전) 이다. 따라서 이 둘의 비율로 heritability를 추정할 수 있다. 인위선택뿐 아니라 자연선택에서도 이 논리를 그대로 적용할 수 있다.

나눔바름고딕 CSS 적용

폰트 정의

@font-face { font-family: 'NanumBarunGothic';
src: url('/fonts/NanumBarunGothic.eot');
src: url('/fonts/NanumBarunGothic.eot') format('embedded-opentype'),
url('/fonts/NanumBarunGothic.woff') format('woff');}

폰트사용시

body {font-family: 'NanumBarunGothic', 'serif';}

SRC 부분은 웹프로그램 소스상에서 폰트 파일들이 어디에 위치하는지를 적어줍니다.

NanumBarunGothic.zip

나눔고딕 CSS 적용

폰트 정의

@font-face { font-family: 'NanumGothic';
src: url('/fonts/NanumGothic.eot');
src: url('/fonts/NanumGothic.eot') format('embedded-opentype'),
url('/fonts/NanumGothic.woff') format('woff');}

폰트사용시

body {font-family: 'NanumGothic', 'serif';}

NanumGothic.eot

NanumGothic.woff

NanumGothicBold.eot

NanumGothicBold.woff