博學之 審問之 愼思之 明辨之 篤行之

  다른 연구 자료의 Survival curve 로 좀 고민하다가 찾은 내용이다.(http://www.cbgstat.com/methods/R_maximal_chi/R_maximal_chi.htm)

  특정 조건에 따른 생존 곡선을 작성할 때 cutoff value 값을 정해야 할 필요가 있는데 친절하게도 p-value 가 가장 작아지는 범위를 골라준다. 이 함수의 이름과는 좀 다르다. :)

  maxstat 이라는 라이브러리를 설치해야 한다. survival 라이브러리는 일반적으로 기본적으로 설치되어 있다.

abc <- read.table("GBM_data.csv", header=TRUE, sep=",")

max <- maxstat.test(Surv(Survival, Death) ~ Age, data=abc, smethod="LogRank", pmethod="HL")
다양한 옵션이 있지만 최소한의 옵션 정도는 위와 같다.

Surv(Survival, Death) ~ Age 
  Survival 은 생존 기간 정보가 들어간 항목, Death 는 특정 event가 발생한 항목, Age는 cutoff 값을 알고 싶은 항목이다. Death 부분에는 기본적으로 0은 생존, 1은 event 발생이 들어가게 된다.

smethod는 어떤 방식에 의하여 검증 할 것인가, pmethod는 p.value 검정 방법이다. 이 두 항목은 잘 모르겠다. :)

이런식으로 결과 값이 나온다. estimated cutpoint 에 나오는 값이 cutoff 값이다. 아직 다른 항목으로 확인을 못했봤지만 작은 값부터 시작하여서 p-value 를 찾는 것 같으며, cutpoint 보다 같거나 작은 값과 그 보다 큰 값으로 분류가 되어 있다. 그리고 문제는 저 p-value 인데 이 통계 방법에서 나오는 p-value 와 SPSS에서 나오는 p-value 는 다르게 나온다. 해결하는 방법은 따로 한 번 더 구하면 되는 것이다. 그 과정은 약간의 자료값을 처리해야 한다. 이제부터는 cutpoint 가 정해진 age 값이 있기 때문에 그 보다 같거나 작은 값과 큰 값을 가지는 값을 새롭게 만들어야 한다.

a <- matrix(nrow=22, ncol=3)
우선 새로운 matrix 를 만들어야 하는데, 22는 증례 숫자, 3은 Survival, Event, Age(cutpoint에 의한 재분류) 값을 위한 것이다.

for (X in 1:22){a[X,1] <- abc[X,2];a[X,2] <- abc[X,3]}
첫 번째 column 은 survival 값이 들어가고, 2번째 column 은 Event 정보가 있는 값이 들어간다.

for (X in 1:22){if (abc[X,1] <= max$estimate) {a[X,3] <- 1} else {a[X,3] <- 2}}
abc라는 자료의 첫 번째 column 에는 나이가 있었는데 그 값이 이번에 구한 cutpoint 인 5보다 같거나 작으면 그 항목에는 1이라는 값을 그렇지 않으면 2라는 값을 입력하도록 했다.

plot(survfit(Surv(a[,1],a[,2]) ~ a[,3]), lty=4:5)
이렇게 하면, survival curve 를 구할 수 있고

survdiff(Surv(a[,1],a[,2]) ~ a[,3])
이렇게 하면, p-value를 알 수가 있다.

Creative Commons License

이 저작물은 크리에이티브 커먼즈 코리아 저작자표시-비영리-변경금지 2.0 대한민국 라이선스에 따라 이용하실 수 있습니다.

  오늘 아침에 Journal 발표 시간에 갑자기 생각난 아이디어를 바탕으로 heatmap 을 다시 만들어 보았다. 만들면서 느낀 건데 이런 chip 기반의 data는 해석하는데 항상 주의가 필요하다. 원 데이터를 조작하지 않아도 중간의 과정을 저자가 원하는 방향으로 이끌게 되면 그렇게 보이는 자료를 만들 수 있다. :) 이 heatmap 은 9, 10 번과 그 이외의 자료가 유의한 차이를 나타나도록 만든 것이다. 하지만, 마음만 먹으면 다른 자료가 유의하게 나오도록 한 다음 다른 설명을 하는 것도 가능하다. -_-

  그림 만드는 것은 어느정도 해결이 되었으니 이제는 원자료를 해석하는데에 집중해야겠다.

  그렇게 깔끔한 구문이라고는 못하지만 문법적인 문제는 없는 것 같고, 컴퓨터도 그렇게 느린 편은 아니라서 좀 길게 만들어도 신경을 쓰지 않았다.

2010/05/24 - [공부해 봅시다/R-Project] - Green-Black-Red heatmap
2010/05/26 - [공부해 봅시다/R-Project] - Array data 살펴 보기 - Simple


library(gplots)
data <- read.table("Book1-1-3.csv", header=T, sep=",")

CONF <- 0.0001 #유의수준

COUNT <- 0  #유의한 값을 가지는 항목 숫자세기
for (X in 1:5467) {
x <- c(data[X,2],data[X,3],data[X,4],data[X,5], data[X,6],data[X,7],data[X,8],data[X,9])
y <- c(data[X,10],data[X,11])
z1 <- t.test(x,y, var.equal=FALSE, paired=FALSE, conf.level=CONF)     
z2 <- t.test(x,y, var.equal=TRUE, paired=FALSE, conf.level=CONF)
if (z1$p.value < CONF | z2$p.value < CONF) {COUNT <- COUNT +1}
}

newdata <- matrix(nrow=COUNT, ncol=10) #유의한 값을 가지는 것을 위한 새로운 matrix 만들기
COUNT1 <- 0
for (X in 1:5467) {
x <- c(data[X,2],data[X,3],data[X,4],data[X,5], data[X,6],data[X,7],data[X,8],data[X,9])
y <- c(data[X,10],data[X,11])
z1 <- t.test(x,y, var.equal=FALSE, paired=FALSE, conf.level=CONF)     
z2 <- t.test(x,y, var.equal=TRUE, paired=FALSE, conf.level=CONF)
if (z1$p.value < CONF | z2$p.value < CONF) {COUNT1 <- COUNT1 +1
newdata[COUNT1, 1] <- data[X,2];newdata[COUNT1, 2] <- data[X,3];newdata[COUNT1, 3] <- data[X,4];newdata[COUNT1, 4] <- data[X,5];
newdata[COUNT1, 5] <- data[X,6];newdata[COUNT1, 6] <- data[X,7];newdata[COUNT1, 7] <- data[X,8];newdata[COUNT1, 8] <- data[X,9];
newdata[COUNT1, 9] <- data[X,10];newdata[COUNT1, 10] <-data[X,11]
}
}

my.hclust <- function(x){hclust(x, method="complete")}
my.dist <- function(x){dist(x, method="euclidean")}
 
heatmap.2(newdata, col=redgreen, trace="none", density.info="none", hclustfun=my.hclust, distfun=my.dist, main=CONF)

Creative Commons License

이 저작물은 크리에이티브 커먼즈 코리아 저작자표시-비영리-변경금지 2.0 대한민국 라이선스에 따라 이용하실 수 있습니다.

for (X in 1:5467) {
x <- c(data[X,2],data[X,3],data[X,4],data[X,5], data[X,6],data[X,7],data[X,8],data[X,9])
y <- c(data[X,10],data[X,11])                       
if (t.test(x,y)$p.value < 10e-07) print(data[X,])
}
  예전에 사놓고 존재를 망각해버린 책을 뒤져 보다가 if 문은 저렇게 쓴다는 것을 알았다. ㅡㅡ;;

2010/05/26 - [공부해 봅시다/R-Project] - Array data 살펴 보기 - 삽질편

Creative Commons License

이 저작물은 크리에이티브 커먼즈 코리아 저작자표시-비영리-변경금지 2.0 대한민국 라이선스에 따라 이용하실 수 있습니다.

  Array 값을 살펴 보다가 의미 있는 차이를 보이는 그룹을 찾는 것도 좋을 것 같아서 방법을 시도해 보았다. 문법을
잘 몰라 한참을 삽질하다가 우여곡절 끝에 대충대충 해결할 수 있었다.

data <- read.table("Book1-1-3.csv", header=T, sep=",")

a <- matrix(nrow=5467, ncol=2)

for (X in 1:5467) {
x <- c(data[X,2],data[X,3],data[X,4],data[X,5], data[X,6],data[X,7],data[X,8],data[X,9])
y <- c(data[X,10],data[X,11])                       
a[X,1] <- X
a[X,2] <- t.test(x,y)$p.value
}

z <- which(a[,2] < 5e-02)
a[z,]

  삽질을 시작하기 전에 우선 5467x2 의 크기를 같은 a 라는 matrix 를 만들어 둔다.
a <- matrix(nrow=5467, ncol=2)

  첫 row 와 column 에는 문자값이 있기 때문에 그것을 제외하면 5467x10 의 값을 가지게 된다. 우선 1~8 column 의 값과 9, 10 column 의 값이 서로 유의한 차이가 있는지 알아보고 싶었다. 그리고 이렇게 구한 t-test 를 5467 번 반복해야 하는 문제도 있었다.
 
  우선 X 를 1 부터 5467 까지 순환하도록 했다.
for (X in 1:5467) { --- }

  그러한 X 값에 따라서 1~8 column 에 포함된 값을 x 라는 항목에 입력을 하고 9, 10 column 값을 y 라는 항목에 입력했다.
x <- c(data[X,2],data[X,3],data[X,4],data[X,5], data[X,6],data[X,7],data[X,8],data[X,9])
y <- c(data[X,10],data[X,11])        

  이러한 x 와 y 를 t-test 를 시행하였을 때의 p.value 만을 알고 싶었고, 이 값과 이 것이 몇 번째 X 값인지를 알아야 할 필요가 있었다.
a[X,1] <- X
a[X,2] <- t.test(x,y)$p.value

  여기까지 한 번에 구해지기 때문에 나중에 a 라고 입력된 값을 보면 다음과 같다.


  이 a 라는 matrix 에서 2 column 에 포함된 값중 p-value 가 0.05 이하인 값을 찾고 싶었기 때문에 다음과 같이 했다. which 를 사용해서 구하기는 하는데 매번 row 값을 출력하는 문제(??)가 있엇다.
z <- which(a[,2] < 5e-02)

  그래서 그냥 이 문제를 안고 살기로 했다. p.value가 0.05 이하인 값을 가지는 row 가 z 로 지정하도록 하고 이 것을 그냥 사용하는게 내가 알고 있는 지식으로는 최고의 결론이었다. ㅡㅡ;;
a[z,]

  유의 수준을 10e-07 까지 올리면 다음과 같이 나온다.

Creative Commons License

이 저작물은 크리에이티브 커먼즈 코리아 저작자표시-비영리-변경금지 2.0 대한민국 라이선스에 따라 이용하실 수 있습니다.

  Heatmap 은 역시 녹색, 검정, 빨강으로 표시가 되어야 제맛이다. :)
  R-Project에서 제공하는 Heatmap 으로는 안되고 별도의 패키지를 설치하여야 하며, 이 경우에는 gplots 패키지를 설치해야 한다. 그러면 heatmap.2 라는 명령어를 사용할 수 있으며, 이 명령어를 통해서 쉽게 만들 수 있다.

library(gplots)
heatmap.2(abc, col=redgreen(75), trace="none", density.info="none", hclustfun=my.hclust, distfun=my.dist)

아직 trace, density.info 에 대한 것은 확인하지 못했다. 없어도 늘 보는 것들과 차이가 없다는 점에서 필요 없을 것 같기는 하다. :)

Creative Commons License

이 저작물은 크리에이티브 커먼즈 코리아 저작자표시-비영리-변경금지 2.0 대한민국 라이선스에 따라 이용하실 수 있습니다.

  작년쯤에 겨우 했는데 기록을 해두지 않았더니 다시 기억을 끄집어 내는데 2시간쯤 걸렸다. 그래서 까먹지 않도록 여기에 기록을 해두어야 겠다는 생각이 들었다. 예전에 기록한 내용이라서 될 줄 알고 그냥 했더니 오류가 있어서 수정을 했다. 기본은 Data 파일을 불러와서, Matrix 로 변환한 다음 Heatmap 작성이다. 말로는 간단하고, 실제로도 간단하다. :)

arraydata <- read.table("b08.csv", header=T, sep=",")
설명) b08.csv 에 저장한 파일을 불러와서 arraydata 에 입력한다. 보통 자료를 관리할 때에 첫 열에 기본 항목을 입력하는 경향이 있는데 이러한 양식을 그대로 사용할 경우 Heatmap 을 만들 수가 없다. 물론, 첫 열을 인식시키지 않는 방법이 있을 것 같기는 하지만, 아직 방법을 찾지 못했다. 첫 줄의 항목은 Header 항목으로 인식하고, 자료 사이의 값은 ",(쉼표)" 에 의하여 분리되어 있다는 의미임. 엑셀에서 파일을 CSV 로 저장할 때에 'CSV (쉼표로 분리)' 라는 항목이 있는 것을 알고 있어야 함.

matrixarraydata <- as.matrix(arraydata)
설명) Distance 를 계산하기 위해서는 자료를 matrix 로 불러와야 한다고 함. 기본적인 개념을 알고 있어야 할 것 같기는 한데.. 아무튼 'as.matrix' 로 하면 한 방에 해결됨. 이러면 b08.csv 에서 불러온 파일은 matrixarraydata 라는 matrix 로 입력된다.

heatmap(matrixarraydata, hclustfun=my.hclust, distfun=my.dist)
설명) 결과값을 heatmap 으로 작성해 줌. Clustering 방법은 my.hclust 에서 지정된 방법, distance 는 my.dist 에서 지정된 방법으로 구함. 기본적으로 euclidean distance 에 complete method 로 clustering 이 되어 있어서 뒷 부분은 빼줘도 되지만, 다른 방법으로 해볼 것이라면 넣어 주는게 좋다.

my.hclust <- function(x){hclust(x, method="complete")}
설명) 대부분의 논문에서 생략하는 부분인데 이런 data 는 Chip 에 DNA 나 RNA 등을 반응시키는 방법부터 Scan 은 어떻게 하였으며, 어떻게 Normalization 을 시행 했는데 기술을 해줘야 한다. Normalization 은 아직 책에서 몇 번 보고 넘어간 부분이라서 생략하도록 하고, 외부에 의뢰를 하면 Normalization 까지는 되어져서 오기 때문에 구할 필요까지는 없다. Complete method를 사용할 경우에는 저렇게 표현하고 Ward method 일 경우에는 "ward"라고 하면 됨.

my.dist <- function(x){dist(x, method="euclidean")}
설명) 위에걸 보면 이것도 쉽게 알 수 있다.



 약 5000개의 값이 10개 있는 data로 만든 표이다. 처음에 중복된 값이 3개씩 있다는 것을 몰라서 15000개 정도의 값을 10개 있는 걸로 했더니 메모리 부족이라며 안되었는데.. 그래도 5000개 정도이면 시간이 좀 걸려서 그렇지 나오기는 한다. 시간을 측정하지는 않았지만 5~10분 정도 소요된다.

2009/06/21 - [공부해 봅시다/R-Project] - Heatmap
2009/06/23 - [공부해 봅시다/R-Project] - 함수 지정하기..ㅡㅡ

2008/05/03 - [지름신의 영접] - Do more in less time

Creative Commons License

이 저작물은 크리에이티브 커먼즈 코리아 저작자표시-비영리-변경금지 2.0 대한민국 라이선스에 따라 이용하실 수 있습니다.

  지난번 Heatmap 결과(2010/05/01 - [공부해 봅시다/R-Project] - Heatmap 만들기) 를 통하여 앞의 2개의 증례 가 다른 증례와는 다소 차이가 나는 것을 확인해 보았다. 하지만 Heatmap 을 만드는 과정에서 자료의 위치가 이동하기 때문에 원자료를 찾아보면 그 2개의 증례는 10번, 11번째 column 에 위치하는 것을 확인하였다.

  여기서 부터가 삽질이었다. 특정 질환에서 발현이 증가하고 다른 질환에서 발현이 감소하는 그러한 값을 찾는것이 기본이 되기 때문에 이 미묘한 자료를 찾는 것은 엑셀에서 순차 정렬을 11번 반복하게 된다.  OTL
  게다가 cut-off 값을 바꾸면 바꿀 때 마다 이 짓을 계속 해야하기 때문에 인터넷 바다를 검색한 끝에 해결책을 찾았다. 비록 결과물이 좀 보기 흉하기는 하지만, 내가 할 삽질을 컴퓨터가 대신 하도록 하는 방법을 찾을 수는 있었다.

data <- read.table("b08.csv", header=T, sep=",")
설명) 자료를 불러오는 것. 이 자료는 지난 자료와는 달리 첫 column 에 chip의 설명을 담고 있다. 특정 유전자를 찾기 위해서는 필요한 정보이다. 문자값이 있기 때문에 Heatmap 으로 진행은 할 수가 없다.

x <- which(data[,2] < X & data[,3] < X & data[,4] < X & data[,5] < X & data[,6] < X & data[,7] < X & data[,8] < X & data[,9] < X & data[,10] < X & data[,11] >= X & data[,12] >= X)
설명) 식이 좀 지저분하지만 Ctrl+C 와 Ctrl+V를 반복하면 쉽게 만들 수 있다. data 의 2~10 column 까지는 X 값보다 작은 것을 선택하고, 11, 12 column 은 X 값보다 같거나 큰 값을 선택한다. 아직 정확한 사용법을 알고 있지는 않지만 이렇게 한 후에 x 값을 살펴 보면 row 번호만 나왔다. 그래서 다음과 같은 식을 입력해주면 그 row 에 해당하는 값을 볼 수가 있다.

data[x,]

X 값을 1부터 대략 20 정도까지 일일이 입력하고 그 결과를 확인하면 되지만.. 그래도 반복되는 구문이라서 해결책을 찾을 수 있었다.

for (X in 1:20) {
x <- which(data[,2] < X & data[,3] < X & data[,4] < X & data[,5] < X & data[,6] < X & data[,7] < X & data[,8] < X & data[,9] < X & data[,10] < X & data[,11] >= X & data[,12] >= X)
print(X);print(data[x,])}
설명) 역시 R-project의 문법을 정확하게 이해하고 있지 않아서 생긴 지저분한 문장이다. {} 괄호로는 하나의 식을 입력할 수 있고, 엔터키로 문장을 바꾸면 실행을 하게 되는 것 같기는 한데.. ㅡㅡ;;
X가 1에서 20까지 1씩 증가하는 값으로 되어 있고 그에 따라서 x 값이 나온다. 결과를 보여주기 위해서 어떤 X 값인지 보여주는 print(X) 와 그 X 값에 따른 결과물을 볼 수가 있다.

Creative Commons License

이 저작물은 크리에이티브 커먼즈 코리아 저작자표시-비영리-변경금지 2.0 대한민국 라이선스에 따라 이용하실 수 있습니다.

  SPSS로 McNemar Test를 하려고 하는데 힘들었다. Matrix로 실행하는 방법이 없는 것 같은데 못찾은 건지는 모르겠다. R이 새삼스럽게 고마워지는 순간이었다. R의 장점 중 하나는 p=0.000 이라는 표현 대신에 지수를 사용해서 작은 값을 표현해준다. :)
  Data가 없이 표로만 자료가 주어지는 상황이라면 다음과 같은 방법으로 구할 수가 있다. Data 입력 순서와 보여지는 방식을 잘 보면 응용하기는 쉽다.

> drug <- matrix(c(55,9,43,77), nrow=2, dimnames=list("Before program"=c("No", "Yes"), "After program"=c("No", "Yes")));drug
              After program
Before program No Yes
           No  55  43
           Yes  9  77
> mcnemar.test(drug)

        McNemar's Chi-squared test with continuity correction

data:  drug
McNemar's chi-squared = 20.9423, df = 1, p-value = 4.733e-06

Creative Commons License

이 저작물은 크리에이티브 커먼즈 코리아 저작자표시-비영리-변경금지 2.0 대한민국 라이선스에 따라 이용하실 수 있습니다.

 Odds ratio를 R을 사용하여 구하기 위해서 정말 노력해 보았다. 고민하다가 Appleforum 게시판에서 R을 사용한다는 분을 찾아서 문의 메일을 보냈고 그 답장이 왔다. 그대로 해보니 그 동안 구글링한 내가 바보 같이 느껴질 정도로 쉽게 구해졌다.  OTL

  우선 Package 중에서 epitools를 설치한다. 그 후에

library(epitools)
oddsratio(htn.rawdata[,9], htn.rawdata[,1], conf.level=0.95) 이렇게 구하면 된다. 하지만.. 일반적으로 책에서는 좌측 상단에 우리의 목표값(??)이 위치하게 되는데 도움말을 참고해보면 우리의 함수님께서는 우측 하단에 목표값(??)이 위치해야 하므로 columns을 바꾸어 주어야 한다.

> library(epitools)
> oddsratio(htn.rawdata[,9], htn.rawdata[,1], conf.level=0.95, rev="columns")
$data
         Outcome
Predictor   2   1 Total
    0     601 819  1420
    1      55 114   169
    Total 656 933  1589

$measure
         odds ratio with 95% C.I.
Predictor estimate    lower    upper
        0 1.000000       NA       NA
        1 1.518482 1.087011 2.144556

$p.value
         two-sided
Predictor midp.exact fisher.exact chi.square
        0         NA           NA         NA
        1 0.01401735   0.01637873 0.01464472

$correction
[1] FALSE

attr(,"method")
[1] "median-unbiased estimate & mid-p exact CI"

Creative Commons License

이 저작물은 크리에이티브 커먼즈 코리아 저작자표시-비영리-변경금지 2.0 대한민국 라이선스에 따라 이용하실 수 있습니다.

이번 숙제에 R을 사용해 보기 위해서 계속 공부중이다. Odds ratio 를 구하는 함수는 없는 것 같지만 계속 찾아볼 예정이고 평균치 비교는 비교적 쉽게 찾을 수 있었다. 이건 비통계전공자의 전투족에 해당하는 내용이다. :)

기본적으로 함수는 t.test 이다. 간단하다. 그래서 어렵다. 아파트나 자동차와 마찬가지로 옵션을 넣어주어야 한다.
사용하는 기본 옵션은 다음과 같다.
var.equal = T(RUE) or F(ALSE); 두 집단이 등분산이면 TRUE를 아니면 FALSE를 선택한다.

Creative Commons License

이 저작물은 크리에이티브 커먼즈 코리아 저작자표시-비영리-변경금지 2.0 대한민국 라이선스에 따라 이용하실 수 있습니다.