## ## t-test 이해 확인 ## remember t test = diff / rand error ## set.seed(101) rnorm2 <- function(n,mean,sd){ mean+sd*scale(rnorm(n)) } A <- rnorm2(16, 26, sqrt(1160/15)) B <- rnorm2(16, 19, sqrt(1000/15)) A <- c(A) B <- c(B) # we know sqrt(1160/15) is A's sdev # hence, A's var is sqrt(1160/15)^2 # hence, A's SS is sqrt(1160/15)^2 * 15 # this is 1160 # from the above, # the difference between the A and B means # remember we try to find # difference due to the treatment / # / random chance of error diff <- mean(A) - mean(B) # for se # we know that the situation refers to # #2 se two independent samples t-test # which is sqrt(pooled.var/na + pooled.var/nb) SSa <- 1160 SSb <- 1000 n.a <- 16 n.b <- 16 df.a <- n.a - 1 df.b <- n.b - 1 # 확률과통계 수업 학생은 아래에서 # se값을 구한 결과만 보고 이것이 # se = sqrt(sigma^2/n)에 # 대응하는 값으로 생각하세요 # we know that we are testing the difference # between two independent sample means. # Hence, we need to use poole variance between # the two group. See # http://commres.net/wiki/t-test#t-test_%EB%B9%84%EA%B5%90 pooled.var <- (SSa + SSb) / (df.a + df.b) se <- sqrt( (pooled.var/n.a) + (pooled.var/n.b) ) # Remember t test calculation is based on # diff / random error t.calculated <- diff / se pooled.var diff se t.calculated # Now use t.test function for two group # (independent sample) t-test # with an assumption that variances of # the two gorup are the same. t.result <- t.test(A, B, var.equal = T) t.result # t.result$statistic = t.calculated # t.result$p.value = probability level of # wrong decision with the t calculated value str(t.result) t.result$statistic t.result$p.value # the above p.value can be obtained with # pt function p.value.calculated <- 2*pt(t.calculated, df = df.a + df.b, lower.tail=FALSE) p.value.calculated t.result$p.value ## ## different approach (분산으로 분석) ## # A combined group with group A and B # We call it group total # we can obtain its mean, variance, ss, df, etc. # A B dat <- c(A,B) dat mean.total <- mean(dat) var.total <- var(dat) # variance를 ms라고 부르기도 한다 # ms = mean square (square 값을 n 으로 나눈 값이 분산이므로) ms.total <- var.total df.total <- length(dat) - 1 ss.total <- var.total * df.total # ss.total 은 아래처럼 구해도 된다 ss.total.check <- sum((dat-mean(dat))^2) mean.total var.total ms.total df.total ss.total ss.total.check # Now for each group mean.a <- mean(A) mean.b <- mean(B) mean.a mean.b # 그룹 간의 차이에서 나타나는 분산 # 수업시간에 설명을 잘 들을 것 # mean.total 에서 그룹a의 평균까지의 차이를 구한 후 # 이를 제곱하여 그룹 A 멤버의 숫자만큼 곱한다 # 이를 그룹b의 평균과도 다시 계산하고 # 이 둘을 더하여 # ss.between에 저장 length(A) * ((mean.total - mean.a)^2) length(B) * ((mean.total - mean.b)^2) ss.between <- length(A)*((mean.total - mean.a)^2) + length(B)*((mean.total - mean.b)^2) ss.between # df between group은 연구에 사용된 # 그룹의 숫자에서 1을 뺀 숫자 df.between <- 2 - 1 # 이 그룹 간 차이에 기인하는 분산 값은 ms.between <- ss.between / df.between # 한편 아래 ss.a 와 ss.b는 각 그룹 내의 # 분산을 알아보기 위한 방법 ss.a <- var(A) * df.a ss.b <- var(B) * df.b ss.within <- ss.a + ss.b df.a <- length(A)-1 df.b <- length(B)-1 df.within <- df.a + df.b ms.within <- ss.within / df.within # 여기까지 우리는 # 전체분산 # 그룹간분산 # 그룹내분산 값을 # 구한 것 # ms.between은 그룹의 차이때문에 생긴 # 분산으로 IV 혹은 treatment 때문에 생기는 # 차이에 기인하는 분산이다. # ms.within은 각 그룹 내부에서 일어나는 분산이므로 # (variation이므로) 연구자의 관심사와는 상관이 없이 # 나타나는 random한 차이를 알려주는 분산이라고 하면 # t test 때와 마찬가지로 # 그룹의 차이 / 랜덤 차이를 (에러 -> 분산은 에러라고도 했다) # 구해볼 수 있다. # 즉, 그룹갑분산은 사실 = diff (between groups) # 그리고 그룹내 분산은 사실 = re # 따라서 우리는 위 둘 간의 비율을 t test와 같이 # 살펴볼 수 있다 # 이것을 f.calculated 이라고 하고 f.calculated <- ms.between / ms.within # 이 값을 출력해 본다 f.calculated # 이 계산은 차이와 랜덤에러의 비율이 # df에 따라서 얼마나 되어야 그 차이가 # 충분히 큰 것인지를 판단하기 위해서 # 쓰인다. 여기서 df에는 두 가지 종류가 # 있다. df.between 그리고 df.within # 이 정보를 이용하여 f.calculated 값에 대한 # probability를 구할 수 있다 (t-test때 처럼) # percentage of f distribution with # df1 and df2 option # 이는 그림의 왼쪽을 나타내므로 # 차이가 점점 커지게 되는 오른쪽을 # 계산하기 위해서는 1-x를 취한다 f.calculated.pvalue <- 1-pf(f.calculated, df1=df.between, df2=df.within) f.calculated.pvalue # 한편, t test를 했었을 때 (A, B 그룹을 가지고 independent # samples t-test를) 아웃 풋은 t.result # 그리고 f 계산에서의 p value는 t test에서의 p.value와 같다 f.calculated.pvalue t.result$p.value # 또한 # 여기엣 t 값은 t.result$statistic 으로 프린트아웃할 수 있다 # 이 값이 2.33333 이었다 t.result$statistic # 혹은 우리가 계산한 값이었던 # t.calculated t.calculated # 그런데 위의 값을 제곱한 값이 바로 f.calculated 값 f.calculated t.calculated t.calculated^2 # 혹은 f.calculated 값을 제곱근한 값이 t.calculated f.calculated sqrt(f.calculated) t.calculated # 또한 우리는 계산한 ms.total의 분자 부분인 ss.total은 # 독립변인이 영향을 주는 ss.between 값과 # 독립변인이 있음에도 불구하고 여전히 남는 ss.within # 으로 이루어짐을 아래처럼 확인한다. ss.total ss.between ss.within ss.total ss.between + ss.within # 한편 df는 # df.total 30 - 1 df.total df.between df.within df.total df.between + df.within # 한 편 # 아래는 위의 분산을 쪼개서 분석하는 펑션인 # aov를 이용해서 통계분석을 한 결과이다. A B comb <- stack(list(a=A, b=B)) comb colnames(comb)[1] <- "values" colnames(comb)[2] <- "group" comb a.res <- aov(values ~ group, data=comb) a.res.sum <- summary(a.res) a.res.sum # 위에서 F value는 5.444 # 그리고 전체적인 아웃풋을 보면 # Df group 과 Df Residuals # Sum Sq group 과 Residuals # Mean Sq (MS) group 과 MS residuals # MS.group = ms.between # MS.within = ms.within # F value ms.between / ms.within f.calculated # 아래는 기존의 아웃풋에서 확인하는 것 str(a.res.sum) a.res.sum[[1]][1,4] sqrt(a.res.sum[[1]][1,4]) t.result$statistic