En este manual veremos los aspectos básicos de una regresión lineal, quizá el análisis estadístico más utilizado en ecología.

Una regresión lineal simple establece una relación lineal entre dos variables, donde la ecuación que la define es:

\(y=\alpha +\beta x+\epsilon\)

Donde:

  • \(y\) es la variable respuesta o dependiente.
  • \(x\) es la variable independiente o predictora.
  • \(\alpha\) es el intercepto y describe el valor que toma \(y\) cuando \(x\) es igual a cero. Normalmente carece de significado ecológico, por lo que nos centraremos en \(\beta\).
  • \(\beta\) es la pendiente. Describe en cuantas unidades aumenta \(y\) cuando \(x\) aumenta en una unidad.
  • \(\epsilon\) es el error y sigue una distribución normal \(N(0,\sigma ^2)\). Representa la varianza de la variable \(y\) para un valor de \(x\).

Para entender mejor el error (también conocido como residuales), podemos decir que para cierto valor de \(x\), el correspondiente valor de \(y\) va a estar distribuido con media \(\alpha + \beta x\) y con varianza \(\sigma ^2\). Mientras \(\sigma ^2\) sea mas grande, el ajuste de la línea de regresión será peor.

Si \(\beta >0\), entonces las variables \(x\) e \(y\) tienen una relación positiva, si \(\beta <0\), entonces las variables \(x\) e \(y\) tienen una relación negativa, y si \(\beta =0\), entonces las variables \(x\) e \(y\) no están relacionadas.

Regresión lineal simple. Línea roja es la línea de regresión y las líneas azules es el intervalo de confianza. Puntos son observaciones.

Regresión lineal simple. Línea roja es la línea de regresión y las líneas azules es el intervalo de confianza. Puntos son observaciones.

0.1 Ajuste de una línea de regresión simple

Normalmente se realiza por el método de mínimos cuadrados. Si asumimos:

\(S=\sum_{i=1}^n d_i^2\)

Donde \(d_i\) es la distancia entre observaciones \((x_i, y_i)\) y puntos sobre la línea de regresión \((x_i,\hat{y_i})\). Entonces, el objetivo es minimizar \(S\).

Asumamos que \(a\) y \(b\) son los mejores estimados que podemos obtener para \(\alpha\) y \(\beta\), entonces:

\(\hat{y_i} = a+bx_i\)

Vamos a realizar paso a paso el ajuste de una línea de regresión simple. Definamos los siguientes términos:

\(L_{xx}=\sum_{i=1}^n(x_i-\bar{x})^2\)

\(L_{yy}=\sum_{i=1}^n(y_i-\bar{y})^2\)

Estos son la suma de cuadrados corregidos. Ahora:

\(L_{xy}=\sum_{i=1}^n(x_i-\bar{x})(y_i-\bar{y})\)

Esta es la suma de productos cruzados corregidos. Nos damos cuenta que cuando \(L_{xy}\) es positivo, entonces \(\beta >0\), y cuando \(L_{xy}\) es negativo, entonces \(\beta <0\).

Se puede probar (dejaremos la prueba de esto para una clase de estadística) que las siguientes estimaciones minimizan la suma de cuadrados de los residuales:

\(b=\frac{L_{xy}}{L_{xx}}\)

\(a=\bar{y}-b\bar{x}\)

Donde \(\bar{y}\) y \(\bar{x}\) representa la media de los valores de \(y\) y \(x\), respectivamente. Con estos valores, podemos obtener los predichos:

\(\hat{y}=a+bx\)

Vemos que el punto \((x,a+bx)\) debe estar siempre sobre la línea de regresión. Ahora, vemos que:

\(y=a+bx = \bar{y}-b\bar{x}+bx=\bar{y}+b(x-\bar{x})\)

Entonces:

\(y-\bar{y}=b(x-\bar{x})\)

Podemos identificar dos componentes en la línea de regresión:

  1. Componente residual: \(y_i - \hat{y_i}\)
  2. Componente de regresión: \(\hat{y_i}-\bar{y}\)

Para tener un buen ajuste, siempre esperamos que el componente residual sea menor que el componente de regresión. Ahora definamos tres suma de cuadrados:

  1. Suma de cuadrados total (\(TSS\)): \(\sum_{i=1}^n(y_i - \bar{y})^2\)
  2. Suma de cuadrados de la regresión (\(RSS\)): \(\sum_{i=1}^n(\hat{y_i} - \bar{y})^2\)
  3. Suma de cuadrados de los residuos (\(ESS\)): \(\sum_{i=1}^n(y_i - \hat{y})^2\)

Y vemos que se cumple que \(TSS = RSS + ESS\).

0.1.1 Prueba F

Podemos implementar una prueba F con \(RSS\) y \(ESS\). Aquí \(H_0: \beta = 0\) y \(H_A:\beta \neq 0\), y el estadístico es:

\(F = \frac{RMS}{EMS}\),

Donde \(RMS=RSS/k\), y \(k\) es el número de variables independientes. Además, \(EMS =ESS/(n-k-1)\), y \(n\) es el número de observaciones.

En este caso, para un nivel de significancia \(\alpha\):

  • \(F>F_{1,n-2,1-\alpha}\) rechazamos \(H_0\)
  • \(F \leq F_{1,n-2,1-\alpha}\), entonces fallamos en rechazar \(H_0\)

0.1.2 One-way ANOVA

Esta es una prueba de hipótesis que se utiliza cuando se quiere comparar la media de alguna variable entre diferentes grupos o tratamientos. Se asume que la distribución por tratamiento de esta variable tiene una distribución normal y las mismas varianzas.

Ejemplo: Se tiene la variable longitud total de una muestra de peces que han sido sometidas a tres tratamientos de temperatura: baja, media, alta (3 tratamientos o grupos). La pregunta que se quiere responder con un ANOVA es si las medias de longitud total entre estos grupos es estadísticamente diferente o no. \(H_0: \mu_1=\mu_2=...=\mu_k\) y \(H_A: \mu_1\neq\mu_2\neq...\neq\mu_k\).

El procedimiento para realizar esta prueba de hipótesis es similar a la que acaba de ser descrita con la prueba F.

0.1.3 Prueba t

Esta es una prueba alternativa a la prueba F. Aquí también: \(H_0:\beta =0\) y \(H_0:\beta \neq0\). El estadístico es:

\(t = \frac{b}{(EMS/L_{xx})^2}\)

Para un nivel de significancia de \(\alpha\):

  • \(t>t_{n-2,1-\alpha/2}\), rechazamos \(H_0\)
  • \(-t_{n-2,1-\alpha /2} \leq t\leq t_{n-2,1-\alpha /2}\), entonces fallamos en rechazar \(H_0\).

Los resultados de esta prueba t es la que observamos para cada parámetro cuando hacemos summary(modelo), donde modelo representa la regresión lineal en R.

También podemos definir \(R^2\), un indicador muy utilizado en una regresión lineal que da información acerca de la bondad de ajuste del modelo. Es la proporción de la varianza de \(y\) que es explicada por \(x\).

\(R^2=\frac{RSS}{TSS}\)

Donde \(R^2\) toma valores entre 0 y 1.

El error estándar (una medida que nos dice que tan lejos están los parámetros estimados del valor real en la población) de los parámetros:

\(se(b)=\sqrt{\frac{EMS}{L_{xx}}}\)

\(se(a)=\sqrt{EMS(\frac{1}{n}+\frac{\bar{x}^2}{L_{xx}})}\)

Donde los intervalos de confianza para los parámetros estimados son:

\(b \pm t_{n-2,1-\alpha /2}se(b)\)

\(a \pm t_{n-2,1-\alpha /2}se(a)\)

Los intervalos de confianza (al \(100(1-\alpha)\%\)) para los valores sobre la línea de regresión estimada:

\(se_1(\hat{y}) = \sqrt{EMS[\frac{1}{n}+\frac{(x-\bar{x})^2}{L_{xx}}]}\)

\(\hat{y} \pm t_{n-2,1-\alpha /2}se_1(\hat{y})\)

Los intervalos de predicción (para valores que no fueron tomados al momento de hacer la regresión, es decir, valores de \(x\) donde se quiere predecir \(y\)) es (al \(100(1-\alpha)\%\)):

\(se_2(\hat{y})=\sqrt{EMS[1+\frac{1}{n}+\frac{(x-\bar{x})^2}{L_{xx}}]}\)

\(\hat{y} \pm t_{n-2,1-\alpha /2}se_2(\hat{y})\)

Es importante mencionar que los intervalos de predicción no son los mismos que los intervalos de confianza. Los intervalos de predicción define el rango en que un nuevo individuo (o valor de \(x\)) donde caerá. Por otro lado, el intervalo de confianza define el rango probable de valores asociados con los parámetros del modelo.

0.1.4 Supuestos de una regresión lineal

Existen tres supuestos importante que debemos tomar en cuenta al momento de hacer una regresión lineal:

  1. Para cualquier valor de \(x\), el valor correspondiente de \(y\) tiene un valor promedio \(\alpha + \beta x\), la cual es una función lineal de \(x\).
  2. Para cualquier valor de \(x\), el valor correspondiente de \(y\) esta normalmente distribuido alrededor de \(\alpha + \beta x\) con la misma varianza \(\sigma ^2\) para cualquier \(x\).
  3. Para dos puntos cualquiera \((x_1,y_1)\), \((x_2,y_2)\), los términos de error \(\epsilon_1, \epsilon_2\) son independientes uno del otro.
#Vamos a observar paso a paso como obtenemos todos los valores que
# acabamos de describir

# Activamos las librerias que vamos a utilizar
require(e1071)
require(MASS)
require(pracma)
require(maps)
require(mapdata)
require(corrplot)
require(PerformanceAnalytics)
require(ggplot2)
require(dplyr)

# Simular datos:
set.seed(100)
x = runif(100)
n = length(x)
y = 3*x+2+0.5*rnorm(n)

# Observemos lo que hemos simulado:
plot(x, y)

# Veamos la distribucion de la variable respuesta:
hist(y)

# Comenzemos a calcular los parametros paso por paso:
var_x = var(x)
var_y = var(y)

cov_xy = cov(x,y)

sd_x = sd(x)
sd_y = sd(y)

Lxx = sum((x-mean(x))^2)
Lyy = sum((y-mean(y))^2)
Lxy = sum((x-mean(x))*(y-mean(y)))

# Calculemos la pendiente y intercepto:
slope = Lxy/Lxx # pendiente
intercept = mean(y)-slope*mean(x) # intercepto

# Predicciones:
y_pred = intercept + slope*x

# Calcular residuales:
res = y - y_pred

# calculamos TSS:
tss = sum((y-mean(y))^2)

# Calcular RSS: (componente de regresion)
rss = sum((y_pred-mean(y))^2)

# Calcular ESS: (componente residual)
ess = sum(res^2)

k = 1 # Numero de variables independientes
# Calcular RMS:
rms = rss/k

# Calcular EMS:
ems = ess/(n-k-1)

# Corremos prueba F para evaluar la pendiente
# Pendiente de referencia:
beta_0 = 0 # H0: Beta = 0

# Estadistico:
F_stat = rms/ems

# Encontramos valor critico
critical_value = qf(p = 1-0.05, df1 = k, df2 = n-k-1)

F_stat > critical_value # Rechazamos la hipotesis nula?
## [1] TRUE
# Calculamos el p-value:
p_value = 1 - pf(q = F_stat, df1 = k, df2 = n-k-1) # acumulada (p-value)

# --------------
# Error estandar de pendiente
se_slope = sqrt(ems/Lxx)

# Error estandar de intercepto
se_intercept = sqrt(ems*(1/n + (mean(x)^2)/Lxx))


# Corremos prueba t para evaluar la pendiente (metodo alternativo)
# Pendiente de referencia:
beta_0 = 0 # H0: Beta = 0

# Estadistico:
t_stat_slope = (slope-beta_0)/se_slope

# Encontramos valor critico
critical_value = qt(p = 1-0.05/2, df = n-k-1)

abs(t_stat_slope) > critical_value # Rechazamos la hipotesis nula?
## [1] TRUE
# Calculamos el p-value:
p_value = 2*(1-pt(q = t_stat_slope, df = n-k-1)) # acumulada (p-value)

# De igual forma lo podemos hacer para el intercepto
# Pendiente de referencia:
intercept_0 = 0 # H0: Alpha = 0

# Estadistico:
t_stat_int = (intercept-intercept_0)/se_intercept

# Encontramos valor critico
critical_value = qt(p = 1-0.05/2, df = n-k-1)

abs(t_stat_int) > critical_value # Rechazamos la hipotesis nula?
## [1] TRUE
# Calculamos el p-value:
p_value = 2*(1-pt(q = t_stat_int, df = n-k-1)) # acumulada (p-value)

# Tambien podemos calcular los CI de los parametros estimados:
# CI de la pendiente:
slope - qt(p = 1 - 0.05/2, df = n-k-1)*se_slope # limite inferior
## [1] 2.582239
slope + qt(p = 1 - 0.05/2, df = n-k-1)*se_slope # limite superior
## [1] 3.353076
# CI de la intercepto:
intercept - qt(p = 1 - 0.05/2, df = n-k-1)*se_intercept # limite inferior
## [1] 1.748796
intercept + qt(p = 1 - 0.05/2, df = n-k-1)*se_intercept # limite superior
## [1] 2.197097
# Terminamos con la prueba de hipotesis

# Calculamos R2
R_sq = rss/tss

# Correlacion:
r = cov_xy/(sd_x*sd_y)
r^2 # Debe darnos R_sq
## [1] 0.7043571
# Tambien podemos calcular R_sq ajustado:
r_adj = sqrt(1-((1-r^2)*(n-1)/(n-2)))
Rsq_adj = r_adj^2

# Ahora calculamos el CI de los valores y usados en la regresion
se1_y_pred = sqrt(ems*((1/n)+((x-mean(x))^2)/Lxx))
y_up = y_pred + qt(p = 1-0.05/2, df = n-k-1)*se1_y_pred
y_low = y_pred - qt(p = 1-0.05/2, df = n-k-1)*se1_y_pred

# Hacemos una figura para observar CI
plot(x,y,main='Regresión lineal simple')
lines(sort(x),y_pred[order(x)], col = 'red')
lines(sort(x),y_up[order(x)],lty=2, col = 'blue')
lines(sort(x),y_low[order(x)],lty=2, col = 'blue')

# Ahora calculamos el error estandar de los valores y
# predichos basados en NUEVOS valores x (intervalos de prediccion)
preddata = data.frame(x = seq(from = 0, to = 1, by = 0.1))
se2_y_pred = sqrt(ems*(1+(1/n)+((preddata$x-mean(x))^2)/Lxx))
y_new_pred = intercept + slope*preddata$x
y_up_pred = y_new_pred + qt(p = 1-0.05/2, df = n-k-1)*se2_y_pred
y_low_pred = y_new_pred - qt(p = 1-0.05/2, df = n-k-1)*se2_y_pred

# Observemos intervalos de confianza y prediccion juntos:
plot(x,y,main='Intervalo de prediccion anadido')
lines(sort(x),y_pred[order(x)], col = 'red')
lines(sort(x),y_up[order(x)],lty=2, col = 'blue')
lines(sort(x),y_low[order(x)],lty=2, col = 'blue')
lines(preddata$x,y_up_pred,lty=2, col = 'green')
lines(preddata$x,y_low_pred,lty=2, col = 'green')

Lo mismo que hemos hecho hasta el momento, lo podemos obtener usando la función lm:

# Ahora usamos la funcion implementada en R para hacer exactamente 
# lo mismo que hemos hecho paso a paso:
lm1 = lm(y~x)
summary(lm1)
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.09375 -0.31748 -0.02187  0.26767  1.33092 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   1.9729     0.1130   17.47   <2e-16 ***
## x             2.9677     0.1942   15.28   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5063 on 98 degrees of freedom
## Multiple R-squared:  0.7044, Adjusted R-squared:  0.7013 
## F-statistic: 233.5 on 1 and 98 DF,  p-value: < 2.2e-16
anova(lm1) # Prueba F
## Analysis of Variance Table
## 
## Response: y
##           Df Sum Sq Mean Sq F value    Pr(>F)    
## x          1 59.859  59.859  233.48 < 2.2e-16 ***
## Residuals 98 25.125   0.256                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
predict(lm1, se.fit = TRUE) # arroja error estandar 
## $fit
##        1        2        3        4        5        6        7        8 
## 2.886291 2.737630 3.612050 2.140272 3.363440 3.408612 4.383879 3.071931 
##        9       10       11       12       13       14       15       16 
## 3.594945 2.478226 3.827722 4.590911 2.804940 3.155522 4.235937 3.958374 
##       17       18       19       20       21       22       23       24 
## 2.580165 3.033958 3.039745 4.021492 3.563050 4.082369 3.570581 4.195639 
##       25       26       27       28       29       30       31       32 
## 3.219664 2.481663 4.258938 4.590282 3.602477 2.797135 3.422071 4.728431 
##       33       34       35       36       37       38       39       40 
## 3.007745 4.804560 4.036282 4.612540 2.508333 3.840763 4.909634 2.359599 
##       41       42       43       44       45       46       47       48 
## 2.954233 4.540328 4.280551 4.428100 3.763406 3.430754 4.288783 4.597029 
##       49       50       51       52       53       54       55       56 
## 2.589370 2.884272 2.953846 2.562558 2.672406 2.788716 3.727785 2.724923 
##       57       58       59       60       61       62       63       64 
## 2.339414 2.655228 3.746345 2.600334 3.349053 3.893321 4.823598 3.980264 
##       65       66       67       68       69       70       71       72 
## 3.293993 3.034696 3.325401 3.294782 2.700297 4.033541 3.196325 2.945524 
##       73       74       75       76       77       78       79       80 
## 3.672122 4.842668 3.936880 3.826835 4.515199 4.272225 4.448053 2.244517 
##       81       82       83       84       85       86       87       88 
## 3.336661 3.751755 4.702366 4.889632 2.085131 3.688067 4.149172 2.711129 
##       89       90       91       92       93       94       95       96 
## 2.865429 4.149624 4.664476 2.595611 3.035777 3.303345 4.662910 3.128669 
##       97       98       99      100 
## 3.508590 2.344613 2.062409 4.263401 
## 
## $se.fit
##   [1] 0.06527416 0.07181178 0.05102425 0.10328042 0.05160521 0.05111658
##   [7] 0.07610311 0.05837265 0.05089810 0.08470035 0.05459504 0.08668825
##  [13] 0.06875821 0.05585254 0.06917555 0.05833426 0.07945240 0.05964818
##  [19] 0.05944883 0.06049188 0.05072804 0.06276044 0.05076050 0.06740631
##  [25] 0.05421447 0.08452014 0.07020983 0.08665485 0.05095066 0.06910482
##  [31] 0.05100323 0.09414031 0.06057250 0.09837700 0.06102686 0.08784106
##  [37] 0.08312915 0.05491992 0.10433284 0.09104231 0.06256314 0.08402346
##  [43] 0.07119701 0.07828731 0.05316444 0.05093805 0.07157680 0.08701355
##  [49] 0.07898910 0.06535761 0.06257805 0.08034376 0.07489920 0.06948098
##  [55] 0.05250067 0.07240387 0.09214309 0.07573142 0.05283473 0.07843969
##  [61] 0.05179531 0.05634144 0.09944734 0.05905877 0.05267215 0.05962264
##  [67] 0.05214327 0.05265795 0.07356451 0.06092692 0.05477895 0.06289961
##  [73] 0.05165741 0.10052351 0.05764866 0.05457337 0.08271700 0.07081500
##  [79] 0.07928777 0.09739168 0.05197217 0.05293694 0.09270676 0.10319025
##  [85] 0.10644037 0.05187422 0.06543875 0.07305191 0.06614425 0.06545752
##  [91] 0.09063971 0.07867606 0.05958532 0.05250683 0.09055468 0.05661673
##  [97] 0.05063564 0.09185902 0.10775075 0.07041249
## 
## $df
## [1] 98
## 
## $residual.scale
## [1] 0.5063347
# Calculamos los valores predichos con nuevos valores de x:
predict(lm1, newdata = preddata, interval = 'predict')
##         fit       lwr      upr
## 1  1.972946 0.9434434 3.002449
## 2  2.269712 1.2470114 3.292413
## 3  2.566478 1.5491637 3.583792
## 4  2.863244 1.8498777 3.876609
## 5  3.160009 2.1491366 4.170882
## 6  3.456775 2.4469297 4.466620
## 7  3.753541 2.7432523 4.763829
## 8  4.050307 3.0381065 5.062507
## 9  4.347072 3.3315006 5.362644
## 10 4.643838 3.6234489 5.664227
## 11 4.940604 3.9139720 5.967236
# Observemos algunas figuras para analizar las asunciones de este metodo
opar = par(mfrow=c(2,2)) # panel 2 por 2
plot(lm1)

par(opar) # vuelve a la configuracion de plot inicial

# Como interpretamos estos resultados?:
# Intercepto: cuando y toma un valor de 0, que valor toma x?
# Pendiente: cada unidad en incremento en x, cuanto crece (o decrece) y

0.1.5 Caso de estudio

0.1.5.1 Datos de salmones en el Pacífico norte

Vamos a implementar algunas regresiones lineales, donde la variable respuesta es la productividad de un stock de salmón y las variables independientes son la biomasa de desovantes y temperatura.

# Iniciamos con datos aplicados:
chum = read.csv('data/coastwide_chum_data.csv')

# Elegir solo los datos que necesitamos:
stock1 = chum[chum$stock=='SE-SC AKPen',c('stock','entry.yr','spawners','recruits','ln.rs',
                                         'loc.sst','w.sst.1','w.sst.2','npgo2','npgo','pdo2','pdo','long','lat')]
head(stock1)
##           stock entry.yr spawners recruits       ln.rs     loc.sst    w.sst.1
## 764 SE-SC AKPen     1963   238600   576200  0.88166631  0.04182643 -0.2207373
## 765 SE-SC AKPen     1964   263000   251500 -0.04471104 -1.19076329  0.6549483
## 766 SE-SC AKPen     1965   160800   313100  0.66636127 -0.47777858  1.0118060
## 767 SE-SC AKPen     1966   203300   192000 -0.05718735 -1.02832160 -1.1020202
## 768 SE-SC AKPen     1967   354800   657883  0.61747285  0.59927106 -0.3414790
## 769 SE-SC AKPen     1968   132800   846192  1.85190205  0.20264870 -0.5955610
##         w.sst.2        npgo2       npgo   pdo2    pdo long  lat
## 764  0.21710551  0.164053783 -0.2557230 -0.659 -0.472  161 55.5
## 765  0.83337719 -0.007529911  0.5838306 -0.891 -0.846  161 55.5
## 766 -0.04510709 -0.714296310 -0.5988904 -0.735 -0.936  161 55.5
## 767 -0.72174962 -0.809721169 -0.8297022 -0.572 -0.534  161 55.5
## 768 -0.46852000 -0.679363840 -0.7897401 -0.516 -0.610  161 55.5
## 769 -0.32856791 -0.796891715 -0.5689875 -0.653 -0.422  161 55.5
# Exploremos que datos tenemos en el mapa:
plot(-stock1$long,stock1$lat,
     ylab='lat',
     xlab='lon',
     ylim=range(chum$lat),
     xlim=range(-chum$long),
     pch=16,main='Chum salmon',col='red',cex=2)
map("world", fill=T, col="lightblue4", add=T)

# Vamos a implementar muchos modelos lineales:
# La variable Y sera ln.rs (indicador de productividad)
lm_0 = lm(ln.rs ~ spawners, data=stock1)
lm_pdo = lm(ln.rs ~ spawners + pdo, data=stock1)
lm_loc_sst = lm(ln.rs ~ spawners + loc.sst, data=stock1)

# Hagamos resumen de algunos modelos
summary(lm_0)
## 
## Call:
## lm(formula = ln.rs ~ spawners, data = stock1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.08398 -0.26782  0.08663  0.29962  0.76672 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.387e+00  1.675e-01   8.279 1.61e-10 ***
## spawners    -1.771e-06  5.332e-07  -3.322  0.00181 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4526 on 44 degrees of freedom
## Multiple R-squared:  0.2005, Adjusted R-squared:  0.1823 
## F-statistic: 11.03 on 1 and 44 DF,  p-value: 0.001807
summary(lm_pdo)
## 
## Call:
## lm(formula = ln.rs ~ spawners + pdo, data = stock1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.99847 -0.24559 -0.01344  0.25640  0.77313 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.514e+00  1.543e-01   9.813 1.53e-12 ***
## spawners    -2.202e-06  4.929e-07  -4.467 5.67e-05 ***
## pdo          2.341e-01  6.762e-02   3.462  0.00122 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4048 on 43 degrees of freedom
## Multiple R-squared:  0.3748, Adjusted R-squared:  0.3457 
## F-statistic: 12.89 on 2 and 43 DF,  p-value: 4.117e-05
summary(lm_loc_sst)
## 
## Call:
## lm(formula = ln.rs ~ spawners + loc.sst, data = stock1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.87700 -0.27362  0.05539  0.21381  0.68280 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.534e+00  1.531e-01  10.021 8.11e-13 ***
## spawners    -2.273e-06  4.897e-07  -4.641 3.25e-05 ***
## loc.sst      2.456e-01  6.675e-02   3.680 0.000646 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3992 on 43 degrees of freedom
## Multiple R-squared:  0.3919, Adjusted R-squared:  0.3637 
## F-statistic: 13.86 on 2 and 43 DF,  p-value: 2.264e-05
# Observemos algunas figuras para analizar las asunciones de este metodo
opar = par(mfrow=c(2,2)) # panel 2 por 2
plot(lm_0)

plot(lm_pdo)

plot(lm_loc_sst)

par(opar) # vuelve a la configuracion de plot inicial

# Podemos utilizar un indicador para comparar modelos:
AIC(lm_0, lm_pdo, lm_loc_sst) # elegimos el valor mas bajo
##            df      AIC
## lm_0        3 61.55952
## lm_pdo      4 52.24755
## lm_loc_sst  4 50.96795

0.1.5.2 Datos de dieta de aves

Vamos a implementar una regresión lineal donde tendremos que tranformar la variable para que los supuestos se cumplan.

# Leemos los datos:
bird = read.csv("data/birdsdiet.csv") 
head(bird)
##               Family   MaxAbund   AvgAbund       Mass        Diet Passerine
## 1 Hawks&Eagles&Kites   2.985544  0.6740216  716.01894  Vertebrate         0
## 2   Long-tailed tits  37.800000  4.0404300    5.30000      Insect         1
## 3              Larks 241.400000 23.1054000   35.77709 PlantInsect         1
## 4        Kingfishers   4.400000  0.5949060  119.37172  Vertebrate         0
## 5      Auks& Puffins   4.532505  2.9632718  315.49552  InsectVert         0
## 6       Ducks& Geese  23.693635  2.7416729 1144.29003 PlantInsect         0
##   Aquatic
## 1       0
## 2       0
## 3       0
## 4       0
## 5       1
## 6       1
# Resumen rapido de mi base de datos
summary(bird) 
##     Family             MaxAbund          AvgAbund           Mass        
##  Length:54          Min.   :  0.200   Min.   : 0.200   Min.   :   5.30  
##  Class :character   1st Qu.:  5.402   1st Qu.: 1.340   1st Qu.:  20.68  
##  Mode  :character   Median : 24.147   Median : 3.114   Median :  59.18  
##                     Mean   : 44.906   Mean   : 5.686   Mean   : 468.48  
##                     3rd Qu.: 43.581   3rd Qu.: 6.258   3rd Qu.: 461.73  
##                     Max.   :413.600   Max.   :47.598   Max.   :5296.23  
##      Diet             Passerine        Aquatic      
##  Length:54          Min.   :0.000   Min.   :0.0000  
##  Class :character   1st Qu.:0.000   1st Qu.:0.0000  
##  Mode  :character   Median :0.000   Median :0.0000  
##                     Mean   :0.463   Mean   :0.2778  
##                     3rd Qu.:1.000   3rd Qu.:1.0000  
##                     Max.   :1.000   Max.   :1.0000
# Exploremos las variables que queremos modelar:
plot(bird$Mass, bird$MaxAbund)

# Implementamos el modelo (regresion) linear:
lm1 = lm(bird$MaxAbund ~ bird$Mass)

# Exploremos los resultados
opar = par(mfrow=c(2,2)) 
plot(lm1)

par(opar)

# Plot Y ~ X y la linea de regresion:
ggplot(bird, aes(x = Mass,y = MaxAbund))+
  geom_point()+
  geom_smooth(method="lm", colour = 'red') +
  ylab('MaxAbundance') +
  xlab('Mass') +
  ggtitle('Regresion') +
  theme_bw()

# Veamos si los datos estan normalmente distribuidos
hist(bird$MaxAbund,col="coral", main="Untransformed data", 
     xlab="Maximum Abundance")

hist(bird$Mass, col="coral", main="Untransformed data", xlab="Mass")

# Podemos correr una prueba de hipotesis para saber si 
# los datos vienen de una distribucion normal
# H0: la distribucion es normal
shapiro.test(bird$MaxAbund) 
## 
##  Shapiro-Wilk normality test
## 
## data:  bird$MaxAbund
## W = 0.5831, p-value = 3.872e-11
shapiro.test(bird$Mass) 
## 
##  Shapiro-Wilk normality test
## 
## data:  bird$Mass
## W = 0.54667, p-value = 1.155e-11
# Si p < 0.05, la distribucion NO es normal
# Si p > 0.05, la distribucion es normal

# Podemos analizar tambien la simetria de mis datos
skewness(bird$MaxAbund) 
## [1] 3.257217
skewness(bird$Mass) 
## [1] 3.23552
# datos positivos: distribucion hacia la izquiera
# datos negativos: distribucion hacia la derecha

# Vamos a transformar mis datos con log (logaritmo natural)
bird$logMaxAbund = log(bird$MaxAbund)
bird$logMass = log(bird$Mass)
names(bird) 
## [1] "Family"      "MaxAbund"    "AvgAbund"    "Mass"        "Diet"       
## [6] "Passerine"   "Aquatic"     "logMaxAbund" "logMass"
# Veamos los nuevos datos:
hist(bird$logMaxAbund,col="yellowgreen", main="Log transformed", 
     xlab=expression("log"[10]*"(Maximum Abundance)"))

hist(bird$logMass,col="yellowgreen", main="Log transformed",
     xlab=expression("log"[10]*"(Mass)"))

# Correr nuevamente la prueba de hipotesis de normalidad y simetria:
shapiro.test(bird$logMaxAbund)
## 
##  Shapiro-Wilk normality test
## 
## data:  bird$logMaxAbund
## W = 0.96708, p-value = 0.1431
skewness(bird$logMaxAbund)
## [1] -0.6089474
shapiro.test(bird$logMass)
## 
##  Shapiro-Wilk normality test
## 
## data:  bird$logMass
## W = 0.93716, p-value = 0.007107
skewness(bird$logMass)
## [1] 0.4768559
# Corremos nuevamente el modelo:
lm2 = lm(bird$logMaxAbund ~ bird$logMass)

# Veamos las figuras para verificar las asunciones:
opar = par(mfrow=c(2,2))
plot(lm2, pch=19, col="gray")

par(opar)

# Parece todo bien, ahora veamos los resultados del modelo:
summary(lm2)
## 
## Call:
## lm(formula = bird$logMaxAbund ~ bird$logMass)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.4569 -0.9206  0.1263  0.9354  3.7180 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    3.8508     0.5691   6.767 1.17e-08 ***
## bird$logMass  -0.2361     0.1170  -2.019   0.0487 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.602 on 52 degrees of freedom
## Multiple R-squared:  0.07267,    Adjusted R-squared:  0.05484 
## F-statistic: 4.075 on 1 and 52 DF,  p-value: 0.04869
# Podemos obtener informacion a partir de este objeto
summary(lm2)$coefficients 
##                Estimate Std. Error   t value     Pr(>|t|)
## (Intercept)   3.8507680  0.5690882  6.766557 1.166186e-08
## bird$logMass -0.2361498  0.1169836 -2.018658 4.869342e-02
summary(lm2)$r.squared 
## [1] 0.07267022
summary(lm2)$adj.r.squared 
## [1] 0.05483696
# Ahora podemos graficar la regresion:
plot(logMaxAbund ~ logMass, data=bird, pch=19, col="yellowgreen", 
     ylab = expression("log"[10]*"(Maximum Abundance)"), xlab = expression("log"[10]*"(Mass)"))
abline(lm2, lwd=2)
# Anadimos intervalos de confianza:
confit = predict(lm2, interval="confidence")
points(bird$logMass,confit[,2]) 
points(bird$logMass,confit[,3])

# Podemos correr otro modelo excluyendo una parte de los datos:
lm3 = lm(logMaxAbund ~ logMass, data=bird, subset =! bird$Aquatic) # remover aves acuaticas

# Examine the model
opar = par(mfrow=c(2,2))
plot(lm3, pch=19, col=rgb(33,33,33,100,maxColorValue=225))

summary(lm3)
## 
## Call:
## lm(formula = logMaxAbund ~ logMass, data = bird, subset = !bird$Aquatic)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.1053 -0.5327  0.0928  0.5280  3.8709 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   5.2272     0.6750   7.744 2.96e-09 ***
## logMass      -0.6429     0.1746  -3.683 0.000733 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.403 on 37 degrees of freedom
## Multiple R-squared:  0.2682, Adjusted R-squared:  0.2485 
## F-statistic: 13.56 on 1 and 37 DF,  p-value: 0.000733
par(opar)

# Comparamos ambas regresiones
opar = par(mfrow=c(1,2))
plot(logMaxAbund ~ logMass, data=bird, main="All birds", ylab = expression("log"[10]*"(Maximum Abundance)"), 
     xlab = expression("log"[10]*"(Mass)"))
abline(lm2,lwd=2)

plot(logMaxAbund ~ logMass, data=bird, subset=!bird$Aquatic, main="Terrestrial birds",
     ylab = expression("log"[10]*"(Maximum Abundance)"), xlab = expression("log"[10]*"(Mass)"))
abline(lm3,lwd=2)

par(opar)

0.2 Correlación

La covarianza entre dos variables (\(X\) e \(Y\)) es una medida para cuantificar la relación entre ellas:

\(Cov(X,Y) = E[(X-\mu_x)(Y-\mu_y)] = E[XY] - \mu_x\mu_y\)

Si \(X\) e \(Y\) son independientes, entonces \(Cov(X,Y) = 0\).

Coeficiente de correlación se calcula de la siguiente manera:

\(\rho = Corr(X,Y) = \frac{Cov(X,Y)}{\sigma_X \sigma_Y}\)

Este coeficiente de correlación no tiene dimensión (no tiene unidades) y toma valores entre -1 y 1. El coeficiente de correlación de Pearson es simplemente un estimado de \(\rho\):

\(r=\frac{L_{xy}}{\sqrt{L_{xx}L_{yy}}}\)

Si \(r\) toma valores cercanos a -1, decimos que hay una correlación negativa muy fuerte. Si toma valores cercanos a 1, decimos que hay una correlación positiva muy fuerte. Si \(r\) toma valores cercanos a cero, decimos que las variables no están correlacionadas.

0.2.1 Prueba t

También se puede realizar una prueba t para observar la significancia de \(r\) y sus intervalos de confianza. Aquí \(H_0: \rho = 0\) y \(H_A:\rho \neq 0\). El estadístico es:

\(t=\frac{r(n-2)^{1/2}}{(1-r^2)^{1/2}}\)

Para un nivel de significancia \(\alpha\):

  • Si \(t>t_{n-2,1-\alpha /2}\) o \(t<-t_{n-2,1-\alpha /2}\), entonces se rechaza \(H_0\).
  • Si \(-t_{n-2,1-\alpha /2}<t<t_{n-2,1-\alpha /2}\), entonces fallamos en rechazar \(H_0\)

0.2.2 Correlación de rango (Rank correlation)

Este tipo de correlación se emplea para variables donde no se tiene seguro que siguen una distribución normal. El coeficiente de correlación de Spearman (\(r_s\)):

\(r_s=\frac{L_{xy}}{\sqrt{L_{xx}L_{yy}}}\)

Donde los \(L\) son calculados a partir de los rangos en vez de los valores como lo hicimos anteriormente.

Vamos a ver algunas funciones para observar la correlación entre dos variables:

# Prueba de hipotesis para correlacion:
cor.test(x, y)
## 
##  Pearson's product-moment correlation
## 
## data:  x and y
## t = 15.28, df = 98, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.7697282 0.8891116
## sample estimates:
##       cor 
## 0.8392598
# Podriamos tambien correr la prueba para uun tipo Spearman
# Esto no asume normalidad
cor.test(x, y, method=c("spearman"))
## 
##  Spearman's rank correlation rho
## 
## data:  x and y
## S = 25024, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##      rho 
## 0.849841
# Ahora veamos un plot de correlacion entre variables ambientables 
# tomados de la base de datos de salmon
# Hacemos un daata frame de solo variables ambientales:
envdata = chum[,c('loc.sst','w.sst.1','w.sst.2','npgo2','npgo','pdo2','pdo')]
corMat = cor(envdata)
corrplot(corMat, method="circle")

corrplot(corMat, method="number")

chart.Correlation(envdata, histogram=TRUE, pch=19)

0.3 Transformación de una variable

Se utilizan transformaciones cuando la variable que queremos incluir en una regresión lineal no tiene una distribución cercana a la normal. Tenemos que tener cuidado para estos casos, ya que la interpretación de los parámetros estimados va a ser diferente. Tenemos los siguientes casos:

  1. La variable dependiente es log-transformada: Para este caso, \((e^\beta - 1)100\%\) es el porcentaje de incremento (o reducción) en la variable respuesta por cada unidad de incremento en la variable independiente.
  2. La variable independiente es log-transformada: El 1% de incremento en la variable independiente incrementa o reduce la variable dependiente en \(\beta /100\) unidades.
  3. Ambas variables son log-transformadas: \(\beta\) es el incremento (o reducción) en porcentaje en la variable dependiente para cada 1% de incremento en la variable independiente.

0.4 Regresión múltiple

Para casos donde tenemos más de una variable independiente:

\(y=\alpha + \sum_{j=1}^k \beta _jx_j+\epsilon\)

Donde \(\beta_j\) son conocidos como coeficientes de regresión parcial. Las variables independientes pueden ser continuas o categóricas. Para este caso, la interpretación de los coeficientes no cambia. Internamente, R asigna valores de 0, 1, 2, … a cada nivel de la variable categórica, y debemos tener presente qué categoría tiene cada número para poder hacer una interpretación adecuada.

Existen casos donde se especifica un efecto de interacción, por ejemplo:

\(y=\alpha + \beta _1 x_{1,continua} + \beta _2 x_{2,categorica}+\beta _3 x_{1,continua}*x_{2,categorica}\)

El símbolo \(*\) significa un efecto de interacción. Para este caso, el efecto de \(x_{1,continua}\) sobre \(y\) va a ser diferente para distintos niveles de \(x_{2,categorica}\). Por lo tanto, no solo tenemos que ver el valor de \(\beta_1\), si no también el de \(\beta_3\). Por ejemplo, el efecto de \(x_{1,continua}\) sobre \(y\) va a ser \(\beta_1+\beta_3x_{2,categorica}\), luego reemplazamos \(x_{2,categorica}\) por sus respectivos valores 0,1,2,.. y procedemos a obtener conclusiones.

Para estos casos, \(\beta_2\) no tiene una interpretación lógica y se suelen ignorar.

En los ejercicios anteriores ya vemos visto algunos casos de regresión múltiple.

0.5 Validación de un modelo

En esta sección vamos a realizar un ejercicio simple de como validar un modelo mediante validación cruzada. Este método se basa en implementar un modelo estadístico (e.g. regresión lineal, modelo lineal generalizado, modelo aditivo generalizado, etc.) con un porcentaje de datos observados, y luego hacer predicciones usando las observaciones restastes, comparando los valores predichos con las observaciones de la variable respuesta.

#Vamos a usar las variables simuladas anteriormente:
obs_data = data.frame(x=x, y=y)

#Seleccionamos 80 observaciones y reservamos 20 observaciones al azar:
n = 100 # numero de observaciones
ind = sample(c(TRUE, FALSE), n, replace=TRUE, prob=c(0.8, 0.2))
data1 = obs_data[ind, ]
data2 = obs_data[!ind, ]

# Construimos la regresion con el 80%
lm1 = lm(y~x, data = data1)

# Predecimos valores usando las observaciones del 20% reservado
pred_data = predict(lm1, newdata = data.frame(x = data2$x), interval = 'predict') 

# Encontramos la correlacion entre lo predicho y observado en el 20% reservado:
# Debe estar muy cercano a 1:
cor(pred_data[,1], data2$y)
## [1] 0.9305692
# Estos pasos se deben hacer muchas veces (iteraciones). La correlacion debe estar
# cercana a 1 para que el modelo sea robusto.