Introducción.

En esta práctica vamos a trabajar con datos procedentes del libro

Environmental and Ecological Statistics with R, Second Edition

Concretamente vamos a utilizar este fichero csv

https://raw.githubusercontent.com/songsqian/eesR/master/R/Data/laketrout2.csv

Este fichero de datos contiene observaciones sobre el contenido en PCBs (bifenilos policlorados) de truchas pescadas en el Lago Michigan. Se sabe, desde hace décadas, que el consumo de pescado que contiene niveles altos de PCBs puede resultar perjudicial, especialmente en el caso de niños y mujeres embarazadas. Una de las preocupaciones de los investigadores al analizar ese problema era establecer una posible relación entre el tamaño de la trucha (variable length en la tabla) y el contenido en PCBs (variable pcb). En los siguientes ejercicios vamos a explorar esa relación, pero usaremos además algunas otras variables de la tabla.

Ejercicio 1

Descarga de los datos.

  • Selecciona la carpeta de trabajo adecuada para esta sesión con R.
  • Descarga el fichero y guárdalo en la subcarpeta datos de tu carpeta de trabajo**.
  • Explora el fichero con un editor de texto como el Bloc de Notas.
  • Lee ese fichero con R, usando el comando read.table y guárdalo en un data.frame con el nombre laketrout.
  • Comprueba que la lectura ha sido correcta.

Ejercicio 2

Exploración inicial de los datos.

  • ¿Cuántes filas y columnas tiene la tabla?

  • Utiliza summary para ver información básica sobre las variables de esa tabla. ¿Hay datos ausentes en alguna variable?

summary(laketrout)
##      length           pcb               n             lgpcb        
##  Min.   : 0.00   Min.   : 0.090   Min.   :1.000   Min.   :-2.4079  
##  1st Qu.:22.00   1st Qu.: 1.200   1st Qu.:1.000   1st Qu.: 0.1823  
##  Median :25.20   Median : 2.300   Median :1.000   Median : 0.8329  
##  Mean   :24.61   Mean   : 4.132   Mean   :1.019   Mean   : 0.8942  
##  3rd Qu.:28.00   3rd Qu.: 4.320   3rd Qu.:1.000   3rd Qu.: 1.4632  
##  Max.   :34.64   Max.   :43.800   Max.   :4.000   Max.   : 3.7796  
##  NA's   :15                                                        
##       year            y       
##  Min.   :1974   Min.   : 0.0  
##  1st Qu.:1983   1st Qu.: 9.0  
##  Median :1985   Median :11.0  
##  Mean   :1985   Mean   :11.4  
##  3rd Qu.:1989   3rd Qu.:15.0  
##  Max.   :2000   Max.   :26.0  
##
  • Utiliza complete.cases (ya apareció en la anterior práctica) para comprobar si hay filas con datos ausentes en esa tabla. Si las hay, fabrica una nueva tabla que no contenga esas filas y llámala igual que la tabla original (laketrout).
all(complete.cases(laketrout))
## [1] FALSE

Es decir, hay filas incompletas. Pasamos a eliminarlas

laketrout = laketrout[complete.cases(laketrout), ]
# comprobar que ya no hay filas incompletas
all(complete.cases(laketrout))
## [1] TRUE
  • Dibuja un diagrama de caja (boxplot) de la variable length. ¿Hay valores atípicos? ¿Cuánto valen? ¿Qué crees que debes hacer con esos valores (es posible que la respuesta depende del valor)?
boxplot(laketrout$length)

Parece que hay una trucha de 0 cm (o casi), que no interesa en el estudio (ni nos creemos)

(minLength = which.min(laketrout$length))
## [1] 130
laketrout[minLength, ]
##     length pcb n    lgpcb year y
## 130      0  12 1 2.484907 1982 8
laketrout = laketrout[-minLength, ]
  • Haz lo mismo con la variable pcb, dibujando un boxplot para empezar a analizar la variable.
boxplot(laketrout$pcb)

Si se eliminan las truchas con mayor contenido de pcb (valores atípicos) estamos introduciendo un sesgo en el estudio, porque eliminamos los ejemplares más contaminados.

  • Después, cuando hayas decidido qué hacer con los atípicos, dibuja histogramas y diagramas de densidad de ambas variables, length y pcb. ¿Las distribuciones son simétricas?
par(mfrow = c(1,2)) # gráficos en 1 fila y 2 columnas
plot(density(laketrout$length))
plot(density(laketrout$pcb))

par(mfrow = c(1,1)) # gráficos en 1 fila y 1 columna

Ejercicio 3

Modelo de regresión lineal.

  • Vamos a tratar de explorar la posible relación entre length y pcb. Empieza por dibujar un diagrama de dispersión con length en el eje \(x\) (variable explicativa) y pcb en el eje \(y\) (variable respuesta).
plot(laketrout$length, laketrout$pcb)

  • ¿Crees que estos datos se describen bien mediante una recta?

  • Antes hemos visto que los datos de pcb son asimétricos (a la derecha). Este tipo de asimetría se presenta muy a menudo en datos biológicos. A menudo, en ese tipo de situaciones, recurrimos a transformar los datos a otra escala, como la escala logarítmica, en la que la relación entre las dos variables puede ser más fácil de analizar. El conjunto de datos que estamos usando ya contempla esa transformación: la columna lgpcb contiene los logaritmos (neperianos) de los valores pcb. Dibuja un nuevo un diagrama de dispersión con length en el eje \(x\) (variable explicativa) pero ahora usando lgpcb en el eje \(y\) (variable respuesta).

head(log(laketrout$pcb))
## [1] 3.443618 2.066863 3.284664 2.116256 2.424803 2.533697
head(laketrout$lgpcb)
## [1] 3.443618 2.066863 3.284664 2.116256 2.424803 2.533697
plot(laketrout$length, laketrout$lgpcb)

  • ¿Crees que estos datos transformados se representan mejor con una recta que los datos sin transformar?

  • Efecticamente, parecen más alineados. Hay que detenerse un momento para entender la transformación:
  • Los datos parecen ajustarse a una curva exponencial \[pcb = a*e^{length} \]

  • Al calcular el logaritmo neperiano del PCB lo que se hace es calcular el logratirmo de ambos miembros de la igualdad \[\ln(pcb) = \ln(b_0*e^{b_1*length})=\ln(b_0)+\ln(e^{b_1*length})= \ln(b_0)+b_1*length\] es decir, si llamamos lgpcb a \(\ln(pcb)\) resulta que tenemos una relación lineal \[lgpcb = \ln(b_0)+b_1*length\] a la que podemos aplicar la técnica de regresión lineal.

  • Calcula ahora dos modelos de regresión lineal. Uno para los datos no transformados y otro para los transformados.

(modelo1 = lm(pcb ~length, data = laketrout))
## 
## Call:
## lm(formula = pcb ~ length, data = laketrout)
## 
## Coefficients:
## (Intercept)       length  
##     -8.0492       0.4962
(modelo2 = lm(lgpcb ~length, data = laketrout))
## 
## Call:
## lm(formula = lgpcb ~ length, data = laketrout)
## 
## Coefficients:
## (Intercept)       length  
##     -2.1597       0.1243
  • Vuelve a dibujar los dos diagramas de dispersión, añadiendo la recta de regresión correspondiente a cada uno de ellos.
par(mfrow = c(1, 2))
plot(laketrout$length, laketrout$pcb)
abline(modelo1, col="red", lwd= 2)

plot(laketrout$length, laketrout$lgpcb)
abline(modelo2, col="red", lwd= 2)

par(mfrow = c(1, 1))
  • Calcula los coeficientes (pendiente, ordenada en el origen) de la recta de regresión en ambos modelos.

Ya están calculados

  • Calcula también los coeficientes de correlación \(r\) de ambos modelos. ¿Qué modelo crees que es más adecuado?
(r = cor(laketrout$pcb, laketrout$length))
## [1] 0.3943744
(r_ln = cor(laketrout$lgpcb, laketrout$length))
## [1] 0.5969693

Es preferible el modelo2 al modelo1.

  • Utilizando ese modelo, ¿cuál es tu estimacion del contenido en pcb de una trucha con length=28? ¿Y con length = 11? Como
range(laketrout$length)
## [1] 11.57 34.64

resulta que el valor length = 18 está dentro del recorrido de la variable explicativa, por lo que tiene sentido hacer una estimación

# pcb predicho para trucha de 28cm
(predichoLn = modelo2$coefficients[1] + modelo2$coefficients[2] *28)
## (Intercept) 
##    1.321893

pero ese valor predichoLn es en realidad el logaritmo del valor que buscamos, porque la recta de regresión la hemos calculado después de transformar la variable pcb en su logaritmo lgpcb. Por eso falta calcular su exponencial

exp(predichoLn)
## (Intercept) 
##    3.750514

para obtener el valor predicho en la escala adecuada. Para length = 11 estamos fuera del recorrido de la variable explicativa, y no es admisible usar el modelo de regresión lineal para hacer esa estimación.

Ejercicio 4

  • Vamos a tratar de estudiar cómo ha cambiado el contenido de pcb a lo largo de los años. Para ello disponemos de la variable year, que indica en que año se hizo cada una de las observaciones. Haz una tabla de frecuencias de esa variable.
table(laketrout$year)
## 
## 1974 1975 1976 1977 1978 1979 1980 1982 1983 1984 1985 1986 1987 1989 1990 
##   28   29   27    3   26    3    7    9   47   86  142   25   19   20   29 
## 1991 1992 1994 1996 1997 1998 2000 
##   10   29   33   22   10   19    9
  • Para simplificar el análisis, vamos a agrupar la variable year en intervalos, que correspondan con las tres décadas que abarca el estudio. Usa para esto la función cut, con cortes en los años 1971, 1981, 1991 y 2001. El resultado es un factor que podemos llamar decada. ¿Cuántos niveles tiene? Pídele a tu profesor que te enseñe a usar labels para elegir los nombres de esos niveles y que te enseñe también a añadir a la tabla una nueva columna con los valores del factor decada.
laketrout$decada = cut(laketrout$year, breaks = c(1971, 1981, 1991, 2001), labels = c("70", "80", "90"))
  • Calcula el valor medio de pcb según la década de la observación. Dibuja los boxplots de pcb correspondientes a cada nivel del factor. ¿Qué observas?
aggregate(pcb ~ decada, data= laketrout, FUN = mean)
##   decada      pcb
## 1     70 9.103659
## 2     80 3.272680
## 3     90 2.104992
boxplot(pcb ~ decada, data= laketrout)