Мощность для двух образцов т тест

Я пытаюсь понять расчет мощности для случая двух независимых выборочных t-тестов (не предполагая равных отклонений, поэтому я использовал Satterthwaite).

Вот диаграмма, которую я нашел, чтобы помочь понять процесс:

Итак, я предположил, что, учитывая следующее о двух популяциях и размер выборки:

mu1<-5
mu2<-6
sd1<-3
sd2<-2
n1<-20
n2<-20

Я мог бы вычислить критическое значение под нулем, относящееся к 0,05 вероятности верхнего хвоста:

df<-(((sd1^2/n1)+(sd2^2/n2)^2)^2) / ( ((sd1^2/n1)^2)/(n1-1) + ((sd2^2/n2)^2)/(n2-1)  )
CV<- qt(0.95,df) #equals 1.730018

и затем вычислите альтернативную гипотезу (которая для этого случая, которую я узнал, является «нецентральным t-распределением»). Я рассчитал бета на диаграмме выше, используя нецентральное распределение и критическое значение, найденное выше. Вот полный скрипт на R:

#under alternative
mu1<-5
mu2<-6
sd1<-3
sd2<-2
n1<-20
n2<-20


#Under null
Sp<-sqrt(((n1-1)*sd1^2+(n2-1)*sd2^2)/(n1+n2-2))
df<-(((sd1^2/n1)+(sd2^2/n2)^2)^2) / ( ((sd1^2/n1)^2)/(n1-1) + ((sd2^2/n2)^2)/(n2-1)  )
CV<- qt(0.95,df)


#under alternative
diff<-mu1-mu2
t<-(diff)/sqrt((sd1^2/n1)+ (sd2^2/n2))
ncp<-(diff/sqrt((sd1^2/n1)+(sd2^2/n2)))


#power
1-pt(t, df, ncp)

Это дает значение мощности 0,4935132.

Это правильный подход? Я обнаружил, что если я использую другое программное обеспечение для расчета мощности (например, SAS, которое, как мне кажется, я настроил в соответствии со своей задачей ниже), я получу другой ответ (для SAS это 0,33).

SAS-код:

proc power;
      twosamplemeans test=diff_satt
         meandiff = 1
         groupstddevs = 3 | 2
         groupweights = (1 1)
         ntotal = 40
         power = .
        sides=1;
   run;

В конечном счете, я хотел бы получить понимание, которое позволило бы мне взглянуть на симуляции для более сложных процедур.

РЕДАКТИРОВАТЬ: Я нашел свою ошибку. должны были быть

1-pt (CV, df, ncp) НЕ 1-pt (t, df, ncp)

Вы близки, хотя некоторые небольшие изменения необходимы:

• $\mu_{2} - \mu_{1}$
• $n_{1}+n_{2} - 2$$t$
• SAS может использовать формулу Уэлча или формулу Саттервейта для df с учетом неравных отклонений (вы найдете в этом PDF-файле, который вы цитировали ) - только с 2 значащими цифрами в результате, который невозможно определить (см. Ниже)

С n1, n2, mu1, mu2, sd1, sd2такими , как определено в вашем вопросе:

> alpha   <- 0.05
> dfGP    <- n1+n2 - 2                     # degrees of freedom (used by G*Power)
> cvGP    <- qt(1-alpha, dfGP)             # crit. value for one-sided test (under the null)
> muDiff  <- mu2-mu1                       # true difference in means
> sigDiff <- sqrt((sd1^2/n1) + (sd2^2/n2)) # true SD for difference in empirical means
> ncp     <- muDiff / sigDiff              # noncentrality parameter (under alternative)
> 1-pt(cvGP, dfGP, ncp)                    # power
[1] 0.3348385

Это соответствует результату G * Power, который является отличной программой для этих вопросов. Он также отображает df, критическое значение, ncp, так что вы можете проверить все эти расчеты отдельно.

Изменить: Использование формулы Satterthwaite или Уэлча не сильно меняется (по-прежнему 0,33 *):

# Satterthwaite's formula
> var1  <- sd1^2
> var2  <- sd2^2
> num   <- (var1/n1 + var2/n2)^2
> denST <- var1^2/((n1-1)*n1^2) + var2^2/((n2-1)*n2^2)
> (dfST <- num/denST)
[1] 33.10309

> cvST <- qt(1-alpha, dfST)
> 1-pt(cvST, dfST, ncp)
[1] 0.3336495

# Welch's formula
> denW <- var1^2/((n1+1)*n1^2) + var2^2/((n2+1)*n2^2)
> (dfW <- (num/denW) - 2)
[1] 34.58763

> cvW   <- qt(1-alpha, dfW)
> 1-pt(cvW, dfW, ncp)
[1] 0.3340453

(обратите внимание, что я немного изменил имена некоторых переменных как t, dfи diffони также являются именами встроенных функций, также обратите внимание, что числитель вашего кода для dfнеправильный, он неуместен ^2, и его ^2слишком много, это должно быть ((sd1^2/n1) + (sd2^2/n2))^2)

Если вы в основном заинтересованы в вычислительной мощности (а не обучение через это делать вручную) и вы уже используете R , то смотрите на pwrупаковке и либо pwr.t.testили pwr.t2n.testфункций. (это может быть полезно для проверки ваших результатов, даже если вы делаете это вручную, чтобы учиться).