Potencia para dos muestras de prueba t

Estoy tratando de entender el cálculo de potencia para el caso de la prueba t de dos muestras independientes (no suponiendo variaciones iguales, entonces utilicé Satterthwaite).

Aquí hay un diagrama que encontré para ayudar a comprender el proceso:

Así que supuse que dado lo siguiente sobre las dos poblaciones y los tamaños de muestra:

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

Podría calcular el valor crítico bajo el valor nulo relacionado con tener 0.05 de probabilidad de cola superior:

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

y luego calcule la hipótesis alternativa (que para este caso aprendí es una "distribución t no central"). Calculé beta en el diagrama anterior usando la distribución no central y el valor crítico encontrado anteriormente. Aquí está el script completo en 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)

Esto da un valor de potencia de 0.4935132.

¿Es este el enfoque correcto? Me parece que si uso otro software de cálculo de potencia (como SAS, que creo que he configurado de manera equivalente a mi problema a continuación) obtengo otra respuesta (de SAS es 0.33).

CÓDIGO SAS:

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

En última instancia, me gustaría obtener una comprensión que me permitiera ver simulaciones para procedimientos más complicados.

EDITAR: Encontré mi error. debería haber sido

1 pt (CV, df, ncp) NO 1 pt (t, df, ncp)

Estás cerca, aunque se requieren algunos pequeños cambios:

• La verdadera diferencia en las medias generalmente se toma como , no al revés.$\mu_{2} - \mu_{1}$
• G * Power usa como grados de libertad para la distribución en este caso (diferentes variaciones, mismos tamaños de grupo), siguiendo una sugerencia de Cohen como se explica aquí$n_{1}+n_{2} - 2$$t$
• SAS podría usar la fórmula de Welch o la fórmula de Satterthwaite para las variaciones desiguales dadas por df (que se encuentra en este pdf que citó ), con solo 2 dígitos significativos en el resultado que uno no puede decir (ver más abajo)

Con n1, n2, mu1, mu2, sd1, sd2lo definido en su pregunta:

> 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

Esto coincide con el resultado de G * Power, que es un gran programa para estas preguntas. También muestra el valor crítico, df, ncp, por lo que puede verificar todos estos cálculos por separado.

Editar: el uso de la fórmula de Satterthwaite o la fórmula de Welch no cambia mucho (todavía 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

(tenga en cuenta que cambié ligeramente algunos nombres de variables como t, dfy difftambién son nombres de funciones integradas, también tenga en cuenta que el numerador de su código dfes incorrecto, tiene un lugar incorrecto ^2y uno ^2demasiado, debería ser ((sd1^2/n1) + (sd2^2/n2))^2)

Si está interesado principalmente en calcular el poder (en lugar de aprender haciéndolo a mano) y ya está utilizando R, mire el pwrpaquete y las funciones pwr.t.testo pwr.t2n.test. (estos pueden ser buenos para verificar sus resultados, incluso si lo hace a mano para aprender).