Category: R

주의: 통계 비전공자의 글이므로 용어의 대한 설명이 부적절할 수 있음.

Bootstrapping.

여러 곳에서 사용되는 단어이자 개념인듯하다.

역시 통계에서도 사용된다.

인터넷을 통한 검색에 의하면

표본크기의 한계를 극복하고 모수통계기법의 비현실성을 극복하기 위한

대안으로 널리 사용되는 방법이라고 한다.

나의 실험에 사용한 하나의 실험 자료가 있고, 이 자료에는 n개의 샘플이 있다고 하자.

그 동안의 통계는 이 실험자료는 모집단에서 추출한 자료로 간주를 하고,

그에 따른 통계적인 분석을 하는 것이 일반적이다.

Bootstrapping 방법에서는 나의 실험자료를 하나의 모집단이라고 가정한다.

이 자료에서 임의로 자료를 추출해서 다른 자료를 만들어가는데,

추출 과정에서 중복 추출이 가능하다(복구 replacement 라고도 하는 듯).

이렇게 추출을 하는데 원래 모집단과 같은 숫자인 n개의 자료를 추출한다.

이렇게 만든 표본을 bootstrapping sample 이라고 한다.

중복 추출이 가능하기 때문에

원자료와 똑같은 표본을 만들 수도 있고,

빈도는 낮을지라도 하나의 자료로만 구성된 sample 도 만들어 질 수 있다.

이렇게 표본을 만들어서 그 결과를 바탕으로 통계 분석을 시행하고,

그 분석을 바탕으로 하여 원자료에서는 어떻게 나올지에 대하여 추론하는 것이다.

표본을 만드는 방식에는 가능한 모든 경우를 다 시행하는 경우 (exact version .?)

그리고 임의로 만드는 방법이 있다.

통상적으로 사용하는 방식은 임의적으로 만들어내는 방식이다.

이렇게 반복되는 계산을 바탕으로 하여 원래의 값을 추정해 나가는 방식의

논리적인 전개 방법을 Monte Carlo Algorithm 이라고 하는 듯 함.

---

따라서, 고려해야 하는 부분은 다음과 같다.

1. 임의로 추출하는 과정이 정말로 임의적(random)일 것이냐 하는 것

2. 그러면 몇 번이나 반복을 해야하는가(k, B 등으로 표현됨)

임의성(Randomness)에 관한 것은 

설명을 듣는다고 해서 내가 이해할 수 있는 것이 아니다.

임의의 수는 R에서 기본적으로 Mersenne-Twister 방식에 의해서 만들어진다고 한다.

이 MT 방식은 암호학적으로 안전한 방식은 아니지만,

통상적인 통계 방법에 적용하기에는 충분하다고 한다.

k의 범위에 대해서는 실험의 성격에 따라서 달라진다고 하는데,

통상적인 실험의 경우라면 1,000번이면 충분하다고 한다.

금융쪽에서는 100만번 정도 반복한다고 한다.

물론 요즘은 컴퓨터 계산속도가 빠르고,

일반적인 의학 연구의 샘플의 크기가 큰 편이 아니기 때문에

k를 좀 늘려도 되기는 하지만, 

통계적으로 얼마나 의미가 있는지는 모른다. 

나의 샘플에서는 100번, 1,000번, 10,000번을 시행하였을 경우,

분포 곡선은 크게 차이 나지 않는 것을 확인하였다.

약간 차이가 있다면 반복횟수가 증가하면 분포 곡선이 더 미끈해지는 것 같다고나 할까..

따라서, 반복 횟수를 10,000번 이상 늘릴 필요는 없으며,

다른 아이디어를 적용하는 것이 시간을 더 효율적으로 사용할 수 있을 것으로 생각된다.

이러한 기본 컨셉을 확인한 후에 bootstrapping 방법을 응용하여

Kaplan-Meier survival analysis 에서 p-value 의 분포를 확인해 보기로 하였다.

R에서 기본으로 제공하는 함수 중 

sample 이라는 함수가 random resampling 을 해주며, 

replace=TRUE 옵션을 통해서 replace 여부를 설정해 줄 수 있다.

이 함수를 이용해서 원래 자료 중에서 몇 번째 자료를

추출하여 bootstrapping sample 을 만들지 정할 수 있다.

survdiff 함수를 통한 생존의 차이를 평가하는 함수는

p.value 를 직접 추출할 수 없기 때문에 간접적으로 추출하여 

직접 계산을 해주어야 한다고 한다.

인터넷을 통하여 방법을 확인할 수 있었다.

그 동안의 R 삽질로 배운 것이 있다면,

반복 구문에서 사용되거나 계산하는 값들은 한 번 실행이 끝날 때마다 

계속 지우고 다시 만들어야

문제가 생겼을 때 어디쯤에서 잘못 했는지 확인하기가 쉽다.

bstr.k <- 10000; 몇 번 반복할지를 미리 지정해주고..

bstr.matrix <- tmp4[,c(tmp01, tmp02, tmp03)]; survdiff 함수에서 필요한 3개의 값들만을 모은 후에

bstr.n <- nrow(bstr.matrix)

tmp1 <- matrix(nrow=bstr.k, ncol=1); p.value 의 분포를 확인하기 위하여 미리 값을 저장할 곳을 만들고

for (Y in 1:bstr.k){

  bstr.order <- sample(c(1:bstr.n), replace=TRUE)

  bstr.temp.matrix <- matrix(nrow=bstr.n, ncol=ncol(bstr.matrix)); bootstrapping sample 을 위한 빈 matrix 만들고

  for (X in 1:bstr.n){bstr.temp.matrix[X,] <- bstr.matrix[bstr.order[X],]}

  tmp2 <- survdiff(Surv(bstr.temp.matrix[,1], bstr.temp.matrix[,2], type=”right”) ~ bstr.temp.matrix[,3])

  p.val <- 1 – pchisq(tmp2$chisq, length(tmp2$n) – 1)

  tmp1[Y, 1] <- p.val

  rm(bstr.order, bstr.temp.matrix, tmp2, p.val)

}

hist(tmp1[,1], breaks=20, freq=TRUE); Histogram

length(which(tmp1[,1] < 0.05))/bstr.k; 0.05 미만인 것들의 빈도

10,000번을 반복했을 경우, p-value 가 0.05 미만일 가능성은 30% 정도 된다.(붉은 색)

원자료만을 바탕으로 생존 분석을 시행하였을 때에는

통계적으로 의미있는 차이를 보여주지는 않았다.

하지만 bootstrapping 을 통하여 분석을 시행한 결과를 언급하면 

조금은 더 그럴듯한 결론을 만들어 낼 수 있다.

즉, 나의 실험 결과를 바탕으로 모집단에서의 빈도를 추정하여 보았을 때

30% 정도는 통계적으로 유의한 결과를 보여준다.

뭔가 그럴듯한 결론을 이끌어 낼 수 있다는 점에서

이 통계기법은 조금은 더 깊게 알아볼만한 가치가 있는 것 같다.

SPSS로 생존자료를 확인해보다가 R이라면 다른 방법이 있을 줄 알고

시도해보기로 했다.

우선 CSV 로 파일을 저장한다음 R에서 불러왔다.

제대로 불러왔는지 확인해보기 위하여

생존 곡선을 그려보기로 하였다.

다행히 SPSS와 같은 모습의 생존 곡선이 나왔다.

이 중에서 내가 관심있는 부분은 파란선과 빨간선으로 표시된 그룹이다.

얼핏 보기엔 유의한 생존곡선의 차이를 보이는 것 같아서 확인을 해보았다.

그런데 실제 결과는 아닌 것으로 나왔다.

아마도 통계적으로 유의한 차이를 보이기에는

증례수가 부족하여 검정력이 안생기는 것으로 보였다. 

원자료를 다시 샘플링 하는 방법을 기반으로 하는 Bootstrapping 을 이용하면

뭔가 나올 줄 알고 

먼가 해볼려고 했는데

censboot 나 bootkm 같은걸 이용해볼려고 했는데

도저히 이해할 수 없는 영역이라는 것을 알았다.

그래서 통계적인 방법이라고는 할 수 없지만

기존 증례를 2배로 불려서 검증을 해보았다.

통계적 유의성이 보여진다.

그래서 TMA 를 한 set 더 염색을 해서 생존 자료 분석을

해보기로 했다. 

뭔가 좋은 결과가 나왔으면 좋겠다.

---

이리저리 찾아보던 중 발견한 것은

1. SPSS에서와 마찬가지로 median 값이 안나오는 이유

인터넷을 찾아본 결과

절반이 죽지 않아서 안나온다.

2. mean survival time 구하는 방법

median survival time 과 비슷한 이유이기는 한데

모든 증례에서 생존값이 나와있지 않기 때문에

mean 값도 구할 수 없는 값이다.

SPSS 에서 KM 분석을 하면 늘 mean survival time 에서 

뭔가 주석을 달아서 표시하는 내용이 있는데

그게 그런 이유였다. ㅡ_ㅡ;;

사진 촬영할 때 마다 미묘하게 어긋날 수 있다는 문제를 해결하려면..

결국 어떤 기준이 되는 점을 정해서 그 점을 기준으로 회전을 하던가 해야할 것이다.

Optic disc의 중간 쯤 되는 어떤 부위가 가장 좋을 것 같아서

이 점을 찾는 어떤 알고리즘을 구현하기로 하였다.

RGB 3개의 색중에서는 Red가 그나마 좀 Optic disc를 표현하는데 적절하지 않을까 생각된다.

Cut off value 는 기억이 안나는데 아무튼 꽤 높은 cut off value 를 선택해서

그 이상이 되는 빨간색만 남게 하면 대충 이런 식이 되게 된다.

우선 optic disc는 기본적으로 타원형의 모양을 하고 있는 것 같다.

안과 전문의의 옥음에 의하면 모양도 다양하다고는 하지만..

그런거 일일이 고려하면 끝도 없으니 우선 타원형을 기준으로 하였다.

Optic disc 이외의 부분에서 굉장히 R값이 높은 부분이 관찰된다.

하지만, 이런 부분은 그 면적이 작기 때문에

해당 부분은 고려는 하되, 중요하지는 않도록 처리를 해야한다.

원점을 기준으로 하는 타원의 공식은 다음과 같다.

타원 방정식에서 a, b 값에 근사할 만한 수치를 지정해준다.

타원이 그림 파일을 쭈~욱 훑으면서 움직이게 하면서

타원 내부에 진한 빨간점이 몇 개인지 혹은 어느 정도인지 평가한다.

이 방법이라면

내가 원하는 기준점쯤에서는 타원 내부에 들어가 있는 빨간점들이 가장 많을 것이다.

나머지는 검색 시간을 줄여주기 위하여 훑어나가는 범위를 줄여주기만 하면 된다.

내 머리속의 평면 좌표는 (X, Y)로 구성되어 있는데

Matrix에서는 [Y, X] 로 되어 있어서

한 참 헷갈린 끝에 적당한 점을 찾는데 성공했다.

이 방법의 문제가 있다면…

가장 큰 문제인데..

거의 흑백에 가깝게 나오는 Red free filter 사진에 적용시키는게 문제라고..

OTL…

자.. 다시 시작하자. 😉