Rastreando las suposiciones hechas por la función `ttest_ind()` de SciPy

Question

Estoy tratando de escribir mi propio código Python para calcular t-statistics y p-values ​​para una y dos pruebas t independientes con cola. Puedo usar la aproximación normal, pero por el momento estoy tratando de usar la distribución t. No he podido comparar los resultados de la biblioteca de estadísticas de SciPy en mis datos de prueba. Podría usar un par de ojos frescos para ver si estoy cometiendo un error tonto en alguna parte.

Tenga en cuenta que esto es cross-posted from Cross-Validated porque ha estado ahí por un tiempo sin respuestas, así que pensé que no estaría de más recibir algunas opiniones de los desarrolladores de software. Estoy tratando de entender si hay un error en el algoritmo que estoy usando, que debería reproducir el resultado de SciPy. Este es un algoritmo simple, por lo que es desconcertante por qué no puedo encontrar el error.

Mi código:

import numpy as np 
import scipy.stats as st 

def compute_t_stat(pop1,pop2): 

    num1 = pop1.shape[0]; num2 = pop2.shape[0]; 

    # The formula for t-stat when population variances differ. 
    t_stat = (np.mean(pop1) - np.mean(pop2))/np.sqrt(np.var(pop1)/num1 + np.var(pop2)/num2) 

    # ADDED: The Welch-Satterthwaite degrees of freedom. 
    df = ((np.var(pop1)/num1 + np.var(pop2)/num2)**(2.0))/( (np.var(pop1)/num1)**(2.0)/(num1-1) + (np.var(pop2)/num2)**(2.0)/(num2-1)) 

    # Am I computing this wrong? 
    # It should just come from the CDF like this, right? 
    # The extra parameter is the degrees of freedom. 

    one_tailed_p_value = 1.0 - st.t.cdf(t_stat,df) 
    two_tailed_p_value = 1.0 - (st.t.cdf(np.abs(t_stat),df) - st.t.cdf(-np.abs(t_stat),df))  


    # Computing with SciPy's built-ins 
    # My results don't match theirs. 
    t_ind, p_ind = st.ttest_ind(pop1, pop2) 

    return t_stat, one_tailed_p_value, two_tailed_p_value, t_ind, p_ind 

Actualización:

Después de leer un poco más en la prueba t de Welch, vi que debo utilizar la fórmula de Welch-Satterthwaite para calcular los grados de libertad. Actualicé el código anterior para reflejar esto.

Con los nuevos grados de libertad, obtengo un resultado más cercano. Mi valor de p de dos lados está desactivado en aproximadamente 0.008 desde la versión de SciPy ... pero sigue siendo un error demasiado grande, así que todavía debo estar haciendo algo incorrecto (o las funciones de distribución de SciPy son muy malas, pero es difícil de creer solo tienen una precisión de 2 decimales).

Segunda actualización:

Sin dejar de probar cosas, pensé que tal vez la versión de SciPy calcula automáticamente la aproximación normal a la distribución t cuando los grados de libertad son lo suficientemente alta (aproximadamente> 30). Así que volví a ejecutar mi código utilizando la distribución Normal en su lugar, y los resultados calculados están más alejados de SciPy que cuando uso la distribución t.

Bono pregunta :) (teoría estadística relacionada más, no dude en pasar por alto)

Además, el estadístico t es negativo. Me preguntaba qué significa esto para la prueba t unilateral. ¿Esto significa típicamente que debería estar mirando en la dirección del eje negativo para la prueba? En mis datos de prueba, la población 1 es un grupo de control que no recibió un determinado programa de capacitación laboral. La población 2 sí la recibió, y los datos medidos son diferencias salariales antes/después del tratamiento.

Así que tengo alguna razón para pensar que la media para la población 2 será mayor. Pero desde el punto de vista de la teoría estadística, no parece correcto inventar una prueba de esta manera. ¿Cómo podría haber sabido comprobar (en una prueba unilateral) en la dirección negativa sin confiar en el conocimiento subjetivo sobre los datos? ¿O es solo una de esas cosas frecuentistas que, si bien no son filosóficamente rigurosas, deben hacerse en la práctica?

Answer 1

Al utilizar la función integrada SciPy source(), pude ver una impresión del código fuente de la función ttest_ind(). Basado en el código fuente, el SciPy incorporado realiza la prueba t suponiendo que las varianzas de las dos muestras son iguales. No está utilizando los grados de libertad de Welch-Satterthwaite. SciPy asume varianzas iguales pero no establece esta suposición.

sólo quiero señalar que, de manera crucial, Es por eso que no sólo debe confiar en funciones de la biblioteca. En mi caso, realmente necesito la prueba t para poblaciones de varianzas desiguales, y el ajuste de grados de libertad puede ser importante para algunos de los conjuntos de datos más pequeños en los que voy a ejecutar esto.

Como mencioné en algunos comentarios, la discrepancia entre mi código y SciPy es de aproximadamente 0.008 para tamaños de muestra entre 30 y 400, y luego va lentamente a cero para tamaños de muestra más grandes. Este es un efecto del término adicional (1/n1 + 1/n2) en el denominador t-estadístico de varianzas iguales. Exactitud, esto es muy importante, especialmente para muestras pequeñas. Definitivamente me confirma que necesito escribir mi propia función. (Posiblemente haya otras mejores bibliotecas de Python, pero al menos debería ser conocido. Francamente, es sorprendente que esto no esté en ningún lugar desde el principio en la documentación de SciPy para ttest_ind()).

Answer 2

Parece que olvidó ** 2 en el numerador de su df. Los grados de libertad Welch-Satterthwaite.

df = (np.var(pop1)/num1 + np.var(pop2)/num2)/( (np.var(pop1)/num1)**(2.0)/(num1-1) + (np.var(pop2)/num2)**(2.0)/(num2-1)) 

debería ser:

df = (np.var(pop1)/num1 + np.var(pop2)/num2)**2/( (np.var(pop1)/num1)**(2.0)/(num1-1) + (np.var(pop2)/num2)**(2.0)/(num2-1)) 

Answer 3

No está calculando la varianza de la muestra , pero en su lugar se utiliza varianzas poblacionales. La varianza muestral se divide por n-1, en lugar de n. np.var tiene un argumento opcional llamado ddof por razones similares a esta.

Esto debe darle a su resultado esperado:

import numpy as np 
import scipy.stats as st 

def compute_t_stat(pop1,pop2): 

    num1 = pop1.shape[0] 
    num2 = pop2.shape[0]; 
    var1 = np.var(pop1, ddof=1) 
    var2 = np.var(pop2, ddof=1) 

    # The formula for t-stat when population variances differ. 
    t_stat = (np.mean(pop1) - np.mean(pop2))/np.sqrt(var1/num1 + var2/num2) 

    # ADDED: The Welch-Satterthwaite degrees of freedom. 
    df = ((var1/num1 + var2/num2)**(2.0))/((var1/num1)**(2.0)/(num1-1) + (var2/num2)**(2.0)/(num2-1)) 

    # Am I computing this wrong? 
    # It should just come from the CDF like this, right? 
    # The extra parameter is the degrees of freedom. 

    one_tailed_p_value = 1.0 - st.t.cdf(t_stat,df) 
    two_tailed_p_value = 1.0 - (st.t.cdf(np.abs(t_stat),df) - st.t.cdf(-np.abs(t_stat),df))  


    # Computing with SciPy's built-ins 
    # My results don't match theirs. 
    t_ind, p_ind = st.ttest_ind(pop1, pop2) 

    return t_stat, one_tailed_p_value, two_tailed_p_value, t_ind, p_ind 

PS: SciPy es de código abierto y en su mayoría implementado con Python. Usted podría haber comprobado el código fuente para ttest_ind y descubrir su error usted mismo.

Por el lado de la bonificación: Usted no elija elija el lado de la prueba de una cola mirando su valor t. Usted lo decide de antemano con su hipótesis. Si su hipótesis nula es que los medios son iguales y su hipótesis alternativa es que la segunda media es más grande, entonces su cola debe estar en el lado izquierdo (negativo). Debido a que los valores suficientemente pequeños (negativos) de su valor t indicarían que la hipótesis alternativa es más probable que sea verdadera en lugar de la hipótesis nula.