Pesca y estrategias de Muestreo y Análisis. Sesión 2
Master en Gestión Sostenible de Recursos Pesqueros
Author
Angelo Santana
Published
October 6, 2022
Mostrar código R
## Librerías necesarias para el análisis de los datos## --------------------------------------------------library(tidyverse)library(rstatix)library(qqplotr)library(flextable)library(DescTools)# Algunas funciones de utilidad# ------------------------------## Media y desviación típica:meanSd <-function(x) sprintf("%.2f (%.2f)", mean(x,na.rm=TRUE),sd(x,na.rm=TRUE))## Mediana y cuartiles:medianIQR <-function(x){ q <-quantile(x,probs=c(0.5,0.25,0.75),na.rm=TRUE)sprintf("%.2f (%.2f; %.2f)",q[1],q[2],q[3])}## Formateo del p-valor para que muestre "<0.0001" cuando p es menor que dicho valor:formatPval=function(p) ifelse(p<0.0001,"<0.0001",sprintf("%.4f",p))
En esta sesión utilizaremos los datos del archivo sargos.csv. A continuación se muestran las primeras filas de este archivo:
Estos datos corresponden a una muestra de peces capturados en el entorno de las distintas islas del archipiélago canario. Para cada ejemplar se ha registrado el sexo y varias medidas biométricas (longitud, longitud hasta la aleta dorsal, etc, peso total, peso de las gónadas y peso del hígado). La variable ptdo vale 1 si el animal tiene parásitos y 0 en caso contrario. Cuando el animal tiene parásitos se cuenta el número de larvas localizadas, que se anotan en la variable larvas.
El modelo de Regresión lineal
El modelo de regresión lineal (simple) es de la forma:
\[y =\beta_0 + \beta_1 x + \varepsilon\]
Este modelo especifica que la variable respuesta\(y\) depende linealmente de una variable explicativa\(x\), siendo \(\beta_0\) la ordenada y \(\beta_1\) la pendiente de la recta que las relaciona. El término \(\varepsilon\) hace referencia a que los valores de \(y\)no pueden predecirse exactamente a partir de los valores de la \(x\). Los valores de \(\varepsilon\) se conocen como residuos del modelo; a su vez, su varianza \(\sigma^2_\varepsilon\) se conoce como varianza residual.
Nótese que la ecuación anterior especifica la relación entre \(x\) e \(y\) en la población donde se miden ambas variables. En las aplicaciones prácticas, normalmente se dispone de una muestra de observaciones \((x_i,y_i)\) de estas variables. El procedimiento lm de R permite obtener la recta ajustada a dicha muestra, de tal forma que:
\[y_i=b_0 +b_1 x_i + e_i\]
El valor de \(b_0\) obtenido en la muestra es una estimación (un valor próximo) del \(\beta_0\) en la población ; a su vez, \(b_1\) es una estimación de \(\beta_1\). Para hacer inferencia en regresión (extender los resultados de la muestra a la población) se asume que los residuos en la población siguen una distribución normal de media 0 y varianza \(\sigma^2_\varepsilon\) constante a lo largo de toda la recta (homoscedasticidad). La varianza de los residuos de la muestra, \(s^2_e\), es una estimación de la varianza residual \(\sigma^2_\varepsilon\).
Debe tenerse en cuenta que los valores estimados dependen de la muestra (muestras distintas darán lugar a diferentes valores estimados). Por ello es conveniente acompañar los valores estimados de sus errores estándar, o intervalos de confianza, para poder tener una idea de la precisión de la estimación. Dichos intervalos se calculan bajo la hipótesis de normalidad y homoscedasticidad de los residuos de la recta. Asimismo suele ser de interés decidir a partir de la muestra si hay evidencia suficiente para asegurar que \(\beta_1\) es distinto de cero, pues si fuese 0, el modelo quedaría reducido a \(y=\beta_0 +\varepsilon\), o lo que es lo mismo, la \(y\) no depende de \(x\) (al menos de forma lineal).
Veremos a continuación un ejemplo para entender bien estas ideas.
La figura siguiente corresponde a la nube de puntos en la que se han representado la longitud (eje \(x\)) y el peso (eje \(y\)) de los sargos de la muestra. La representación gráfica muestra que la relación entre estas variables no es lineal:
Utilizando las propiedades de los logaritmos, a partir de la ecuación anterior se obtiene:
\[e^{log(y)} = e^{b_0 + b_1\cdot log(x)}\]
\[y = e^{b_0}x^{b_1}\]
y llamando \(a=e^{b_0}\) concluimos que el modelo logarítmico es equivalente a:
\[y= a x^{b_1}\]
Esta ecuación es de uso común en el estudio de la alometría (relación entre las distintas dimensiones) del cuerpo de los seres vivos. Cuando, como en nuestro ejemplo, la variable \(y\) es el peso y la \(x\) la talla, puede entenderse que \(x^{b_1}\) viene a ser una medida del volumen del cuerpo del individuo, y por tanto la relación anterior expresa que el peso del cuerpo es proporcional a su volumen (que es de hecho algo que cabe esperar).
Cuando utilizamos R para estimar la ecuación de esta recta obtenemos como resultado los coeficientes siguientes:
El primer p-valor es el contraste para \(\beta_0\). Como es menor que el nivel de significación habitual \(\alpha=0.05\), ello quiere decir que la muestra contiene evidencia suficiente para asegurar que \(\beta_0\) (la ordenada en la población) es distinta de 0. Asimismo, el segundo p-valor corresponde al contraste para \(\beta_1\). Como también es menor que 0.05, podemos concluir que \(\beta_1\) es distinto de cero en la población. Cuando se puede concluir que un coeficiente es distinto de cero, se dice que dicho coeficiente es significativo. En este caso ambos coeficientes han resultado significativos.
El resumen anterior nos muestra también el valor de \(R^2=0.9765\). Ello quiere decir que el ajuste es muy bueno; concretamente, que el 97.65% de la variabilidad en la variable \(y\) (el peso) puede explicarse por la variabilidad en la variable \(x\) (el peso); dicho de otra forma, las diferencias en peso se explican casi exclusivamente por las diferencias en talla.
En la tabla anterior figuran también los errores estándar de los coeficientes de la regresión (0.0870 y 0.0285 respectivamente). Cuando se cumplen las hipótesis para la inferencia en regresión (normalidad y homoscedasticidad de los residuos), puede probarse que con un 95% de confianza los valores en la población de \(\beta_0\) y \(\beta_1\) están a una distancia menor de dos errores estándar de sus valores estimados. La función confint de R muestra dichos intervalos con más precisión:
Así pues, el valor estimado de \(\beta_0\) es -2.87, aunque con un 95% de confianza su valor real en la población puede estar entre -3.047 y -2.704. A su vez, el valor estimado de \(\beta_1\) es 2.587, aunque con un 95% de confianza su valor real en la población puede estar entre 2.531 y 2.643.
Los contrastes de hipótesis e intervalos de confianza anteriores solo son válidos cuando se cumplen las hipótesis la normalidad y homoscedasticidad de los residuos. El qqplot de los residuos indica que su distribución es normal:
El test de Shapiro-Wilk no detecta falta de normalidad:
Mostrar código R
shapiro.test(resids)
Shapiro-Wilk normality test
data: resids
W = 0.98828, p-value = 0.0992
Si representamos los residuos frente a las predicciones, no se aprecia que haya problemas de heteroscedasticidad, ya que la dispersión de los puntos alrededor de la recta es homogénea en todo su recorrido:
La función ncvTest del paquete car permite llevar a cabo el test de Breusch-Pagan que contrasta si los residuos son homoscedaticos (hipotesis nula). En este caso se obtiene:
Mostrar código R
# ver https://fhernanb.github.io/libro_regresion/homo.htmlcar::ncvTest(recta)
Non-constant Variance Score Test
Variance formula: ~ fitted.values
Chisquare = 0.09117028, Df = 1, p = 0.76269
E p-valor 0.7627 es mayor que 0.05 y por tanto puede aceptarse la homoscedasticidad de los residuos de la recta.
Análisis de la covarianza
El análisis de la covarianza es el nombre que recibe el modelo lineal cuando incluye como variables explicativas no solo variables continuas, sino también variables binarias 0-1.
A modo de ejemplo consideremos que queremos decidir si la recta de regresión ajustada a los machos de la muestra anterior coincide con la ajustada a las hembras. Gráficamente, cuando separamos ambos sexos obtenemos las nubes de puntos que se muestran a continuación:
Call:
lm(formula = log(peso) ~ log(long), data = subset(sargos, sexo ==
"Hembra"))
Residuals:
Min 1Q Median 3Q Max
-0.140331 -0.054127 -0.003681 0.045802 0.172192
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -2.84858 0.11140 -25.57 <2e-16 ***
log(long) 2.58016 0.03681 70.09 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.07182 on 111 degrees of freedom
Multiple R-squared: 0.9779, Adjusted R-squared: 0.9777
F-statistic: 4912 on 1 and 111 DF, p-value: < 2.2e-16
:::
Como vemos, los coeficientes de ambas ecuaciones son distintos, aunque muy parecidos. La ordenada es -2.93 para los machos y -2.84 para las hembras. La pendiente es 2.60 para los machos y 2.58 para las hembras. Dado que estos valores se han estimado a partir de esta muestra, si se toma una nueva muestra de sargos seguramente estos valores cambiarán. Por tanto se plantea la pregunta: ¿son tan parecidos los valores para machos y hembras en las muestras que podemos aceptar la hipótesis de que en la población ambas rectas son iguales?; o por el contrario la diferencia observada aunque pequeña ¿es suficientemente grande para asegurar que en la población los machos y las hembras se ajustan a rectas distintas?
Dado que las dos rectas se han estimado por separado, no es posible realizar un test estadístico para decidir si la diferencia observada entre las rectas es significativa (son distintas en la población) o no. Para ello ambas rectas deben estimarse en un único modelo y es aquí donde entra en juego el análisis de la covarianza (ANCOVA).
Para construir el modelo ANCOVA consideramos la variable indicatriz siguiente:
$$I_{M_{i}}=
\[\begin{cases}
0 & \textrm{Si el sujeto i es Hembra}\\
1 & \textrm{Si el sujeto i es Macho}
\end{cases}\]
$$
Esta variable vale 1 para los machos y 0 para las hembras. En R puede codificarse fácilmente mediante:
Así, pues, \(b_0+b_{0M}\) es la ordenada de la recta de los machos y \(b_1+b_{1M}\) es la pendiente de dicha recta. Comparando con la recta de las hembras, podemos concluir que \(b_{0M}\) es la diferencia en ordenadas entre ambas rectas, y \(b_{1M}\) es la diferencia en pendientes.
En R podríamos ajustar este modelo de la forma siguiente:
Call:
lm(formula = log(peso) ~ IM + log(long) + IM * log(long), data = sargos)
Coefficients:
(Intercept) IM log(long) IM:log(long)
-2.84858 -0.08313 2.58016 0.02396
Por tanto:
\(b_0=-2.8485\) es la ordenada de la recta de las hembras.
\(b_{0M}=-0.08313\) es la diferencia entre la ordenada de la recta de las hembras y la de los machos.
\(b_1=2.58016\) es la pendiente de la recta de las hembras.
\(b_{1M}=0.02396\) es la diferencia entre la pendiente de la recta de las hembras y la de los machos.
Si calculamos \(b_0+b_{0M}=-2.8485-0.08313=-2.93163\) obtenemos la ordenada de la recta de los machos, que coincide con la que ya habíamos obtenido más arriba.
Asimismo, si calculamos \(b_1+b_{1M}=2.58016+0.02396=2.60412\) obtenemos la pendiente de la recta de los machos, que también coincide con la que ya habíamos obtenido más arriba.
El modelo del análisis de la varianza nos permite contrastar si los coeficientes \(b_{0M}\) y \(b_{1M}\) son o no distintos de cero:
Mostrar código R
summary(modeloGlobal)
Call:
lm(formula = log(peso) ~ IM + log(long) + IM * log(long), data = sargos)
Residuals:
Min 1Q Median 3Q Max
-0.155351 -0.052499 -0.008052 0.049076 0.261579
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -2.84858 0.11957 -23.824 <2e-16 ***
IM -0.08313 0.17684 -0.470 0.639
log(long) 2.58016 0.03952 65.295 <2e-16 ***
IM:log(long) 0.02396 0.05789 0.414 0.679
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.07709 on 196 degrees of freedom
Multiple R-squared: 0.9766, Adjusted R-squared: 0.9763
F-statistic: 2729 on 3 and 196 DF, p-value: < 2.2e-16
Como vemos, el p-valor para IM (que corresponde al coeficiente \(b_{0M}\)) es 0.639. Esto significa que puede aceptarse que dicho coeficiente es 0 en la población; dicho de otro modo, la diferencia entre las ordenadas de las rectas de machos y hembras no es significativa.
Asimismo el p-valor para IM:log(long) que corresponde al coeficiente \(b_{1M}\) es 0.679, también mayor que 0.05. Por tanto el coeficiente correspondiente (0.02396) es no significativo, lo que quiere decir que puede aceptarse que dicho coeficiente es cero en la población, o lo que es lo mismo, que no hay diferencias entre las pendientes de las rectas de machos y hembras.
De modo equivalente, el modelo ANCOVA puede ajustarse en R sin necesidad de definir la variable indicatriz IM. Para ello, se requiere que la variable que define los grupos (en este caso el sexo) sea un factor; la sintaxis a utilizar es entonces:
En lugar de realizar dos contrastes, podemos decidir si el modelo con dos rectas distintas (una para machos y otra para hembras) muestra diferencias significativas con el modelo con una sola recta (la misma para ambos sexos). Ello puede llevarse a cabo mediante la siguiente sintaxis:
Analysis of Variance Table
Model 1: log(peso) ~ log(long)
Model 2: log(peso) ~ log(long) * sexo
Res.Df RSS Df Sum of Sq F Pr(>F)
1 198 1.1708
2 196 1.1649 2 0.0059321 0.4991 0.6079
El p-valor obtenido (0.6079) es mayor que 0.05, lo que indica que no hay diferencias significativas entre ambos modelos. Por tanto, podemos considerar que la relación talla peso es la misma para ambos sexos.
Podemos utilizar la misma estrategia para decidir si la relación talla-peso difiere entre islas; para ello ajustamos un modelo de regresión distinto para cada isla (mediante un modelo ANCOVA) y lo comparamos con el modelo que presupone una recta común:
Analysis of Variance Table
Model 1: log(peso) ~ log(long)
Model 2: log(peso) ~ log(long) * isla
Res.Df RSS Df Sum of Sq F Pr(>F)
1 198 1.1708
2 186 1.0771 12 0.093736 1.3489 0.1943
También en este caso el p-valor (0.1943) mayor que 0.05 indica que ambos modelo son equivalentes. Por tanto podemos concluir que no hay diferencias significativas en la relación talla peso entre las distintas islas (o, dicho de otra forma, que la relación talla-peso es la misma en todas las islas).
Regresión logística
Abordamos por último el modelo de regresión logística. En este caso la variable respuesta \(y\) es una variable 0/1 (Presencia/Ausencia). El modelo de regresión logística es de la forma:
siendo \(f\left(s\right)=\frac{e^{s}}{1+e^{s}}\). Así pues, éste es un modelo que pretende predecir la probabilidad de que \(Y=1\) a partir de los valores tomados por las variables explicativas \(X_1, X_2,\dots X_p\) mediante una función (función de enlace o link function) que depende de una combinación lineal de los valores de dichas variables. Nótese que tal como se ha definido la función \(f\) sus valores están siempre entre 0 y 1, por lo que pueden interpretarse correctamente como probabilidades.
En nuestro ejemplo, podemos preguntarnos si la probabilidad de que un sargo esté parasitado depende del tamaño del pez, de la isla donde vive y de su sexo. El modelo logístico en este caso sería de la forma:
Aquí la función \(I_{\{isla=i\}}\) es una función que vale 1 si el pez se ha capturado en la isla \(i\) (donde \(i\) va desde la isla 2 hasta la 7) y 0 si se ha capturado en alguna otra isla; asimismo la variable \(I_M\) vale 1 si el pez es macho y 0 si es hembra. Nótese que si el pez tiene una longitud long, se ha capturado en la isla 1 y es hembra, el modelo anterior se reduce a:
La comparación de las dos expresiones anteriores permite intuir que, a igualdad de longitud y sexo, el valor de \(b_2\) evalúa el efecto de estar en la isla 2 con respecto a estar en la isla 1. Si \(b_2\) fuese 0, ello significaría que la probabilidad de tener parásitos en la isla 2 es la misma que en la isla 1 (cuando se mantienen constantes la longitud del pez y su sexo). Idéntico razonamiento puede hacerse con el resto de parámetros.
Este modelo puede estimarse en R mediante la sintaxis:
Call:
glm(formula = ptdo ~ long + isla + sexo, family = "binomial",
data = sargos)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.0786 -0.6090 -0.4724 -0.3251 2.5632
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -4.10625 1.31432 -3.124 0.00178 **
long 0.08222 0.05541 1.484 0.13782
islaGC -0.49849 0.76240 -0.654 0.51321
islaHI -0.79157 1.17804 -0.672 0.50162
islaLG -0.06664 1.20409 -0.055 0.95587
islaLP 0.48427 0.75316 0.643 0.52023
islaLZ 0.64686 0.71812 0.901 0.36771
islaTF 0.65295 0.66921 0.976 0.32921
sexoMacho 0.75041 0.42052 1.784 0.07435 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 169.08 on 199 degrees of freedom
Residual deviance: 157.16 on 191 degrees of freedom
AIC: 175.16
Number of Fisher Scoring iterations: 5
Como vemos, el único parámetro significativo es el \(b_0\); no hay efecto significativo de la longitud, de la isla ni del sexo (todos los p-valores son mayores que 0.05). Podemos concluir, por tanto, que la probabilidad de que un pez tenga parásitos no depende de su longitud, ni de la isla en que vive ni de su sexo.
Del mismo modo que hemos hecho mśa arriba con el ANCOVA, puede realizarse un único contraste para decidir si hay o no efecto de alguna de las variables consideradas. Si definimos el modelo nulo como:
es decir, un modelo en el que la probabilidad de tener parásitos es constante, podemos compararlo con el modelo logístico global anterior mediante (hay que especificar test="Chisq" pues la comparación en este caso se hace mediante un test de la chi-cuadrado):
Mostrar código R
anova(modeloNulo,modeloLogist,test="Chisq")
Analysis of Deviance Table
Model 1: ptdo ~ 1
Model 2: ptdo ~ long + isla + sexo
Resid. Df Resid. Dev Df Deviance Pr(>Chi)
1 199 169.08
2 191 157.16 8 11.925 0.1546
Como vemos, el p-valor del contraste global es 0.1546; esto significa que considerar que la probabilidad de tener parásitos depende del tamaño del pez, de la isla en que vive o de su sexo, no aporta ninguna diferencia significativa con respecto a considerar que dicha probabilidad es constante y no depende de ninguna de esas variables. Aplicando el principio de Ockham (entre dos modelos equivalentes es mejor elegir el más simple), concluímos que ninguna de esas variables influye en la probabilidad de tener parásitos.
