Esta práctica tiene dos partes. En la primera se utilizan los contrastes basados en la distribución \(\chi^2\) para responder algunas preguntas sobre la distribución de supervivientes en el hundimiento del Titanic. En la segunda se investigan algunas de las propiedades del contraste de Kolmogorov-Smirnov.

Contrastes basados en la distribución \(\chi^2\)

El fichero titanic.txt contiene la lista de pasajeros del Titanic junto con otra información acerca de los mismos. Agunos datos son faltantes, pero no en las variables que analizaremos en esta práctica.

Variable Descripción o valores
Name Nombre completo del pasajero
Survived Indicador de si sobrevivió al hundimiento: Dead, Alive
Boarded Puerto de embarque: Belfast, Cherbourg, Queenstown, Southampton
Class 1, 2, 3, Crew
MWC Man, Woman, Child
Age Edad
Adut_or_Chld Adult, Child
Sex Male, Female
Paid Precio del billete
Ticket_No Número del billete
Boat_or_Body Identificador del bote si sobrevivió, número identificativo entre corchetes si se recuperó el cadaver
Job Oficio
Class_Dept
Class_Full

El objetivo de esta sección es contrastar si la distribución de supervivientes fue la misma según clasificaciones de los pasajeros (por sexo o por clase). Para ello, llevaremos a cabo un contraste \(\chi^2\) de homogeneidad.

Para cargar los datos desde la web:

titanic <- read.delim("https://dasl.datadescription.com/download/data/3484")

Pinchando en el data frame titanic en la pestaña Environment (arriba a la derecha) podemos ver la tabla de datos.

Para hacer un contraste de homogeneidad de la distribución de supervivientes según el sexo del pasajero, lo primero que tenemos que calcular son las frecuencias observadas de las dos variables:

# Primero extraemos las observaciones de las variables Survived y Sex:
SurvivedSex = list(Survived = titanic$Survived,Sex = titanic$Sex)
# Tabla de contingencia de las frecuencias observadas:
SurvivedSexT = table(SurvivedSex)
SurvivedSexT
##         Sex
## Survived Female Male
##    Alive    359  353
##    Dead     130 1366

Finalmente llevamos a cabo el contraste de homogeneidad:

chisq.test(SurvivedSexT)
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  SurvivedSexT
## X-squared = 484.84, df = 1, p-value < 2.2e-16

Se observa que el p-valor es muy pequeño, por lo que se rechaza que la distribución de supervivientes fuera la misma entre hombres y mujeres para los niveles de significación habituales. Por defecto se obtiene la información básica. En realidad lo que calcula R es una lista cuya estructura se puede ver de la forma siguiente:

resultado <- chisq.test(SurvivedSexT)
ls.str(resultado)
## data.name :  chr "SurvivedSexT"
## expected :  num [1:2, 1:2] 158 331 554 1165
## method :  chr "Pearson's Chi-squared test with Yates' continuity correction"
## observed :  'table' int [1:2, 1:2] 359 130 353 1366
## p.value :  num 1.9e-107
## parameter :  Named int 1
## residuals :  'table' num [1:2, 1:2] 16.03 -11.06 -8.55 5.9
## statistic :  Named num 485
## stdres :  'table' num [1:2, 1:2] 22.1 -22.1 -22.1 22.1

Si solo nos interesa calcular el estadístico de contraste o necesitamos la matriz de frecuencias esperadas podemos calcularlos de la siguiente forma:

resultado$statistic
## X-squared 
##  484.8352
resultado$expected
##         Sex
## Survived   Female      Male
##    Alive 157.6848  554.3152
##    Dead  331.3152 1164.6848

Ejercicios

  1. Contrasta si la distribución de supervivientes fue homogénea según la clase (Class). Ya nos imaginamos que no, ¿verdad?

Contraste de Kolmogorov-Smirnov

En primer lugar vamos a escribir una función (ksnoest) que calcula el estadístico de Kolmogorov-Smirnov cuando estamos interesados en contrastar si nuestros datos siguen una distribución normal estándar:

ksnoest <- function(datos){
    y <- ks.test(datos,pnorm)$statistic
    return(y)
}

Para comprobar si funciona podemos generar una muestra normal de tamaño 20 y aplicar después la función:

ksnoest(rnorm(20))
##         D 
## 0.1807957

Supongamos que queremos contrastar la hipótesis nula de que los datos son normales (con valores arbitrarios de la media y la desviación típica). Una posibilidad (que, como veremos, no funciona) es estimar los parámetros de la normal y comparar la función de distribución empírica \(F_n\) con la función de distribución de una variable \(\mbox{N}(\hat{\mu},\hat{\sigma}^2)\). La siguiente función calcula el correspondiente estadístico de Kolmogorov-Smirnov:

ksest <- function(datos){
    mu <- mean(datos)
    stdev <- sd(datos)
    y <- ks.test(datos, pnorm, mean=mu, sd=stdev)$statistic
    return(y)
}

Para ver la diferencia entre las dos situaciones, generamos 1000 muestras de tamaño 20 de una N(0,1) y calculamos ambos estadísticos para cada una de ellas:

B <- 1000
n <- 500
datos <- matrix(rnorm(n*B), n)
test <- apply(datos, 2, ksest)
tnoest <- apply(datos, 2, ksnoest)

Mediante diagramas de cajas, podemos comparar las distribuciones del estadístico en cada caso:

boxplot(test,tnoest, names=c('Estimando','Sin estimar'))

Ejercicios

  1. Claramente las distribuciones de test y de tnoest son diferentes, por lo que no podemos usar las mismas tablas para hacer el contraste en las dos situaciones. ¿En cuál de los dos casos se obtienen en media valores menores? ¿Podrías dar una razón intuitiva?

  2. Imagina que estimamos los parámetros y usamos las tablas de la distribución del estadístico de Kolmogorov-Smirnov para hacer el contraste a nivel \(\alpha\). El verdadero nivel de significación, ¿es mayor o menor que \(\alpha\)?

  3. Para resolver el problema se ha estudiado la distribución del estadístico de Kolmogorov-Smirnov en el caso de muestras normales con parámetros estimados. Es lo que se conoce como contraste de normalidad de Kolmogorov-Smirnov-Lilliefors (KSL) (véase, por ejemplo, pag. 471 y Tabla 9 de Peña 2008). Según la tabla del estadístico KSL, el nivel crítico para \(\alpha=0.05\) y \(n=20\) es 0.192. Esto significa que el porcentaje de valores de test mayores que 0.192 en nuestra simulación debe ser aproximadamente del 5%. Compruébalo haciendo sum(test > 0.192)/B. Haz una pequeña simulación similar a la anterior para aproximar el nivel de significación del contraste KSL cuando se utiliza un valor crítico 0.12 para muestras de tamaño 40.

  4. Genera \(B=10000\) muestras de tamaño \(n=30\) de una distribución exponencial de media 1 y utilízalas para determinar en este caso la potencia aproximada del test de Kolmogorov-Smirnov con \(\alpha=0.05\) para \(H_0: X\sim \mbox{N}(1,1)\). (El comando rexp() puede utilizarse para generar los datos exponenciales).

Referencias

Peña, D. (2008). Fundamentos de Estadística. Alianza Editorial.