Presentación del problema.

Los datos de este fichero se refieren a una muestra de gatos domésticos. La tabla contiene las variables

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 y lectura de los datos.

gatos = read.table(file = "Practica09-cats.csv", header = TRUE, sep = ";")

Ejercicio 1

Comprueba que la lectura ha sido correcta: visualiza la cabecera de la tabla, obtén un resumen de las variables y calcula las dimensiones de la tabla

Para comprobar que la lectura ha sido correcta vamos a ver las primeras lineas de la tabla

head(gatos)
##   Sex Bwt Hwt
## 1   F 2.0 7.0
## 2   F 2.0 7.4
## 3   F 2.0 9.5
## 4   F 2.1 7.2
## 5   F 2.1 7.3
## 6   F 2.1 7.6

las dimensiones de la tabla y usamos summary para obtener una exploración inicial de las variables:

dim(gatos)
## [1] 144   3
summary(gatos)
##      Sex                 Bwt             Hwt       
##  Length:144         Min.   :2.000   Min.   : 6.30  
##  Class :character   1st Qu.:2.300   1st Qu.: 8.95  
##  Mode  :character   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.50

Para reaprovechar el código de la pantilla, renombramos las variables

x = gatos$Bwt
y = gatos$Hwt

Ejercicio 2

Visualiza la nube de puntos para determinar de forma visual si tiene sentido utilizar un modelo lineal para relacionar las variables Bwt y Hwt. 

plot(x, y, 
     xlab = "Peso cuerpo", ylab = "Peso corazón", 
     main = "Gatos")

Ejercicio 3

Estudia la existencia de datos atípicos calculando con el boxplot de cada una de las variables.

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

bpWwt = boxplot(x, xlab = "Peso cuerpo")
stripchart(x, method = "jitter", add = TRUE, 
           vertical = TRUE, pch=19, 
           col= "red", cex=0.5)

bpHwt = boxplot(y, ylab = "Peso corazón")
stripchart(y, 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.

En este caso, son valores plausibles y vamos a conservarlos.

Ejercicio 4

Estima los coeficientes de la recta de regresión tanto puntualmente como por intervalos de confianza (al 90% de confianza). Recuerda que queremos estudiar la posible relación entre Bwt y Hwt. Empezamos dibujando un diagrama de dispersión.

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(y ~ x))
## 
## Call:
## lm(formula = y ~ x)
## 
## Coefficients:
## (Intercept)            x  
##     -0.3567       4.0341

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

modelo$coefficients
## (Intercept)           x 
##  -0.3566624   4.0340627

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

b0 = modelo$coefficients[1] 
b1 = modelo$coefficients[2]

La recta de regresión de Hwt frente a Bwt es, aproximadamente: \[ Hwt = -0.3567 + 4.034 \cdot Bwt \] Esto quiere decir que por cada kilogramo que aumenta el peso del gato, el peso de su corazón aumenta 4 gramos.

Sin embargo, esto son estimaciones puntuales que cambian con la muestra. Para determinar los valores de los coeficientes de la recta con un nivel de confianza del 90%, hacemos

(ci_coef = confint(modelo, level = 0.9))
##                   5 %      95 %
## (Intercept) -1.502834 0.7895096
## x            3.619716 4.4484094
signif(unname(ci_coef[ 2, ]), 2)
## [1] 3.6 4.4

y concluimos que, con un nivel de confianza del 90%, por cada kilo que aumenta el peso del cuerpo del gato el peso de su corazón aumenta entre 3.6 y 4.4 gramos.

# summary(modelo)
# cor(x, y)^2

Ejercicio 5

Añade la recta de regresión al diagrama de dispersión

Podemos usar el comando abline

plot(x, y, xlab = "Peso cuerpo", ylab = "Peso corazón", 
     main = "Gatos", col="blue", pch=19)
abline(modelo, col="red", lwd=3)

Ejercicio 6

Calcula los valores predichos de la variable Hwt para Bwt = 2.75 y 3.25. Haz tanto la predicción puntual como por intervalo de predicción y confianza al 90%.

  • Para hacer la predicción puntual basta con sustituir en la ecuación de la recta

    round(modelo$coefficients[1] + modelo$coefficients[2] * c(2.75, 3.75), digits = 1)

    Como se obtiene el peso del corazón en gramos, basta con usar un decimal.

    Sabemos que las estimaciones puntuales son endebles, y que lo más fiable es hacer estimaciones por intervalos:

  • Intervalo de predicción: reajustar el código de la plantilla; donde pone newdata = data.frame(x = 30) ponemos newdata = data.frame(x = c(2.75, 3.25)), es decir, cambiamos 30 por los valores que queremos predecir ahora, además de ajustar el nivel de confianza:

    # intervalo de prediccion valor predicho 
(predVal = round(
  predict(modelo, newdata = data.frame(x = c(2.75, 3.25)),    interval = "prediction", level = 0.9), 
  digits = 1))
    ##    fit  lwr  upr
## 1 10.7  8.3 13.1
## 2 12.8 10.3 15.2

    donde, de nuevo, hemos redondeado para eliminar decimales supérfluos. Ahora podemos afirmar que, con una probabilidad de 0.9, el valor esperado para los pesos del cuerpo de 2.75 y 3.25 estará, respectivaente, en los intervalos (8.3, 13.1) y (10.3, 15.2).

  • En cuanto a los intervalos de confianza, en este caso lo son para la media de los valores predichos. Es decir, en intervalo de predicción se refiere a valores individuales (para cierto valor fijo de la variable explicativa) y los intervalos de confianza se refieren al valor medio de dichos valores individuales. Son, por tanto, intervalos más pequeños:

    # intervalo de prediccion valor predicho 
(predVal2 = round(
  predict(modelo, newdata = data.frame(x = c(2.75, 3.25)),    interval = "confidence", level = 0.9), 
  digits = 1))
    ##    fit  lwr  upr
## 1 10.7 10.5 10.9
## 2 12.8 12.5 13.1

    Ahora podemos afirmar que, con una probabilidad de 0.9, el valor promedio esperado para los pesos del cuerpo de 2.75 y 3.25 estará, respectivaente, en los intervalos (10.5, 10.9) y (12.5, 13.1).

intervalo de confianza valor promedio predicho

predict(lmXY, newdata = data.frame(x = 30), interval = “confidence”, level = 0.95)

Ejercicio 7

Calcula e interpreta el coeficiente de correlación. A continuación calculamos el coeficiente de correlación \(r\) de este modelo.

cor(x, y)
## [1] 0.8041274

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.

Ejercicio 8

Calcula e interpreta el coeficiente de determinación.

Vale

cor(x, y)^2
## [1] 0.6466209

e indica qué porcentaje de la variabilidad del peso del corazón de los gatos se explica a través del peso de su cuerpo. El resto (de la variabilidad) se explica con otros factores no considerados en el modelo. Podrían añadirse variables, es lo que se conoce como regresión lineal múltiple, una herramienta algo más avanzada.