Category: R

  다른 연구 자료의 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를 알 수가 있다.

  오늘 아침에 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(dat, col=rainbow(20*10, start = 0/6, end = 4/6), hclustfun=my.hclust, distfun=my.dist, dendrogram=”both”, trace=”none”, density.info = “none”, key = TRUE)