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.
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
Class
). Ya nos imaginamos que no, ¿verdad?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'))
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?
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\)?
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.
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.