Francisco Guijarro

Universidad Politécnica de Valencia

Creative Commons Attribution 4.0 International License (CC BY 4.0)

1 Primeros pasos

1.1 Una primera aproximación a través de los datos del naufragio del Titanic

En esta sección se realiza una breve introducción a la programación en el entorno R a partir de un pequeño ejemplo, que servirá para familiarizarse con algunas funciones básicas. Para ello, se leerá una primera fuente de datos a partir de una librería de R.

Las librerías de R son conjuntos de funciones y datos que facilitan la programación de diferentes tareas y análisis. A través de las librerías, los usuarios de R pueden acceder a código que ya está previamente programado, por lo que se agiliza la utilización de este software tanto para su uso directo como para su modificación si se requiere incorporar nuevas funcionalidades.

Emplear estas librerías tiene al menos dos ventajas. Por una parte, aceleramos el proceso de programación, ya que nos ahorramos tener que programar algunas funciones que otros han programado por nosotros. Por otro lado, las librerías siempre van a contener funciones que han sido previamente depuradas por el equipo de desarrollo de R. Es decir, no van a contener error de programación, lo que también acelerará la fase de desarrollo de cualquier programa que queramos codificar.

Cuando las librerías no están previamente instaladas en el ordenador, cosa que ocurre habitualmente con la mayoría de ellas, dichas librerías deben instalarse en el disco duro del ordenador con el comando install.packages, para después cargarlas en memoria con el comando library. Por lo tanto, tener instalado R no significa que se haya instalado todas las librerías disponibles. La instalación de una librería sólo debe hacerse una vez (salvo actualizaciones), pero para que sus funciones y datos estén disponibles al programador es necesario cargarlas en memoria. Para ello se emplea el comando library. Veamos un primer ejemplo con la librería carData:

install.packages("carData")
library(carData)

Ahora ya tenemos acceso al contenido de la librería. Por ejemplo, podemos leer en memoria una de las tablas que incluye, denominada TitanicSurvival. Para ello se puede emplear la función data. Tras ejecutar la siguiente línea de código, el panel Environment se actualizará con los datos correspondientes al naufragio del célebre trasatlántico Titanic:

data(TitanicSurvival)

La función head permite conocer las primeras observaciones de una estructura de datos (el encabezado), mientras que tail muestra las últimas observaciones de los datos:

head(TitanicSurvival)
##                                 survived    sex     age passengerClass
## Allen, Miss. Elisabeth Walton        yes female 29.0000            1st
## Allison, Master. Hudson Trevor       yes   male  0.9167            1st
## Allison, Miss. Helen Loraine          no female  2.0000            1st
## Allison, Mr. Hudson Joshua Crei       no   male 30.0000            1st
## Allison, Mrs. Hudson J C (Bessi       no female 25.0000            1st
## Anderson, Mr. Harry                  yes   male 48.0000            1st
tail(TitanicSurvival)
##                           survived    sex  age passengerClass
## Yousseff, Mr. Gerious           no   male   NA            3rd
## Zabour, Miss. Hileni            no female 14.5            3rd
## Zabour, Miss. Thamine           no female   NA            3rd
## Zakarian, Mr. Mapriededer       no   male 26.5            3rd
## Zakarian, Mr. Ortin             no   male 27.0            3rd
## Zimmerman, Mr. Leo              no   male 29.0            3rd

En nuestro caso, puedes ver que los registros de la base de datos se han ordenado alfabéticamente por el nombre del pasajero.

Con la función str podemos conocer cómo se estructuran los datos:

str(TitanicSurvival)
## 'data.frame':    1309 obs. of  4 variables:
##  $ survived      : Factor w/ 2 levels "no","yes": 2 2 1 1 1 2 2 1 2 1 ...
##  $ sex           : Factor w/ 2 levels "female","male": 1 2 1 2 1 2 1 2 1 2 ...
##  $ age           : num  29 0.917 2 30 25 ...
##  $ passengerClass: Factor w/ 3 levels "1st","2nd","3rd": 1 1 1 1 1 1 1 1 1 1 ...

La función str nos dice que los datos de TitanicSurvival se han agrupado en lo que denomina un data.frame; esto es, una estructura de datos. Veremos, más adelante, éste y otro tipo de datos que maneja R.

La estructura de los datos nos permite ver el número de observaciones o filas, así como el número de columnas o variables. En el caso de los datos del Titanic, podemos ver que tenemos 1309 filas (registros, observaciones) y 4 variables. Junto al nombre de cada variable aparece el tipo: tenemos una variable numérica (age) y 3 variables tipo factor (survived, sex, passengerClass), que analizaremos en detalle más adelante.

Para acceder a las variables del dataframe debemos emplear el símbolo $:

head(TitanicSurvival$age)
## [1] 29.0000  0.9167  2.0000 30.0000 25.0000 48.0000

Otras funciones de interés son summary, que obtiene el resumen de la estructura de datos, y help (también accesible meidante el símobolo ?), que permite consultar la ayuda de cualquier función de R:

summary(TitanicSurvival)
##  survived      sex           age          passengerClass
##  no :809   female:466   Min.   : 0.1667   1st:323       
##  yes:500   male  :843   1st Qu.:21.0000   2nd:277       
##                         Median :28.0000   3rd:709       
##                         Mean   :29.8811                 
##                         3rd Qu.:39.0000                 
##                         Max.   :80.0000                 
##                         NA's   :263
help(TitanicSurvival)
?TitanicSurvival

Ejercicio 1

¿Cuántas observaciones (filas) tiene la tabla? ¿Cuántas variables?

Ejercicio 2

En la tabla de supervivientes del Titanic aparecen diferentes tipos de variables. ¿Puedes diferenciarlos? ¿Cómo se tratarán desde un punto de vista estadístico ¿Qué tipo de preguntas nos podríamos hacer a partir de la información recogida en la tabla? Para responder a esta pregunta puedes ver el resultado obtenido con las funciones head y str

También podemos obtener un resumen visual de las variables utilizando la función plot:

plot(TitanicSurvival)

plot(TitanicSurvival$age)

Ejercicio 3

Consulta la ayuda de la función plot y dibuja la supervivencia de los pasajeros en función de la clase en la que viajaban:

Puedes comprobar que el ancho de los 3 niveles correspondientes a la clase (medido a través del eje X) es proporcional al número de pasajeros que viajaban en cada una de esas clases. Es decir, la clase con mayor número de pasajeros era tercera clase, mientras que la clase con menos afluecia de pasajeros fue segunda clase.

Ejercicio 4

¿Crees que la probabilidad de supervivencia pudo estar relacionada con la clase en la que se viajaba?

Lo observado en el gráfico, puede ser también analizado numéricamente. Por ejemplo, podemos conocer el número de personas que sobrevivieron frente a las que no para cada una de las 3 clases de billetes del Titanic a través de una tabla de contingencia, con la función table:

table(TitanicSurvival$passengerClass, TitanicSurvival$survived)
##      
##        no yes
##   1st 123 200
##   2nd 158 119
##   3rd 528 181

A tenor de estos datos, puede deducirse que los viajeros de primera clase tuvieron mayores probabilidades de sobrevivir que los de segunda o tercera. Pero en estadística este tipo de afirmaciones sólo pueden realizarse tras comprobarlo mediante el correspondiente test (y con el consiguiente nivel de confianza).

Los datos de la anterior tabla de contingencia son los valores observados. Para analizar si existe una relación entre ambas variables, llevaremos a caba un análisis Chi-cuadrado (\(\chi^2\)), donde contrastaremos las diferencias entre el número de pasajeros observados en la tabla y los valores teóricos que deberíamos haber observado suponiendo que ambas variables fueran independientes.

Si asumimos que las variables son independientes, ¿cuántos de los 1.309 pasajeros de la tabla anterior deberíamos tener en cada una de las celdas de la tabla de contingencia?

No Yes
1st 218,17 218,17
2nd 218,17 218,17
3rd 218,17 218,17

En esta tabla hemos supuesto que los 1.309 pasajeros se distribuyen de forma homogénea entre las 6 celdas: \(\frac{1309}{6}=218,17\).

Sin embargo, enseguida podemos darnos cuenta de que este enfoque tiene un sesgo respecto de la tasa de supervivencia. Si en la tabla de contingencia sumamos las personas que perecieron en el naufragio (columna No) contabilizamos un total de 809 pasajeros (\(123+158+528=809\)), mientras que las personas que sobrevivieron ascienden a 500 pasajeros (\(200+119+181=500\)). Es decir, que en el naufragio del Titanic fue más probable perecer que sobrevivir.

Si lo traducimos en porcentaje, los supervivientes representan un \(\frac{500}{809+500}=38,2%\) frente a los no supervivientes del \(\frac{809}{809+500}=61,8%\). Si balanceamos los datos de la tabla para respetar estos porcentajes, manteniendo la suma por filas (\(218,17+218,17=436,33\)), los nuevos valores serían los siguientes:

No Yes
1st 269,67 166,67
2nd 269,67 166,67
3rd 269,67 166,67

Los valores de esta tabla se obtienen multiplicando las tasas de no supervivencia / superviviencia por el valor 218,17: \(61,8\% \times 436,33= 269,67\) y \(38,2\% \times 436,33 = 166,67\).

Pero de nuevo podemos pensar que estos valores vuelven a estar sesgados, ya que los pasajeros no se distribuyeron de manera uniforme entre las 3 clases.

En primera clase viajaron el 24,7% (\(\frac{323}{1309}=24,7\%\)), mientras que en segunda viajaron el 21,2% (\(\frac{277}{1309}=21,2\%\)) y en tercera el 54,2% (\(\frac{709}{1309}=54,2\%\)). Esto se traduce en que por filas no deberíamos tener 436,33 pasajeros, como hemos supuesto en el caso anterior, sino los siguientes valores:

Pasajeros
1st 323
2nd 277
3rd 709

Y multiplicando estos valores por 61,8% y 38,2% tenemos:

No Yes
1st 199,62 123,38
2nd 171,19 105,81
3rd 438,18 270,82

Que serán los valores teóricos que compararemos con los valores observados para realizar el contraste.

Los valores de la última tabla se podrían haber obtenido de forma directa aplicando las siguientes ecuaciones:

\(e_{ij}=p_{ij} \times n\)

\(p_{ij}=\frac{n_{i*}}{n} \times \frac{n_{*j}}{n}\)

donde \(n\) es el número de pasajeros totales, \(n_{i*}\) es la suma de pasajeros en la fila i-ésima (clase), \(n_{*j}\) es la suma de pasajeros en la columna j-ésima (supervivencia), \(p_{ij}\) denota la probabilidad de pertenecer a la celda \(ij\) de la tabla, mientras que \(e_{ij}\) denota el número de pasajeros esperado (valor teórico) para la celda \(ij\).

El estadístico \(\chi^2\) se calcula a partir de la siguiente fórmula:

\(\chi_{(k-1) \times (m-1)}^2 = \sum_{i=1}^k \sum_{j=1}^m \frac{(n_{ij}-e_{ij})^2}{e_{ij}}\)

donde \(k\) es el número de niveles en las clases de pasajeros (3) y \(m\) el número de niveles en la variable de superviviencia (2). Por lo tanto, el valor de la \(\chi^2\) a buscar en tablas es con 2 y 1 grados de libertad, respectivamente:

\(\chi_{2}^2 = \frac{(123-199,62)^2}{199,62} + \frac{(200-123,38)^2}{123,38} + \frac{(158-171,19)^2}{171,19} + \frac{(119-105,81)^2}{105,81} + \frac{(528-438,18)^2}{438,18} + \frac{(181-270,82)^2}{270,82} = 127,86\)

Mientras que el valor de la distribución en tablas para un nivel de confianza del 95% es:

qchisq(.95, df = 2)
## [1] 5.991465

Como el valor observado (\(127,86\)) es mayor que el valor en tablas 5.99, podemos concluir que se rechaza la hipótesis nula de igualdad en las distribuciones. Es decir, que el la supervivencia y la clase de los pasajeros no son variables independientes.

Como puedes ver, la probabilidad de encontrar un valor superior a \(127,86\) es prácticamente 0 (p-value):

curve(dchisq(x, df = 2), from = 0, to = 20)

curve(dchisq(x, df = 5), from = 0, to = 20)

Podríamos llegar al mismo resultado utilizando la función chisq.test que trae incorporada R:

analisis <- chisq.test(table(TitanicSurvival$passengerClass, TitanicSurvival$survived))
analisis
## 
##  Pearson's Chi-squared test
## 
## data:  table(TitanicSurvival$passengerClass, TitanicSurvival$survived)
## X-squared = 127.86, df = 2, p-value < 2.2e-16

Estos resultados confirman que hubo diferencias en la probabilidad de supervivencia según la clase en la que se viajaba. Ello no implica que las diferencias fueran significativas para cada una de las 3 clases, sino que al menos para una de ella la probabilidad de sobrevivir fue diferente al resto.

Ejercicio 5

Comprueba que la variable analisis es una variable tipo lista1, y que se puede acceder a los elementos que la componen mediante el símbolo $. Por ejemplo, analisis$statistic ofrece el valor calculado para el estadístico \(\chi^2\). Comprueba el resto de elementos de analisis, y a qué valores o tablas de las elaboradas anteriormente se corresponden.

Ejercicio 6

Una vez analizadas visualmente las diferencias en la tasa de supervivencia de los pasajeros en función de la clase en la que viajaban, y su posterior análisis estadístico para conocer si esas diferencias son estadísticamente significativas, ¿podrías repetir el análisis para comprobar si hubo diferencias en la tasa de supervivencia según el género de los pasajeros? Deberías obtener unos resultados como los siguientes:

## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  table(TitanicSurvival$sex, TitanicSurvival$survived)
## X-squared = 363.62, df = 1, p-value < 2.2e-16

Para terminar con el análisis básico de la tabla de supervivientes, e introducir otras funciones genéricas de R, vamos a calcular algunos estadísticos básicos de la variable edad (TitanicSurvival$age) y cómo se relaciona con otras variables. Podemos, por ejemplo, obtener una estadística básica de la edad de los pasajeros del Titanic, que incluya la media junto con los cuartiles de su distribución:

summary(TitanicSurvival$age)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##  0.1667 21.0000 28.0000 29.8811 39.0000 80.0000     263

Vemos que la media no coincide exactamente con la mediana, lo que nos puede estar informando sobre una ligera asimetría en la distribución. Lo comprobamos gráficamente mediante un histograma:

hist(TitanicSurvival$age)

Al hacer el resumen de la distribución con la función summary, se ha podido comprobar que no se disponía de la edad de todos los pasajeros. Hay 263 pasajeros cuya edad aparece como NA (Not Available). Los valores faltantes representan un problema cuando queremos evaluar algunas funciones o test estadísticos. Por ejemplo, si queremos conocer la edad media a través de la función mean, comprobamos que dicha función no puede calcular la media cuando aparece algún valor NA:

mean(TitanicSurvival$age)
## [1] NA

Sin embargo, podemos excluir aquellos pasajeros con edad NA y calcular la media sobre las observaciones restantes a través de la función na.exclude:

mean(na.exclude(TitanicSurvival$age))
## [1] 29.88113

Comprobamos cómo el resultado obtenido coincide con el que nos proporcionaba la función summary. Esto es, que la función summary obtiene la media empleando el mismo código que nosotros (excluyendo las observaciones NA). También podemos completar el análisis obteniendo otros estadísticos de dispersion:

sd(na.exclude(TitanicSurvival$age))
## [1] 14.4135
var(na.exclude(TitanicSurvival$age))
## [1] 207.749
var(TitanicSurvival$age, na.rm = TRUE)
## [1] 207.749
range(TitanicSurvival$age, na.rm = TRUE)
## [1]  0.1667 80.0000

Observa cómo el parámetro na.rm tiene el mismo efecto que aplicar la función na.exclude sobre los datos originales.

1.2 Diferencias de superviviencia en el naufragio del Titanic: adulto vs. niño

Ya hemos comprobado anteriormiente que la tasa de supervivencia en el naufragio estuvo relacionada con la clase en la que viajaban los pasajeros y con el género de los mismos. Pero, ¿qué ocurrió con los niños? ¿se cumplió aquello de “las mujeres y los niños primero”?

Titanic - Women and children first - Scene

En aquél viaje se consideró menor de edad a todo pasajero con menos de 12 años, por lo que conviene entonces analizar la tasa de superviviencia diferenciando entre pasajeros por debajo de los 12 años (niños) frente a pasajeros con 12 o más años (considerados adultos). Los niños contaban con un descuento del 50% sobre los adultos. Además, en el caso de bebés (menos de 1 año) su billete era gratis. Nuestro análisis simplemente diferenciará adultos (12 años o más) del resto, para lo que crearemos una nueva variable en el data.frame a la que denominaremos adult:

TitanicSurvival$adult <- TitanicSurvival$age >= 12
table(TitanicSurvival$adult)
## 
## FALSE  TRUE 
##    91   955

La primera línea de código crea la mencionada variable adult. Para ello, evalúa la condición TitanicSurvival$age >= 12. Esto es, analiza cada una de las observaciones o filas de TitanicSurvival y comprueba si la edad es 12 años o más, en cuyo caso asigna el valor TRUE; si la edad está por debajo de los 12 años, el resultado de TitanicSurvival$age >= 12 es FALSE. Una cuestión a tener en cuenta es que los pasajeros para los que no figura su edad tendrán un valor NA en la nueva variable creada. La segunda línea de código emplea la función table para obtener una estadística o tabla con los dos posibles valores de la nueva variable recién creada. También podemos ver el encabezado de la tabla para comprobar que los valores obtenidos son coherentes con el objetivo de diferenciar adultos de niños:

head(TitanicSurvival)
##                                 survived    sex     age passengerClass adult
## Allen, Miss. Elisabeth Walton        yes female 29.0000            1st  TRUE
## Allison, Master. Hudson Trevor       yes   male  0.9167            1st FALSE
## Allison, Miss. Helen Loraine          no female  2.0000            1st FALSE
## Allison, Mr. Hudson Joshua Crei       no   male 30.0000            1st  TRUE
## Allison, Mrs. Hudson J C (Bessi       no female 25.0000            1st  TRUE
## Anderson, Mr. Harry                  yes   male 48.0000            1st  TRUE

Ejercicio 7

Una vez creada la variable adult, comprueba si existen diferencias en la superviviencia a tenor del valor de la misma. Utiliza la función table para conocer el cruce entre las variables adult y survived:

##        
##          no yes
##   FALSE  40  51
##   TRUE  579 376
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  table(TitanicSurvival$adult, TitanicSurvival$survived)
## X-squared = 8.882, df = 1, p-value = 0.00288

El resultado parece confirmar que hubo diferencias en la tasa de superviviencia entre adultos y niños (a favor de estos últimos), aunque el p-valor obtenido indica que la significación estadística de este resultado no es tan grande como en el caso de mujeres frente a hombres.

El último análisis consistirá en tener en cuenta todos los factores que pudieron incidir en la tasa de supervivencia de los pasajeros, en lugar de analizarlos de forma aislada. Para ello emplearemos un análisis de regresión. La variable dependiente será survived, mientras que utilizaremos sex, age y passengerClass como variables independientes. Puesto que la varaible dependiete sólo puede tomar dos valores, sobrevivió (1) o no sobrevivió (0), llevaremos a cabo una regresión lineal generalizada donde la familia de la variable dependiente será la binomial:

regresion <- glm(survived ~ sex+age+passengerClass, data = TitanicSurvival, family = "binomial")
summary(regresion)
## 
## Call:
## glm(formula = survived ~ sex + age + passengerClass, family = "binomial", 
##     data = TitanicSurvival)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.6399  -0.6979  -0.4336   0.6688   2.3964  
## 
## Coefficients:
##                    Estimate Std. Error z value Pr(>|z|)    
## (Intercept)        3.522074   0.326702  10.781  < 2e-16 ***
## sexmale           -2.497845   0.166037 -15.044  < 2e-16 ***
## age               -0.034393   0.006331  -5.433 5.56e-08 ***
## passengerClass2nd -1.280570   0.225538  -5.678 1.36e-08 ***
## passengerClass3rd -2.289661   0.225802 -10.140  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1414.62  on 1045  degrees of freedom
## Residual deviance:  982.45  on 1041  degrees of freedom
##   (263 observations deleted due to missingness)
## AIC: 992.45
## 
## Number of Fisher Scoring iterations: 4

Observa cómo la primera línea de código calcula la regresión entre la variable dependiente survived y el resto de variables independientes. Al ejecutar esta línea de código, se creará una nueva variable en la pestaña Environmentcon el nombre regresion. La función summary nos indicará la significación estadística de cada variable.

Ejercicio 8

¿Te parece que todas las variables son igual de importantes en la explicación de la supervivencia del naufragio? ¿Cuál te parece más relevante y cuál menos?

En este tipo de regresiones no tenemos el estadístico \(R^2\) que nos informa sobre la bondad del modelo, debido a que la variable dependiente es dicotómica (survived). Sin embargo, podemos utilizar una matriz de confusión para tener una idea del “grado de acierto” del modelo de regresión:

library(regclass)
confusion_matrix(regresion)
##            Predicted no Predicted yes Total
## Actual no           520            99   619
## Actual yes          126           301   427
## Total               646           400  1046

De la anterior matriz podemos concluir lo siguiente: De las 619 personas que no sobrevivieron, el modelo acierta en 520 casos (el 84.0%), y falla en el resto. Respecto de las personas que sobrevivieron (427), el modelo acierta en 301; esto es, en el 70.5% de los casos. Si tenemos en cuenta el número de casos acertados 520 + 301 sobre el total 1046, el porcentaje general de acierto es 78.5%.

También podemos visualizar la tabla en términos relativos, lo que clarifica los resultados:

confusion <- confusion_matrix(regresion)
confusion/confusion[3, 3]
##            Predicted no Predicted yes     Total
## Actual no     0.4971319    0.09464627 0.5917782
## Actual yes    0.1204589    0.28776291 0.4082218
## Total         0.6175908    0.38240918 1.0000000

Ejercicio 9

En el anterior ejemplo, hemos calculado la matriz de confusión en términos relativos. Para ello, primero se ha guardado el resultado de confusion_matrix(regresion) en la nueva variable confusion, para después dividir todos los elementos de dicha variable por el valor que ocupa la posición [3, 3]. Otra alternativa sería haber calculado el mismo resultado con confusion_matrix(regresion) / confusion_matrix(regresion)[3, 3]. ¿Qué opción te parece mejor y por qué?

2 Ejercicios

En este apartado vamos a repasar las funciones vistas durante el capítulo con la ayuda de un nuevo conjunto de datos: Arrests:

data(Arrests)

Según la ayuda de este dataframe help(Arrests), la base de datos incluye “Data on police treatment of individuals arrested in Toronto for simple possession of small quantities of marijuana. The data are part of a larger data set featured in a series of articles in the Toronto Star newspaper.”; con 5226 observaciones y las siguientes variables: released, colour, year, age, sex, employed, citizen, checks. Nuestro objetivo será explicar la variable released: “Whether or not the arrestee was released with a summons”. La variable indica si el detenido fue liberado o no.

Ejercicio 10

Obtén el encabezado, la estructura de datos y un resumen de Arrests. ¿Qué variable intentaremos explicar a partir del resto? Puedes hacer hipótesis sobre la relación (positiva o negativa) entre esta variable y el resto.

##   released colour year age    sex employed citizen checks
## 1      Yes  White 2002  21   Male      Yes     Yes      3
## 2       No  Black 1999  17   Male      Yes     Yes      3
## 3      Yes  White 2000  24   Male      Yes     Yes      3
## 4       No  Black 2000  46   Male      Yes     Yes      1
## 5      Yes  Black 1999  27 Female      Yes     Yes      1
## 6      Yes  Black 1998  16 Female      Yes     Yes      0
## 'data.frame':    5226 obs. of  8 variables:
##  $ released: Factor w/ 2 levels "No","Yes": 2 1 2 1 2 2 2 2 2 2 ...
##  $ colour  : Factor w/ 2 levels "Black","White": 2 1 2 1 1 1 2 2 1 2 ...
##  $ year    : int  2002 1999 2000 2000 1999 1998 1999 1998 2000 2001 ...
##  $ age     : int  21 17 24 46 27 16 40 34 23 30 ...
##  $ sex     : Factor w/ 2 levels "Female","Male": 2 2 2 2 1 1 2 1 2 2 ...
##  $ employed: Factor w/ 2 levels "No","Yes": 2 2 2 2 2 2 1 2 2 2 ...
##  $ citizen : Factor w/ 2 levels "No","Yes": 2 2 2 2 2 2 2 2 2 2 ...
##  $ checks  : int  3 3 3 1 1 0 0 1 4 3 ...
##  released     colour          year           age            sex      
##  No : 892   Black:1288   Min.   :1997   Min.   :12.00   Female: 443  
##  Yes:4334   White:3938   1st Qu.:1998   1st Qu.:18.00   Male  :4783  
##                          Median :2000   Median :21.00                
##                          Mean   :2000   Mean   :23.85                
##                          3rd Qu.:2001   3rd Qu.:27.00                
##                          Max.   :2002   Max.   :66.00                
##  employed   citizen        checks     
##  No :1115   No : 771   Min.   :0.000  
##  Yes:4111   Yes:4455   1st Qu.:0.000  
##                        Median :1.000  
##                        Mean   :1.636  
##                        3rd Qu.:3.000  
##                        Max.   :6.000

Ejercicio 11

Dibuja un histograma (hist) de la variable checks. De esta forma podemos observar si las personas a las que intervino la policía contaban con precedentes.

Ejercicio 12

Ahora dibuja la variable released frente a colour, para conocer si la liberación de las personas arrestadas pudo estar relacionada con su raza. Repite el análisis con las variables sex y employed. ¿Ves alguna relación? Por otra parte, ¿se distribuyeron por igual los arrrestos entre hombres y mujeres, empleados y no empleados? En la función plot puede ponerse un título para el gráfico a través del parámetro main.

Ejercicio 13

A simple vista parece que las mayores diferencias se ven por la variable employed. Obtén una tabla de contingencia entre released y employed, y entre los antecedentes previos checks y employed.

##      
##         No  Yes
##   No   349  543
##   Yes  766 3568
##    
##       No  Yes
##   0  201 1650
##   1  145  709
##   2  172  617
##   3  304  649
##   4  244  399
##   5   46   81
##   6    3    6

Ejercicio 14

Comprueba si existen diferencias en la liberación en función de los antecedentes previos a través del estadístico \(\chi^2\).

## Warning in chisq.test(table(Arrests$released, Arrests$checks)): Chi-squared
## approximation may be incorrect
## 
##  Pearson's Chi-squared test
## 
## data:  table(Arrests$released, Arrests$checks)
## X-squared = 345.21, df = 6, p-value < 2.2e-16

Verás que la función ha lanzado un warning sobre los resultados. En Google puedes encontrar una respuesta a este mensaje: “It gave the warning because many of the expected values will be very small and therefore the approximations of p may not be right”.

Ejercicio 15

Realiza ahora un análisis de regresión que explique las diferencias en la variable released utilizando el resto de variables. ¿Son todos los coeficientes significativos? ¿Contar con antecedentes aumenta o disminuye la probabilidad de terminar el arresto con una liberación? ¿Podríamos decir que las personas con antecedentes tienen más o menos probabilidad de terminar siendo liberados que el resto? Explica cómo influye cada variable en la posibilidad de terminar el arresto con una liberación, y si el comportamiento de las variables coincide con el que habrías predicho.

## 
## Call:
## glm(formula = released ~ ., family = "binomial", data = Arrests)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.3909   0.3579   0.4320   0.6047   1.7067  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  9.371821  56.717803   0.165    0.869    
## colourWhite  0.389109   0.085663   4.542 5.56e-06 ***
## year        -0.004218   0.028379  -0.149    0.882    
## age          0.002236   0.004631   0.483    0.629    
## sexMale      0.007317   0.150189   0.049    0.961    
## employedYes  0.757302   0.084735   8.937  < 2e-16 ***
## citizenYes   0.576519   0.104246   5.530 3.20e-08 ***
## checks      -0.364101   0.025984 -14.013  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 4776.3  on 5225  degrees of freedom
## Residual deviance: 4299.1  on 5218  degrees of freedom
## AIC: 4315.1
## 
## Number of Fisher Scoring iterations: 5

Ejercicio 16

Obtén la matriz de confusión de la anterior regresión, en términos absolutos y relativos, y comenta los resultados obtenidos.

##            Predicted No Predicted Yes Total
## Actual No            55           837   892
## Actual Yes           62          4272  4334
## Total               117          5109  5226
##            Predicted No Predicted Yes    Total
## Actual No    0.01052430     0.1601607 0.170685
## Actual Yes   0.01186376     0.8174512 0.829315
## Total        0.02238806     0.9776119 1.000000

Como puedes observar, el modelo se ha centrado en resolver fundamentalmente el caso en que ha habido liberaciones, ya que es el que representa un mayor porcentaje del total.

Como muchas de las variables del modelo han resultado no significativas (¿qué significa esto?), vamos a repetir el análisis de regresión pero guardando únicamente aquellos coeficientes que hayan resultado significativos. Para ellos podemos hacer uso de la regresión paso a paso (stepwise regression) de la librería MASS:

library(MASS)
step.regresion <- stepAIC(regresion, trace = FALSE)
summary(step.regresion)
## 
## Call:
## glm(formula = released ~ colour + employed + citizen + checks, 
##     family = "binomial", data = Arrests)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.3580   0.3579   0.4316   0.6061   1.6982  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  1.00474    0.12741   7.886 3.12e-15 ***
## colourWhite  0.38915    0.08521   4.567 4.95e-06 ***
## employedYes  0.75367    0.08398   8.974  < 2e-16 ***
## citizenYes   0.56839    0.09916   5.732 9.91e-09 ***
## checks      -0.36283    0.02573 -14.101  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 4776.3  on 5225  degrees of freedom
## Residual deviance: 4299.3  on 5221  degrees of freedom
## AIC: 4309.3
## 
## Number of Fisher Scoring iterations: 5

Ejercicio 17

Repite la matriz de confusión sobre la nueva regresión y comenta si hay diferencias respecto de la primera regresión.

3 Resumen de funciones empleadas en el capítulo

  • ? (help)
  • confusion_matrix
  • curve
  • glm
  • head
  • hist
  • install.packages
  • library
  • mean, sd, var, range
  • plot
  • qchisq, dchisq, chisq.test
  • stepAIC
  • summary
  • table
  • tail

  1. Los tipos de variables en R se estudian con detalle en un capítulo posterior↩︎