Práctica 3. Primera parte de la sesión.
Universidad de Alcalá. Curso 2017-18.
Estadística (650008). Grado en biología sanitaria.
Actualizado: 2017-10-05
Introducción.
Este fichero es un documento reproducible en formato Rmarkdown. Puedes abrirlo en RStudio y procesarlo usando el botón Knit para obtener un documento html (que puedes ver en tu navegador).
Presentación del problema.
En esta práctica vamos a empezar trabajando con los datos de este fichero:
http://www3.uah.es/marcos_marva/sanitaria1718/datos/cats.csv
Los datos del fichero proceden de uno de los conjuntos de datos de ejemplo que incluye R. Puedes leer más sobre este conjunto de datos en este enlace:
https://stat.ethz.ch/R-manual/R-devel/library/MASS/html/cats.html
Vamos a utilizar esos datos para explorar el posible grado de asociación entre dos variables relacionadas con la anatomía de los gatos: su peso corporal y el peso del corazón del gato.
Descarga de los datos.
Recuerda antes de seguir que debes seleccionar la carpeta de trabajo adecuada para esta sesión con R.
Descargamos el fichero y lo guardamos en la subcarpeta
datosde tu carpeta de trabajo**.Exploramos el fichero con un editor de texto como el Bloc de Notas, para ver su estructura.
Leemos ese fichero con R, usando el comando
read.table.
cats = read.table(file = "datos/cats.csv", header = TRUE, sep = ";")Para comprobar que la lectura ha sido correcta vamos a ver las dimensiones de la tabla y usamos summary para obtener una exploración inicial de las variables:
dim(cats)## [1] 144 3summary(cats)## Sex Bwt Hwt
## F:47 Min. :2.000 Min. : 6.30
## M:97 1st Qu.:2.300 1st Qu.: 8.95
## Median :2.700 Median :10.10
## Mean :2.724 Mean :10.63
## 3rd Qu.:3.025 3rd Qu.:12.12
## Max. :3.900 Max. :20.50Utilizamos complete.cases y all para comprobar si hay filas con datos ausentes en esa tabla.
Dibujamos diagramas de caja (boxplot) de las variables numéricas Bwt y Hwt, con los valores superpuestos.
bpWwt = boxplot(cats$Bwt)
stripchart(cats$Bwt,
method = "jitter", add = TRUE, vertical = TRUE,
pch=19, col= "red", cex=0.5)
bpHwt = boxplot(cats$Hwt)
stripchart(cats$Hwt,
method = "jitter", add = TRUE, vertical = TRUE,
pch=19, col= "red", cex=0.5)
Como ves en este último diagrama, hay valores atípicos. El tratamiento de esos valores atípicos es distinto en cada caso y no se pueden establecer reglas generales. A veces un valor atípico se debe simplemente a un error al anotar los datos y, en tal caso, lo más sensato es eliminar esa observación antes de seguir adelante con el análisis. Pero en otras ocasiones un valor atípico apunta a un fenómeno interesante y eliminarlo sería un error.
Pero, como hemos dicho, a veces la decisión que tomaremos será la de eliminar esos valores. Sin entrar a discutir si es lo más correcto o no, vamos a usar estos datos para ver como podemos eliminar los valores atípicos de Hwt. El primer paso es identificarlos.
(atipicos = bpHwt$out)## [1] 17.2 20.5Localizamos las posiciones que ocupan usando %in%:
(filasAtipicos = which(cats$Hwt %in% atipicos))## [1] 135 144Ahora seleccionamos todas las filas de la tabla menos las que contienen atípicos. Esto se hace así:
catsSinAtipicos = cats[-filasAtipicos, ]Fíjate en que hemos vuelto a ponerle el nombre cats a la tabla que se obtiene eliminando esos dos valores atípicos de Hwt.
Hemos hecho esto para eliminar los valores atípicos de Hwt. Pero si vuelves a dibujar el boxplot de esa variable:
bpHwt = boxplot(catsSinAtipicos$Hwt)
verás que lo único que hemos conseguido es convertir en atípico el siguiente valor más grande de la variable. Insistimos: eliminar los atípicos por principio es, a menudo, un error. Y como muestra este ejemplo, a veces ni siquiera sirve de gran cosa.
En el resto del análisis vamos a mantener la tabla cats original.
Modelo de regresión lineal.
Recuerda que queremos estudiar la posible relación entre Bwt y Hwt. Empezamos dibujando un diagrama de dispersión.
plot(cats$Bwt, cats$Hwt, xlab="Bwt", ylab="Hwt", col="blue", pch=19)
¿Hay algún detalle que te llame la atención en este gráfico?
El diagrama muestra indicios de una correlación entre esas dos variables. Vamos a construir el modelo de regresión lineal correspondiente usando la función lm .
(modelo = lm(Hwt ~Bwt, cats))##
## Call:
## lm(formula = Hwt ~ Bwt, data = cats)
##
## Coefficients:
## (Intercept) Bwt
## -0.3567 4.0341y extraemos los coeficientes de la recta de regresión con:
modelo$coefficients## (Intercept) Bwt
## -0.3566624 4.0340627Como ves, la pendiente de la recta es (aproximadamente) 4.034, mientras que la ordenada en el origen es -0.3567.
Recuerda que hemos escrito la ecuación genérica de la recta de regresión así: \[
y = b_0 + b_1 \cdot x
\] siendo \(x\) la variable predictora (en nuestro caso Bwt) e \(y\) la variable respuesta (en nuestro caso Hwt). Como se obtiene:
(b = modelo$coefficients)## (Intercept) Bwt
## -0.3566624 4.0340627b0 = b[1] # cuidado con confundir los nombres teóricos b0, b1 con los índices de R
b1 = b[2]La recta de regresión de Hwt frente a Bwt es: \[
Hwt = -0.357 + 4.034 \cdot Bwt
\]
Puedes obtener mucha más información sobre el modelo usando summary. No te preocupes si en este momento no entiendes mucha de esta información, quedará más claro más adelante.
summary(modelo)##
## Call:
## lm(formula = Hwt ~ Bwt, data = cats)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.5694 -0.9634 -0.0921 1.0426 5.1238
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.3567 0.6923 -0.515 0.607
## Bwt 4.0341 0.2503 16.119 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.452 on 142 degrees of freedom
## Multiple R-squared: 0.6466, Adjusted R-squared: 0.6441
## F-statistic: 259.8 on 1 and 142 DF, p-value: < 2.2e-16Vamos a añadir el dibujo de la recta de regresión al diagrama de dispersión usando abline
plot(cats$Bwt, cats$Hwt, xlab="Bwt", ylab="Hwt", col="blue", pch=19)
abline(modelo, col="red", lwd=3)
A continuación calculamos el coeficiente de correlación \(r\) de este modelo.
cor(cats$Bwt, cats$Hwt)## [1] 0.8041274Teniendo en cuenta ese valor, el tipo de datos de que se trata y a la vista del gráfico de dispersión, podemos considerar que el ajuste de la recta de regresión a los datos es razonablemente bueno.
Predicción usando el modelo.
Vamos a utilizar ese modelo para estimar o predecir el valor de Hwt en un gato con un valor de Bwt (peso corporal) igual a 2.55. Para eso basta con sustituir en la ecuación de la recta:
peso = 2.55
(valorHwt = b0 + b1 * peso)## (Intercept)
## 9.930197El valor que hemos usado no aparece en la muestra, pero está dentro del recorrido de la variable explicativa:
range(cats$Bwt)## [1] 2.0 3.9Recuerda que debes evitar la extrapolación, que se produce al usar la recta de regresión para predecir a partir de un valor fuera de este recorrido.
El modelo que obtenemos con lm permite también obtener los valores predichos por la recta para los 144 puntos de la muestra, usando:
modelo$fitted.values## 1 2 3 4 5 6 7
## 7.711463 7.711463 7.711463 8.114869 8.114869 8.114869 8.114869
## 8 9 10 11 12 13 14
## 8.114869 8.114869 8.114869 8.114869 8.114869 8.518276 8.518276
## 15 16 17 18 19 20 21
## 8.518276 8.518276 8.518276 8.518276 8.921682 8.921682 8.921682
## 22 23 24 25 26 27 28
## 8.921682 8.921682 8.921682 8.921682 8.921682 8.921682 8.921682
## 29 30 31 32 33 34 35
## 8.921682 8.921682 9.325088 9.325088 9.325088 9.325088 9.728494
## 36 37 38 39 40 41 42
## 9.728494 10.131901 10.131901 10.131901 10.535307 10.535307 10.535307
## 43 44 45 46 47 48 49
## 11.342119 11.342119 11.342119 11.745526 11.745526 7.711463 7.711463
## 50 51 52 53 54 55 56
## 8.114869 8.518276 8.518276 8.518276 8.518276 8.518276 8.518276
## 57 58 59 60 61 62 63
## 8.518276 8.518276 8.921682 9.325088 9.325088 9.325088 9.325088
## 64 65 66 67 68 69 70
## 9.325088 9.728494 9.728494 9.728494 9.728494 9.728494 9.728494
## 71 72 73 74 75 76 77
## 9.728494 9.728494 10.131901 10.131901 10.131901 10.131901 10.131901
## 78 79 80 81 82 83 84
## 10.131901 10.535307 10.535307 10.535307 10.535307 10.535307 10.535307
## 85 86 87 88 89 90 91
## 10.535307 10.535307 10.535307 10.938713 10.938713 10.938713 10.938713
## 92 93 94 95 96 97 98
## 10.938713 10.938713 10.938713 11.342119 11.342119 11.342119 11.342119
## 99 100 101 102 103 104 105
## 11.342119 11.745526 11.745526 11.745526 11.745526 11.745526 11.745526
## 106 107 108 109 110 111 112
## 11.745526 11.745526 11.745526 12.148932 12.148932 12.148932 12.148932
## 113 114 115 116 117 118 119
## 12.148932 12.148932 12.552338 12.552338 12.552338 12.552338 12.552338
## 120 121 122 123 124 125 126
## 12.552338 12.955744 12.955744 12.955744 12.955744 12.955744 13.359151
## 127 128 129 130 131 132 133
## 13.359151 13.359151 13.359151 13.359151 13.762557 13.762557 13.762557
## 134 135 136 137 138 139 140
## 13.762557 13.762557 14.165963 14.165963 14.165963 14.165963 14.569370
## 141 142 143 144
## 14.972776 14.972776 15.376182 15.376182Estos son los valores que hemos llamado \({\hat y}_i\) en el libro. Para que quede más claro: el primer valor de la muestra tiene este valor de Bwt:
cats$Bwt[1]## [1] 2Y si sustituyes esto en la recta \(y = b_0 + b_1 x\):
b0 + b1 * cats$Bwt[1]## (Intercept)
## 7.711463puedes comprobar que este es el primer valor de los que hemos obtenido con modelo$fitted.values.
Regresión para un subconjunto de los datos.
Vamos a terminar proponiéndote que repitas el análisis anterior, pero sólo para las gatas. Tendrás que empezar seleccionando esas filas de la tabla, con algo como:
gatas = [??, ??]Empieza completando este comando. Luego dibuja el diagrama de dispersión y usa usa lm para analizar la relación entre Bwt y Hwt en el caso de las gatas. ¿Obtienes valores muy distintos de los del conjunto completo de observaciones?