---
title: "Regresión lineal simple"
author: "Un primer ejemplo de estadística reproducible con [R markdown](http://rmarkdown.rstudio.com/)"
output: html_document
---

Vamos a estudiar la relación entre 
el peso del corazón y el peso del cuerpo del gato doméstico. 
Disponemos de los datos de la libreria `MASS` de R, de modo que lo 
primero es cargar dicha libreria (si no la has instalado, tendrás que hacerlo antes de cargarla en memoria) y averiguar en qué formato se presentan los datos

```{r, echo =FALSE}
library(MASS)
```

es decir, se trata de un `data.frame` llamado `cats`. Veamos la cabecera de los datos 

```{r, echo=FALSE}
head(cats)
```

Hay tres variables, `sex` (female, male), `Bwt` (iniciales de body weight, peso del cuerpo, en kilogramos) y `Hwt` (heart weilgt, peso del corazón, en gramos), cada una de ellas en una columna 

Empezamos representando las variables `Bwt` (explicativa) frente a `Hwt` (respuesta)

```{r, echo = FALSE, fig.align='center', fig.height=3, fig.width= 4}
x = cats$Bwt
y = cats$Hwt
plot(x, y,  xlab = "Bwt (kg)", ylab = " Hwt (g)", cex = .5)
```

A la vista del diagrama, parece que cuanto mayor es el peso del gato, mayor es el peso su corazón (aunque no es imprescindible, puedes usar colores para distingir hembras (F) de machos (M))

```{r,  echo = FALSE, fig.align='center', fig.height=3, fig.width= 4}
plot(x, y, col=c("red","blue")[cats$Sex], xlab = "Bwt (kg)", ylab = " Hwt (g)", cex = .5)
```

A partir del resumen del modelo de regresión lineal, obtenemos el coeficiente de correlación de Pearson

```{r, echo = FALSE, fig.align='center', fig.height=3, fig.width= 4}
lmXY = lm(y ~ x)
resumen_cats = summary(lmXY)
resumen_cats$r.squared
```

que no es muy bueno. Aún así, vamos a seguir adelante para analizar el resto de información que podemos obtener y los procedimientos para obtenerla. 

#### Representar la recta de regresión  

junto con la nube de dispersión

```{r, echo = FALSE, echo = FALSE, fig.align='center', fig.height=3, fig.width= 4}
plot(x, y, xlab = "Bwt (kg)", ylab = " Hwt (g)", cex = .5)
abline(lmXY)
```

Además, es fundamental comprobar que se cumplen las hipótesis necesarias 
para que funcione la técnica de regresión lineal:  

   * Normalidad en y para cada valor de x
   * Homocedasticidad en los residuos


```{r, echo = FALSE, echo = FALSE, fig.align='center', fig.height=4, fig.width= 7}
par(mfrow = c(1,2)) # gráficos de tabla 1 fila x 2 columnas
## HOMOCEDASTICIDAD
# Obtener el residuo de cada observación
cats.res = resid(lmXY)
# Representar esosresidios
plot(cats$Bwt, cats.res, cex =0.5, main = "Residuos")
abline(0, 0)

# NORMALIDAD
qqnorm(cats.res, main = "QQ-plot")
qqline(cats.res)
par(mfrow = c(1,1)) # gráficos de tabla uno en uno
```

   * **Normalidad**: gráfico izquierda, se observan ligeras desviaciones respecto de la normalidad (falta de simetría respecto del cero en algunas columnas), sobre todo para los individuos más grandes. Gráficode la derecha, el qq-plot muestra un ajuste bastante bueno.
   * **Homocedasticidad**: gráfico izquierda, no se infringe severamente, aunque convendría eliminar algunos outlayers. 


**Calcular el peso del corazón predicho** por el modelo para el peso del cuerpo de 
2, 2.25, 2.5, 2.75 y 3kg:

```{r, echo = F}
predict(lmXY, newdata = data.frame(x=seq(from = 2, to = 3, length.out = 5)))
```



#### Cálculo del intervalo de confianza para los coeficientes de la recta de regresión 

Nos referimos a la recta teórica 
```{r}
confint(lmXY)
```

La pendiente de la recta es significativamente distintade cero.


#### Representar las bandas de confianza

La idea es  

  *  calcular el valor estimado (esperado) de la variable `y` para muchos valores de la variable explicativa 
  y
  * calcular el intervalo de confianza para cada uno de esos valores predichos
  * al representar conjuntamente cada uno de esos intervalos, el resultado son las bandas de confianza: una superior y otra inferior, al unir todos los extremos superiores e inferiores de los intervalos.
  
```{r, echo =FALSE, echo = FALSE, fig.align='center', fig.height=4, fig.width= 5}
plot(x, y)
abline(lmXY)
# lista de 250 valores de x para los que calcularesmos el valor predicho
pred.x = seq(min(x), max(x), length.out=250)

# la siguiente funcion calcula 
# el intervalo de confianza para los valores predichos para x en pred.x con el nc fijado
confd.intrvl = predict(lm(y ~ x), interval="confidence", newdata=data.frame(x = pred.x), level=0.95)

# la siguiente funcion representa los intervalos de confinza y une los extremos con un segmento
matlines(pred.x, confd.intrvl[ ,2:3], lty=c(1, 1, 1), col="red", lwd=2)
```

Muchos puntos escapan a las bandas de confianza. Esto es coherente con un coeficiente de correlación bajo (poca tendencia de los puntos a estar alineados).
