---
title: "Práctica 3. Primera parte de la sesión."
author: "Estadística (650008). Grado en biología sanitaria."
date: 'Actualizado: `r format(Sys.time(), "%Y-%m-%d")`'
output:
  html_document: default
subtitle: Universidad de Alcalá. Curso 2017-18.
---

# 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](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](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 `datos` de 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`.

```{r}
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:

```{r}
dim(cats)
summary(cats)
```


Utilizamos `complete.cases` y `all` para comprobar si hay filas con datos ausentes en esa tabla. 

```{r echo=FALSE, eval=FALSE}
all(complete.cases(cats))
```

Dibujamos diagramas de caja (boxplot) de las variables numéricas `Bwt` y `Hwt`, con los valores superpuestos. 

```{r}
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.

```{r}
(atipicos = bpHwt$out)
```
Localizamos las posiciones que ocupan usando `%in%`:

```{r}
(filasAtipicos = which(cats$Hwt %in% atipicos))
```

Ahora seleccionamos todas las filas de la tabla **menos las que contienen atípicos**. Esto se hace así:

```{r}
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: 

```{r}
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.


```{r}
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` .

```{r}
(modelo = lm(Hwt ~Bwt, cats))
```

y extraemos los coeficientes de la recta de regresión con:

```{r}
modelo$coefficients
```

Como ves, la pendiente de la recta es (aproximadamente) `r signif(modelo$coefficients[2], 4)`, mientras que la ordenada en el origen es `r signif(modelo$coefficients[1], 4)`. 

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:

```{r}
(b = modelo$coefficients)
b0 = 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.

```{r}
summary(modelo)
```


Vamos a añadir el dibujo de la recta de regresión al diagrama de dispersión usando `abline`

```{r}
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. 

```{r}
cor(cats$Bwt, cats$Hwt)
```
Teniendo 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:

```{r}
peso = 2.55
(valorHwt = b0 + b1 * peso)
```

El valor que hemos usado no aparece en la muestra, pero está dentro del recorrido de la variable explicativa:


```{r}
range(cats$Bwt)
```
Recuerda 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 `r nrow(cats)` puntos de la muestra, usando: 

```{r}
modelo$fitted.values
```

Estos 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`:

```{r}
cats$Bwt[1]
```

Y si sustituyes esto en la recta $y = b_0 + b_1 x$: 

```{r}
b0 + b1 * cats$Bwt[1]
```

puedes 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:

```{r eval=FALSE, comment=NA}
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?

