R에서 쌍을 이룬 t- 검정의 검정력을 수동으로 계산
연습으로 R에서 수동으로 paired t-test를 수행하여 과거에했던 강의를 새로 고치고 싶었습니다. 모든 것이 순조롭게 진행되었지만,이 쌍을 이룬 t- 검정의 검정력을 계산하는 것에 대해 생각해 보았고 여기서 문제가 시작되었습니다.
검정력은 대체 분포 아래 영역에서 유형 II 오류 영역을 뺀 영역이라는 것을 알고 있습니다 ($\beta$)로 구분됩니다. $\alpha$유의 수준. 따라서 기본적으로이 예제에서는$P(X ≤ \alpha)$내가 계산 한 관측 된 평균 차이를 중심으로하는 대체 분포의 경우, 솔직히 그 분포를 구성하는 방법을 잘 모르겠습니다. Null 아래에서 t- 통계에 대해 동일한 절차를 사용하려고했지만 예상 평균과 관측 된 평균이 동일하므로 전체 항이 0 ( 1-pt((expMean - obsMean)*stdError, df) 과 같기 때문에 의미가 없습니다 . 그리고 내가 아는 한 t- 분포는 귀무 가설이 참이라는 가정 하에서 만 사용됩니다. 여기서부터 나는 더 혼란스러워지고 분명한 것을 놓치고 있다고 생각합니다.
결과를 비교하기 위해 pwr 패키지 의 pwr.t.test 함수를 사용했습니다 .
다른 곳에서 찾은 대부분의 솔루션은 내가 수동으로 수행하려는 부분을 건너 뛰고 일종의 전력 계산기를 사용하기 때문에 누군가가 이러한 테스트를 수동으로 수행하도록 도와 줄 수 있다면 매우 유용 할 것입니다.
내가 사용한 코드 :
# data
aP <- c(0.5331039, 0.4578532, 0.3129205, 0.5144858, 0.8149759, 0.4136268)
aM <- c(0.2750040, 0.5056830, 0.4828734, 0.4439654, 0.2738658, 0.3081768)
# difference between P and M
Diff <- aM - aP
# INIT t test
obsMean <- mean(Diff)
expMean <- 0
stdError <- (sqrt(length(Diff))/sd(Diff))
n <- length(aP)
df <- n - 1
alpha = 0.05
# T-statistic
T_stat <- (obsMean-expMean)*stdError; T_stat
# critical value
crit_values <- qt(c(0.025,0.975),df) # lower bound = -2.570582
p_value <- 2*(pt(T_stat, df)); p_value
p_value < alpha
# comparison
t.test(aM, aP, paired = TRUE, alternative = "two.sided")
# INIT power
obsMean <- mean(Diff)
expMean <- mean(Diff)
# power???
power <- 1-pt((expMean - obsMean)*stdError, df); power
# comparison
cohensD <- (mean(aM)-mean(aP))/(sqrt((sd(aM)^2+sd(aP)^2)/2))
pwr.t.test(n = 6,d = cohensD, type = "paired", alternative = "two.sided")
# power = 0.4210006
```
답변
나는 여기서 속이고있다 ... 나는 pwr.t.test에 대한 코드를 찾고 쌍을 이루는 양면 t-test의 힘을 생성하기 위해 관련 부분을 추출했습니다.
귀하의 의견 :
aP <- c(0.5331039, 0.4578532, 0.3129205, 0.5144858, 0.8149759, 0.4136268)
aM <- c(0.2750040, 0.5056830, 0.4828734, 0.4439654, 0.2738658, 0.3081768)
cohensD <- (mean(aM)-mean(aP))/(sqrt((sd(aM)^2+sd(aP)^2)/2))
pwr.t.test(n = length(aP), d = cohensD, type = "paired", alternative = "two.sided", sig.level= 0.05)
# power = 0.4210006
수동으로 재현하려면 :
n <- length(aP)
tsample <- 1 # 1 because type is paired
tside <- 2
sig.level <- 0.05
d <- cohensD
nu <- (n - 1) * tsample
qu <- qt(sig.level/tside, nu, lower = FALSE)
pt(qu, nu, ncp = sqrt(n/tsample) * d, lower = FALSE) +
pt(-qu, nu, ncp = sqrt(n/tsample) * d, lower = TRUE)
# [1] 0.4210006
편집 다음은 위 코드의 주석이 달린 버전입니다.
주어진 유형 1 오류로 쌍을 이룬 t- 검정의 검정력을 계산하려고합니다. $\alpha = 0.05$및 효과 크기 (Cohen의 d)는 샘플 쌍에 의해 결정됩니다 aP, aM. 따라서 입력은 다음과 같습니다.
aP <- c(0.5331039, 0.4578532, 0.3129205, 0.5144858, 0.8149759, 0.4136268)
aM <- c(0.2750040, 0.5056830, 0.4828734, 0.4439654, 0.2738658, 0.3081768)
sig.level <- 0.05
cohensD <- (mean(aM)-mean(aP))/(sqrt((sd(aM)^2+sd(aP)^2)/2))
먼저 5 % 사례에서 귀무 가설을 잘못 받아들이는 t- 통계의 임계 값을 찾아야합니다. 테스트는 양면이므로 다음 값을 찾는 것을 의미합니다.$x$ 아래 그림의 확률 밀도 함수에서 두 개의 음영 처리 된 꼬리를 정의하며, 음영 처리 된 각 영역은 전체 영역의 2.5 %입니다.
이를 위해 분위수 함수를 다음 qt과 같이 사용할 수 있습니다.$n - 1$ 자유도:
df <- (length(aP) - 1)
qu <- qt(sig.level/2, df, lower = FALSE)
# Code for plot
x <- seq(-6, 6, length.out= 100)
y <- dt(x, df= df)
plot(x, y, type= 'l', lwd= 1.5, xlab= 'Value of T', ylab= 'Density')
polygon(c(x[x > qu], qu), c(y[x > qu], 0), col= "grey", border= 'black')
polygon(c(x[x < -qu], -qu), c(y[x < -qu], 0), col= "grey", border= 'black')
-Inf 사이 와 사이 와 Inf 사이에 PDF를 통합 하여 임계 값 qu(및 -qu)이 영역의 2.5 %를 정의 하는지 확인할 수 있습니다 .-ququ
integrate(dt, -Inf, -qu, df= df) # -> 0.025 with absolute error < 6.1e-05
integrate(dt, qu, Inf, df= df) # -> 0.025 with absolute error < 6.1e-05
이제 귀무 가설이 거짓이고 평균 간의 차이가 0이 아니라 원하는 Cohen의 d가 있다고 가정합니다. 그래서 우리는 효과 크기의 방향으로 치우치는 비 중심 모수를 가진 t- 분포를보고 있습니다. R 문서가 NCP를 설명하는 방법은 다음과 같습니다.
가장 많이 사용되는 응용 프로그램은 t- 검정에 대한 검정력 계산입니다. Let T = (mX-m0) / (S / sqrt (n)) 여기서 mX는 '평균'이고 S는 X_1의 표본 표준 편차 ( 'sd')입니다. X_2, ..., X_n iid N (mu, sigma ^ 2) 그러면 T는 'df'= n-1 자유도 및 n on- c entrality p arameter 'ncp'=를 갖는 비 중심 t로 분포됩니다. (mu-m0) * sqrt (n) / sigma.
그래서 우리는 :
ncp <- sqrt(length(aP)) * cohensD
우리는 중요한 값 밖에이 NCP와 자유의 정도와 t 분포의 비율 지역 알고 싶어 -qu하고 qu위의를. 즉, 아래에 음영 처리 된 영역이 필요합니다 (오른쪽 꼬리 영역은 거의 보이지 않음).
right <- pt(qu, df, ncp = ncp, lower = FALSE)
left <- pt(-qu, df, ncp = ncp, lower = TRUE)
right + left
[1] 0.42 # As per pwr.t.test()
# Code for plot
x <- seq(-12, 5, length.out= 200)
y <- dt(x, df= df, ncp= ncp)
plot(x, y, type= 'l', lwd= 1.5, xlab= 'Value of T', ylab= 'Density')
polygon(c(x[x > qu], qu), c(y[x > qu], 0), col= "grey", border= 'black')
polygon(c(x[x < -qu], -qu), c(y[x < -qu], 0), col= "grey", border= 'black')
abline(v= c(-qu, qu), lty= 'dashed', col= 'blue')
다시 PDF를 통합하여 확인할 수 있습니다.
integrate(dt, -Inf, -qu, df= df, ncp= ncp) # -> 0.42 with absolute error < 1.3e-05
integrate(dt, qu, Inf, df= df, ncp= ncp) # -> 6.9e-05 with absolute error < 2.8e-08
이것이 도움이되기를 바랍니다 (올바른지 확인하십시오)!
쌍체 t 검정은 차이에 대한 1- 표본 검정입니다. $D_i = X_i-Y_i,$ ...에 대한 $i=1,2, \dots, n$ 과 $D_i$ 독립적으로 $\mathsf{Norm}(\mu_D, \sigma_D).$
테스트 고려 $H_0:\mu=0$ 대 $H_a:\mu > 0$ 5 % 수준에서 $n = 25.$ 특정 대안에 대한 테스트의 힘을 추구합니다. $\mu = \mu_a = 2 > 0.$
힘을 찾으려면 다음의 가치에 대한 교육적인 추측이 필요합니다. $\sigma.$ 와 $\alpha = 0.05, n = 25, \sigma = 3,$ 찾을 수 있습니다 $P(\mathrm{Rej\;} H_0\,|\, \mu=\mu_a).$[물론, 경우 당신은 정확한 값을 알고 의를$\sigma,$ 그러면 t-test 대신 z-test를하게됩니다.]
Minitab 소프트웨어 : 다음은 Minitab 최신 릴리스의 관련 출력입니다. [R 및 기타 통계 소프트웨어 프로그램에는 유사한 절차가 있습니다. @dariober의 답변 (+1)은 양측 테스트에 대해 간단히 언급합니다.]
지정된 매개 변수의 검정력은 다음과 같습니다. $\pi = 0.944.$ [제 2 종 오류 확률은 $\beta = 1 - \pi = 0.065.]$
Power and Sample Size
1-Sample t Test
Testing mean = null (versus > null)
Calculating power for mean = null + difference
α = 0.05 Assumed standard deviation = 3
Sample
Difference Size Power
2 25 0.944343
시뮬레이션. 100,000 회 반복으로 약 2 자리 정확도를 예상 할 수 있습니다. R에서 다음 시뮬레이션의 대략적인 결과는 다음과 같습니다.$\pi = 0.945.$
set.seed(2020)
pv = replicate(10^5, t.test(
rnorm(25, 2, 3), alt="g")$p.val)
mean(pv <= 0.05)
[1] 0.9449
비 중심 t 분포 사용.
A에 대한 임계 값 (하나는 양면)의 테스트 $ H_0은 : 0 $ = \ MU 대 \ 민> 0 $ : $ H_a가 갖는 5 % 수준에서 $ N = 25 $가 있다 . C = 1.7109 $ $ 즉 , $ T_0 = \ frac {\ bar D-0} {S_D. \ sqrt {n}} \ ge c. $ 이면 $ H_0 $ 을 거부 합니다.
c = qt(.95, 24); c
[1] 1.710882
우리는 추구 = 0.9443, $ - $ P \ 왼쪽 (오른쪽 {S_D / \ SQRT {N}} \ GE의 C \ T_a = \ FRAC {\ mu_a \ 바 D}) 여기서 $ T_a $가 가 noncentral t 분배 의 학위를 자유 $ \ nu = n-1 = 24 $ 및 비중 심성 매개 변수 $ \ delta = \ sqrt {n} (2) / 3 = 10 / 3. $ [R CDF 함수의 세 번째 매개 변수 df가 비중 심성임을 유의하십시오. 매개 변수.]
del = 5(2)/3
1 - pt(c, 24, del)
[1] 0.9443429