Moc dla testu t dwóch próbek

Próbuję zrozumieć obliczenia mocy dla przypadku dwóch niezależnych próbnych testów t (nie zakładając równych wariancji, więc użyłem Satterthwaite).

Oto schemat, który znalazłem, aby pomóc zrozumieć proces:

Więc założyłem, że biorąc pod uwagę następujące informacje o dwóch populacjach i biorąc pod uwagę rozmiary próbek:

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

Mógłbym obliczyć wartość krytyczną poniżej zera odnoszącą się do prawdopodobieństwa górnego ogona 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

a następnie obliczyć alternatywną hipotezę (którą dla tego przypadku nauczyłem się jest „niecentralny rozkład t”). Obliczyłem beta na powyższym schemacie, używając rozkładu niecentralnego i wartości krytycznej podanej powyżej. Oto pełny skrypt w języku 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)

Daje to wartość mocy 0,4935132.

Czy to jest właściwe podejście? Zauważyłem, że jeśli używam innego oprogramowania do obliczania mocy (takiego jak SAS, które, jak myślę, mam skonfigurowany równoważnie do mojego problemu poniżej), otrzymuję inną odpowiedź (z SAS to 0.33).

KOD SAS:

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

Ostatecznie chciałbym uzyskać zrozumienie, które pozwoliłoby mi spojrzeć na symulacje w celu uzyskania bardziej skomplikowanych procedur.

EDYCJA: Znalazłem swój błąd. powinien był być

1-pt (CV, df, ncp) NOT 1-pt (t, df, ncp)

Jesteś blisko, potrzebne są jednak niewielkie zmiany:

• $\mu_{2} - \mu_{1}$
• $n_{1}+n_{2} - 2$$t$
• SAS może użyć formuły Welcha lub formuły Satterthwaite dla df, biorąc pod uwagę nierówne wariancje (znalezione w cytowanym pliku pdf ) - z jedynie 2 cyframi znaczącymi w wyniku, których nie można powiedzieć (patrz poniżej)

Z n1, n2, mu1, mu2, sd1, sd2jak zdefiniowano w pytaniu:

> 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

Odpowiada to wynikowi G * Power, który jest świetnym programem na te pytania. Wyświetla również df, wartość krytyczną, ncp, więc możesz sprawdzić wszystkie te obliczenia osobno.

Edycja: użycie formuły Satterthwaite lub Welcha niewiele się zmienia (wciąż 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

(zwróć uwagę, że nieznacznie zmieniłem niektóre nazwy zmiennych jako t, dfi diffsą to także nazwy wbudowanych funkcji, zauważ również, że licznik twojego kodu dfjest nieprawidłowy, ma niewłaściwe położenie ^2i ^2powinien być o jeden za dużo ((sd1^2/n1) + (sd2^2/n2))^2)

Jeśli jesteś zainteresowany obliczeniem mocy (zamiast uczyć się, jak to robić ręcznie) i już używasz R, spójrz na pwrpakiet i funkcje pwr.t.testlub pwr.t2n.test. (mogą być przydatne do weryfikacji wyników, nawet jeśli robisz to ręcznie, aby się uczyć).