Medias y desviaciones tipicas por fecha y playa

Las Canteras:

date mean sd
2015-09-03 0.003867 0.002811
2015-09-06 0.0094 0.007252
2015-09-23 0.0011 0.0009539
2015-09-28 4.512 2.437
2015-10-05 0.1079 0.1006
2015-10-14 2.576 2.322
2015-10-21 1.444 0.7444
2015-10-27 1.784 2.556
2015-11-05 4.49 3.506
2015-11-12 1.013 0.5085
2015-11-19 0.03503 0.05602
2015-11-26 1.384 1.463
2015-12-05 0.0916 0.03098
2015-12-18 0.07943 0.1079
2015-12-28 0.3656 0.1759
2016-01-14 0.4681 0.3082
2016-01-19 0.3168 0.2585
2016-01-28 0.1507 0.1516
2016-02-22 1.642 0.7351
2016-03-03 0.03807 0.03472
2016-03-24 0.6737 0.7101
2016-04-10 3.779 4.181
2016-04-23 0.4333 0.3228
2016-05-10 1.613 0.1832
2016-05-21 0.9436 0.2536
2016-06-06 22.33 9.997
2016-06-21 0.4057 0.1842
2016-07-04 4.423 0.4383
2016-07-23 1.299 0.7099
2016-08-04 6.398 4.877
2016-08-19 14.26 14.66
2016-09-05 0.1474 0.07182

Famara

  date mean sd
33 2015-08-04 0.1245 0.01586
34 2015-08-08 0.03947 0.04274
35 2015-08-14 0.2088 0.1562
36 2015-08-24 0.1165 0.1616
37 2015-08-28 0.1713 0.2114
38 2015-09-14 0.163 0.1432
39 2015-09-28 1.9 1.752
40 2015-10-14 1.87 0.8129
41 2015-10-28 65.45 17.21
42 2015-11-14 5.5 2.106
43 2015-11-25 5.488 2.268
44 2015-12-12 9.677 4.432
45 2015-12-21 5.218 0.9616
46 2016-01-18 3.093 1.215
47 2016-01-29 5.301 3.034
48 2016-02-09 2.706 1.501
49 2016-02-22 3.995 1.149
50 2016-03-08 2.406 0.4666
51 2016-03-24 2.164 1.231
52 2016-04-09 9.894 2.094
53 2016-04-21 13.23 1.712
54 2016-05-04 10.36 7.183
55 2016-05-22 4.73 1.062
56 2016-06-03 2.078 0.8195
57 2016-06-20 1.286 0.246
58 2016-07-04 0.6516 0.3278
59 2016-07-22 0.32 0.1268
60 2016-08-02 0.741 0.3846
61 2016-08-20 0.6443 0.4106
62 2016-09-02 0.6014 0.5662

Lambra

  date mean sd
63 2015-08-03 1.107 0.01651
64 2015-09-29 9.319 4.376
65 2015-10-29 38.79 9.962
66 2015-11-15 3.2 0.8176
67 2015-11-27 30.31 6.501
68 2015-12-11 20.9 4.757
69 2015-12-22 12.74 4.355
70 2016-01-19 2.102 0.7116
71 2016-02-01 11.33 5.761
72 2016-02-10 5.505 2.062
73 2016-02-23 11.34 3.915
74 2016-03-09 61.09 10.18
75 2016-03-26 7.12 1.944
76 2016-04-11 18.75 4.032
77 2016-04-22 1.23 1.239
78 2016-05-05 1.413 0.8875
79 2016-05-23 1.661 1.865
80 2016-06-04 4.438 4.17
81 2016-06-21 1.047 0.3469
82 2016-07-05 2.653 2.312
83 2016-07-21 0.9474 0.4276
84 2016-08-03 1.447 0.1545
85 2016-08-21 1.955 0.8507
86 2016-09-03 4.623 2.845
87 2016-09-17 4.67 1.637

 

 

Gráficos:

Mostramos la representación gráfica de los datos anteriores:

 

 

Modelo para evaluar la cantidad total de microplásticos según la estación.

 

Las Canteras

Mostramos en primer lugar un resumen de datos (media y desviación típica) por estación en las Canteras:

season mean sd
spring 4.963 8.893
summer 3.368 6.605
autumn 1.46 2.109
winter 0.4969 0.6226

El mismo resumen por mes:

monthYear mean sd ndays
2015-09 1.132 2.288 4
2015-10 1.478 1.772 4
2015-11 1.73 2.389 4
2015-12 0.1789 0.1747 3
2016-01 0.3119 0.2552 3
2016-02 1.642 0.7351 1
2016-03 0.3559 0.5687 2
2016-04 2.106 3.224 2
2016-05 1.278 0.4167 2
2016-06 11.37 13.57 2
2016-07 2.861 1.791 2
2016-08 10.33 10.68 2
2016-09 0.1474 0.07182 1

Obsérvese en la gráfica anterior que de septiembre de 2015 a mayo de 2016 apenas hay cambio en el promedio del volumen total de microplásticos por mes. Se aprecia una entrada fuerte de microplásticos en el mes de junio de 2016, que desciende en julio y vuelve a subir en agosto. Ahora bien, en cada uno de estos tres meses se muestreó sólo dos días; no hay datos de los mismos meses en otros años (en realidad, el muestreo cubre escasamente un año; el único mes muestreado dos veces, en 2015 y 2016 fue septiembre). Por tanto no es posible determinar si hay una tendencia a que haya mayor llegada de microplásticos en estos meses.

Los gráficos siguientes muestran si se percibe alguna asociación del volumen de microplásticos registrado con alguna de las otras variables; para ello, en cada caso se muestra la recta de regresión para cada estación, que permite visualizar las posibles tendencias.

Si eliminamos los outliers:

 

 

Para determinar si el incremento en el volumen total de microplásticos que se registra en primavera-verano de 2016 muestra diferencias significativas con el resto de observaciones aplicaremos un modelo lineal de efectos mixtos. En primer lugar mostramos un boxplot con la cantidad total de microplásticos por estación:

Como vemos, hay bastante asimetría y aparecen los outliers que ya hemos comentado repetidas veces. Para poder aplicar el modelo aplicamos la siguiente transformación a los valores de la variable total (\(y\)):

\[y^{*}=\frac{1}{\left(y+1\right)^{2}}\]

La variable así transformada presenta un aspecto más regular:

El resultado del análisis es el siguiente:

  numDF denDF F-value p-value
(Intercept) 1 64 55.96 2.654e-10
season 3 28 1.013 0.4017

Como puede verse, el p-valor obtenido indica que las diferencias entre estaciones no son significativas (como cabía esperar dado el diseño del muestreo, básicamente que sólo hay datos de un año, y en concreto de un par de días por mes de ese año)

Para validar el modelo comprobamos homoscedasticidad (test de Levene) y normalidad (test de Shapiro-Wilk) en los residuos:

Levene.test.P Shapiro.test.P
0.9442 0.2722

 

 

 

Lambra

Mostramos en primer lugar un resumen de datos (media y desviación típica) por estación en Lambra:

season mean sd
spring 5.769 6.74
summer 2.306 1.92
autumn 20.5 14.48
winter 17.35 20.99

El mismo resumen por mes:

monthYear mean sd ndays
2015-08 1.107 0.01651 1
2015-09 9.319 4.376 1
2015-10 38.79 9.962 1
2015-11 16.76 15.42 2
2015-12 16.82 6.054 2
2016-01 2.102 0.7116 1
2016-02 9.392 4.657 3
2016-03 34.1 30.28 2
2016-04 9.991 9.961 2
2016-05 1.537 1.313 2
2016-06 2.742 3.233 2
2016-07 1.8 1.756 2
2016-08 1.701 0.6135 2
2016-09 4.647 2.076 2

Los datos disponibles son aproximadamente quincenales, salvo las dos primeras observaciones (3 de agosto y 29 de septiembre, entre las que hay casi dos meses) y la tercera (casi un mes):

## Time differences in days
##  [1] 57 30 17 12 14 11 28 13  9 13 15 17 16 11 13 18 12 17 14 16 13 18 13
## [24] 14

Considerando que aproximadamente hay quince días entre observaciones, podemos trazar gráficos de autocorrelación y autocorrelación parcial para determinar la posible estructura temporal de la serie (excluímos las dos primeras observaciones por distar dos meses y un mes, respectivamente, del resto):

Estos gráficos nos indican que no hay correlación temporal entre las sucesivas observaciones.

En esta playa, al contrario de lo que sucede en las Canteras, los meses “estables” ocurren a partir de mayo de 2016. Los meses anteriores presentan mucha más variabilidad. De hecho hay varios meses en los que sólo se muestreó un día, lo que efectivamente hace esperar gran variabilidad entre medidas. Como en el caso anterior, al abarcar el periodo de muestreo sólo un año de datos, no es posible hablar de tendencias estacionales.

Como ya se ha hecho para las Canteras, los gráficos siguientes muestran si se percibe alguna asociación del volumen total de microplásticos registrado con alguna de las otras variables:

La aparente tendencia que parece haber en invierno tiene más que ver con el outlier que con otra cosa; de hecho, si suprimimos el outlier, la tendencia desaparece:

Aquí sí parece haber cierta tendencia creciente (mayor volumen de total de microplásticos con mayor altura de ola) durante el otoño.

Si presentamos un ajuste global a todos los datos sin distinguir estación:

la recta así obtenida presenta una pendiente (\(\hat{\beta_1}=14.98\))) significativa (p=0.00000001)

 

Si ajustamos una tendencia global sin distinguir entre estaciones, observamos que esta es creciente:

El valor estimado de la pendiente de esta recta (\(\hat{\beta_1}=1.12\))) es significativo (p=0.0656)

 

aunque es muy probable que esta asociación se deba más a la presencia de outliers (hay dos valores muy bajos de marea) que a otra cosa.

(la asociación aparente en invierno está relacionada con la presencia de un outlier)

 

 

Las siguientes figuras muestran la variación conjunta del total de microplásticos (promediado para cada día de observación) con la velocidad del viento, la altura de ola y el periodo de ola. Se muestran los datos tanto en escala natural:

y en escala logarítmica:

 

 

Para determinar si existe variación significativa en el volumen total de microplásticos que se registra durante las distintas estaciones en la playa de Lambra aplicaremos nuevamente un modelo lineal de efectos mixtos. En primer lugar mostramos un boxplot con la cantidad total de microplásticos por estación:

Como vemos, hay bastante asimetría y aparecen algunos outliers a los que ya nos hemos referido. Para poder aplicar el modelo aplicamos una transformación logarítmica a la variable total. La variable así transformada presenta un aspecto más regular:

Recordemos que habíamos detectado ya una cierta asociación con el periodo de ola. Representamos ahora todos los valores (todas las réplicas de todos los días, no sólo las medias mensuales como antes):

 

Si analizamos el efecto de la estación sin tener en cuenta el periodo de ola obtenemos:

  numDF denDF F-value p-value
(Intercept) 1 50 62.41 2.377e-10
season 3 21 6.738 0.00233
Levene.test.P Shapiro.test.P
0.8837 0.4221

El p-valor del anova indica que hay diferencias significativas entre estaciones. Identificamos entre qué estaciones se producen dichas diferencias y estimamos sus magnitudes mediante los test a posteriori de Tukey:

## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: Tukey Contrasts
## 
## 
## Fit: lme.formula(fixed = tmic ~ season, data = lambra, random = ~1 | 
##     date, weights = varIdent(form = ~1 | season))
## 
## Linear Hypotheses:
##                      Estimate Std. Error z value Pr(>|z|)    
## summer - spring == 0  -0.4453     0.5483  -0.812  0.84835    
## autumn - spring == 0   1.6523     0.6022   2.744  0.03052 *  
## winter - spring == 0   1.2641     0.5790   2.183  0.12733    
## autumn - summer == 0   2.0976     0.5446   3.851  < 0.001 ***
## winter - summer == 0   1.7094     0.5188   3.295  0.00523 ** 
## winter - autumn == 0  -0.3882     0.5755  -0.675  0.90644    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)

Como vemos, entre el verano y la primavera no existen diferencias significativas, así como tampoco entre otoño e invierno; sí que se detectan diferencias sigificativas entre otoño/invierno y verano. También hay diferencia significativa entre otoño y primavera.

El método de Tukey nos proporciona también intervalos de confianza para estas diferencias:

## 
##   Simultaneous Confidence Intervals
## 
## Multiple Comparisons of Means: Tukey Contrasts
## 
## 
## Fit: lme.formula(fixed = tmic ~ season, data = lambra, random = ~1 | 
##     date, weights = varIdent(form = ~1 | season))
## 
## Quantile = 2.5679
## 95% family-wise confidence level
##  
## 
## Linear Hypotheses:
##                      Estimate lwr     upr    
## summer - spring == 0 -0.4453  -1.8533  0.9628
## autumn - spring == 0  1.6523   0.1059  3.1987
## winter - spring == 0  1.2641  -0.2226  2.7509
## autumn - summer == 0  2.0976   0.6990  3.4961
## winter - summer == 0  1.7094   0.3771  3.0417
## winter - autumn == 0 -0.3882  -1.8660  1.0896

 

 

Efecto combinado estación-periodo de ola

Más arriba hemos visto que tanto la altura de ola como el periodo de ola parecen tener cierta asociación con el volumen total de microplásticos, y que dicha asociación puede depender de la estación. Podemos incluir todos estos factores en un único modelo de efectos mixtos.

Hemos ajustado dos modelos; en el primer modelo incluímos interacción tanto de la altura de ola como del periodo de ola con la estación. El análisis de la varianza de dicho modelo arroja el siguiente resultado (las estaciones se han combinado en primavera-verano y otoño invierno):

  numDF denDF F-value p-value
(Intercept) 1 50 145.8 2.22e-16
periodo_ola 1 19 18.86 0.00035
season2 1 19 27.73 4.402e-05
ola 1 19 8.513 0.008829
periodo_ola:season2 1 19 13.13 0.001807
season2:ola 1 19 0.8476 0.3688

Como vemos, la interacción entre la altura de ola y la estación no resulta significativa una vez que se ha tenido en cuenta la interacción con el periodo de ola. Podemos simplificar el modelo prescindiendo de dicha interacción, con el siguiente resultado:

  numDF denDF F-value p-value
(Intercept) 1 50 142.6 2.22e-16
ola 1 20 25.52 6.096e-05
periodo_ola 1 20 9.282 0.006369
season2 1 20 19.18 0.0002901
periodo_ola:season2 1 20 12.88 0.001833

Entre ambos modelos no hay diferencias significativas, por lo que resulta preferible el segundo que es más simple:

##             Model df      AIC      BIC    logLik   Test   L.Ratio p-value
## lambra.lme      1 12 163.5480 191.3578 -69.77398                         
## lambra.lme2     2 10 160.4545 183.6294 -70.22724 1 vs 2 0.9065188  0.6356

El resultado obtenido nos indica:

  • Que existe un efecto significativo de la altura de ola (p<0.0001)

  • Que existe un efecto significativo del periodo de ola (p = 0.00637)

  • Que existe diferencia significativa en los valores medios totales de microplásticos entre primavera/verano y otroño/invierno (p = 0.00029)

  • Que el efecto del periodo de ola en otoño-invierno es distinto que en primavera-verano (p = 0.00181)

Concretamente Los parámetros estimados (y sus intervalos de confianza al 95%) en el segundo modelo son:

  lower est. upper
(Intercept) -5.48 -3.649 -1.817
ola 0.4043 0.9033 1.402
periodo_ola 0.1262 0.294 0.4619
season2winter 2.996 5.326 7.656
periodo_ola:season2winter -0.6228 -0.3989 -0.1749

Podemos interpretar estos parámetros del siguiente modo:

  • Por cada unidad que aumenta la altura de ola, el log del total de microplásticos aumenta en 0.9032519 unidades.

  • Por cada unidad que aumenta el periodo de ola, el log del total de microplásticos aumenta en 0.2940242 unidades.

  • En otoño-invierno, el log el total de microplásticos es en media 5.3261686 unidades mayor que en primavera-verano.

  • En otoño invierno, por cada unidad que aumenta el periodo de ola, el log total de microplásticos disminuye en 0.3988508- 0.2940242= 0.1048266 unidades.

 

Para validar el modelo comprobamos homoscedasticidad (test de Levene) y normalidad (test de Shapiro-Wilk) en los residuos:

Levene.test.P Shapiro.test.P
0.8764 0.3017

¿Añade algo el efecto del viento?

  numDF denDF F-value p-value
(Intercept) 1 50 96.72 2.811e-13
viento 1 20 0.00188 0.9658
periodo_ola 1 20 13.22 0.001643
season2 1 20 18.94 0.0003087
periodo_ola:season2 1 20 7.842 0.01105

Una vez que se tiene en cuenta el resto de las variables, el viento no aporta nada al modelo. Incluso si se deja el viento solo sin otras variables más que la estación:

  numDF denDF F-value p-value
(Intercept) 1 50 82.27 3.858e-12
viento 1 17 0.00173 0.9673
season 3 17 9.232 0.0007529
viento:season 3 17 1.013 0.4114

 

 

 

 

Análisis con modelos GAM

Ajustamos un modelo GAM con las variables anteriores, teniendo en cuenta que es un modelo de efectos mixtos (el día y las réplicas son efectos aleatorios, el resto de variables son efectos fijos).

## This is mgcv 1.8-17. For overview type 'help("mgcv-package")'.
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## total ~ season + s(ola) + s(periodo_ola_Invierno)
## 
## Parametric coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     2.513      1.723   1.458  0.14939    
## seasonsummer   -8.570      2.744  -3.124  0.00261 ** 
## seasonautumn   27.841      2.997   9.291 8.81e-14 ***
## seasonwinter   21.037      3.061   6.872 2.26e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                           edf Ref.df      F p-value    
## s(ola)                  1.655  1.655  8.613 0.00121 ** 
## s(periodo_ola_Invierno) 1.255  1.255 20.771 6.9e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  -10.2   
##   Scale est. = 1.5369e-05  n = 75

El ajuste anterior nos indica que existen diferencias significativas entre la primavera y el resto de las estaciones, así como que hay efecto de la altura y el periodo de ola. Los valores edf mostrados en la tabla anterior indican que la relación de ambas variables con el total de microplásticos sigue unas curvas cuyos grados son 1.65 y 1.25; téngase en cuenta que el grado 1 corresponde a una recta y el grado dos a una parábola. Los efectos de la altura y periodo de ola (durante el invierno) se recogen en las gráficas siguientes:

Ambos gráficos muestran que en el tramo en que se realiza el ajuste es muy similar a una recta, por lo que en realidad aplicar aquí un modelo GAMM no aporta gran novedad con respecto al modelo lineal mixto ajustado más arriba. De hecho, si comparamos la variabilidad explicada por un modelo con componente curva en altura y periodo de ola, con un modelo en el que el efecto de ambas variables es lineal obtenemos:

##                 Model df      AIC      BIC    logLik   Test      L.Ratio
## lambra.gam$lme      1 14 488.0814 520.5263 -230.0407                    
## lambra.gam2$lme     2 12 484.0814 511.8913 -230.0407 1 vs 2 3.384389e-08
##                 p-value
## lambra.gam$lme         
## lambra.gam2$lme       1

por lo cual podemos quedarnos con el análisis hecho más arriba en el que se ha utilizado un modelo lineal de efectos mixtos.

 

 

Por cierto, que tampoco se detecta efecto del viento con el modelo GAMM:

## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## total ~ s(viento)
## 
## Approximate significance of smooth terms:
##             edf Ref.df     F p-value
## s(viento) 2.131  2.131 3.058  0.0757

 

 

 

FAMARA

Mostramos en primer lugar un resumen de datos (media y desviación típica) por estación en Famara:

season mean sd
spring 6.248 5.217
summer 0.3438 0.3448
autumn 13.59 22.58
winter 3.5 1.801

El mismo resumen por mes:

monthYear mean sd ndays
2015-08 0.1321 0.1319 5
2015-09 1.032 1.463 2
2015-10 33.66 36.49 2
2015-11 5.494 1.958 2
2015-12 7.447 3.767 2
2016-01 4.197 2.395 2
2016-02 3.351 1.389 2
2016-03 2.285 0.8432 2
2016-04 11.56 2.504 2
2016-05 7.543 5.53 2
2016-06 1.682 0.6936 2
2016-07 0.4858 0.2871 2
2016-08 0.6927 0.3597 2
2016-09 0.6014 0.5662 1

Los datos disponibles son aproximadamente quincenales, salvo las cinco primeras observaciones (4, 8, 14, 24 y 28 de agosto de 2015, que fueron semanales.

## Time differences in days
##  [1]  4  6 10  4 17 14 16 14 17 11 17  9 28 11 11 13 15 16 16 12 13 18 12
## [24] 17 14 18 11 18 13

Considerando que aproximadamente hay quince días entre observaciones, podemos trazar gráficos de autocorrelación y autocorrelación parcial para determinar la posible estructura temporal de la serie (excluímos las observaciones segunda y cuarta para mantener el periodo aproximadamente quincenal):

Estos gráficos nos indican que no hay correlación temporal entre las sucesivas observaciones.

En esta playa, los datos son bastante “estables” salvo un outlier en ocurreel 28 de octubre y una subida entre los meses de abril-mayo de 2016. Como en el caso anterior, al abarcar el periodo de muestreo sólo un año de datos, no es posible hablar de tendencias estacionales.

Como ya se ha hecho para las Canteras, los gráficos siguientes muestran si se percibe alguna asociación del volumen total de microplásticos registrado con alguna de las otras variables. En todos los casos, suprimimos el outlier de octubre de 2015 ya que este tiene un efecto importante.

Si presentamos un ajuste global a todos los datos sin distinguir estación:

la recta así obtenida presenta una pendiente (\(\hat{\beta_1}=-0.09\))) no significativa (p=0.26316532)

Si presentamos un ajuste global a todos los datos sin distinguir estación:

la recta así obtenida presenta una pendiente (\(\hat{\beta_1}=-0.69\))) no significativa (p=0.64305141)

 

Si ajustamos una tendencia global sin distinguir entre estaciones, observamos que esta es creciente:

Aunque la pendiente de esta recta (\(\hat{\beta_1}=0.29\))) tampoco es significativa (p=0.2538)

 

 

 

Las siguientes figuras muestran la variación conjunta del total de microplásticos (promediado para cada día de observación) con la velocidad del viento, la altura de ola y el periodo de ola. Se muestran los datos tanto en escala natural:

como en escala logarítmica:

 

 

Para determinar si existe variación significativa en el volumen total de microplásticos que se registra durante las distintas estaciones en la playa de Famara aplicaremos nuevamente un modelo lineal de efectos mixtos. En primer lugar mostramos un boxplot con la cantidad total de microplásticos por estación:

Como vemos, hay bastante asimetría y aparecen algunos outliers a los que ya nos hemos referido. Para poder aplicar el modelo aplicamos la transformación:

\[y=\frac{1}{\left(total+1\right)^{3/4}}\]

a la variable total. La variable así transformada presenta un aspecto más regular:

Recordemos que habíamos detectado ya una cierta asociación con el periodo de ola. Representamos ahora todos los valores (todas las réplicas de todos los días, no sólo las medias mensuales como antes):

Como además solo hay un día con periodo de ola observado mayor que 15 y éste tiene mucha influencia sobre el periodo de ola, lo eliminamos, con lo que obtenemos la siguiente gráfica:

En lo que sigue nos quedamos solo con los días en que en Famara el periodo de ola fue menor que 15 y el total de microplásticos menor que 40 (eliminamos así dos outliers que perturban el análisis)

 

Si analizamos el efecto de la estación sin tener en cuenta ninguna otra variable obtenemos:

  numDF denDF F-value p-value
(Intercept) 1 56 497 0
season 3 24 36.74 3.937e-09
Levene.test.P Shapiro.test.P
0.9081 0.895

El p-valor del anova indica que hay diferencias significativas entre estaciones. Identificamos entre qué estaciones se producen dichas diferencias y estimamos sus magnitudes mediante los test a posteriori de Tukey:

## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: Tukey Contrasts
## 
## 
## Fit: lme.formula(fixed = tmic ~ season, data = famara, random = ~1 | 
##     date, weights = varIdent(form = ~1 | season))
## 
## Linear Hypotheses:
##                       Estimate Std. Error z value Pr(>|z|)    
## summer - spring == 0  0.510528   0.059151   8.631   <1e-06 ***
## autumn - spring == 0  0.006453   0.069445   0.093    1.000    
## winter - spring == 0  0.020201   0.075299   0.268    0.993    
## autumn - summer == 0 -0.504075   0.064390  -7.828   <1e-06 ***
## winter - summer == 0 -0.490326   0.070663  -6.939   <1e-06 ***
## winter - autumn == 0  0.013749   0.079481   0.173    0.998    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)

Como vemos, entre el otoño, invierno y primavera no existen diferencias significativas, siendo estas tres estaciones significativamente distintas del verano, que tiene una media más baja.

El método de Tukey nos proporciona también intervalos de confianza para estas diferencias:

## 
##   Simultaneous Confidence Intervals
## 
## Multiple Comparisons of Means: Tukey Contrasts
## 
## 
## Fit: lme.formula(fixed = tmic ~ season, data = famara, random = ~1 | 
##     date, weights = varIdent(form = ~1 | season))
## 
## Quantile = 2.5629
## 95% family-wise confidence level
##  
## 
## Linear Hypotheses:
##                      Estimate  lwr       upr      
## summer - spring == 0  0.510528  0.358931  0.662125
## autumn - spring == 0  0.006453 -0.171528  0.184434
## winter - spring == 0  0.020201 -0.172782  0.213185
## autumn - summer == 0 -0.504075 -0.669099 -0.339051
## winter - summer == 0 -0.490326 -0.671429 -0.309224
## winter - autumn == 0  0.013749 -0.189952  0.217449

 

 

Efecto combinado estación-periodo de ola

Más arriba hemos visto que tanto la altura de ola como el periodo de ola parecen tener cierta asociación con el volumen total de microplásticos, y que dicha asociación puede depender de la estación. Podemos incluir todos estos factores en un único modelo de efectos mixtos.

Hemos ajustado dos modelos; en el primer modelo incluímos interacción tanto de la altura de ola como del periodo de ola con la estación. El análisis de la varianza de dicho modelo arroja el siguiente resultado:

  numDF denDF F-value p-value
(Intercept) 1 56 699.2 0
periodo_ola 1 16 36.18 1.8e-05
season 3 16 40.74 1.019e-07
ola 1 16 0.7361 0.4036
periodo_ola:season 3 16 0.9687 0.4317
season:ola 3 16 1.211 0.3379

Eliminamos las interacciones, ya que no presentan significación:

  numDF denDF F-value p-value
(Intercept) 1 56 595.1 0
ola 1 22 0.4429 0.5127
periodo_ola 1 22 38.37 3.106e-06
season 3 22 32.11 3.229e-08

La altura de ola tampoco es significativa; la eliminamos del modelo:

  numDF denDF F-value p-value
(Intercept) 1 56 590.4 0
periodo_ola 1 23 30.15 1.391e-05
season 3 23 34.44 1.121e-08

Entre este modelo (que incluye solo periodo de ola y estación), y el primer modelo (con la altura de ola y todas las interacciones) no hay diferencias significativas, por lo que resulta preferible este último modelo, que es más simple:

##              Model df       AIC       BIC  logLik   Test  L.Ratio p-value
## famara.lme       1 18 -88.28959 -44.53489 62.1448                        
## famara.lme2b     2 10 -96.90581 -72.59764 58.4529 1 vs 2 7.383785  0.4958

El resultado obtenido nos indica:

  • Que existe un efecto significativo del periodo de ola (p = 0.51266) y que este efecto no varía significativamente entre estaciones (no hay interacción).

  • Que existe diferencia significativa en los valores medios totales de microplásticos entre estaciones (p<0.0001), una vez ajustado el periodo de ola.

 

Concretamente Los parámetros estimados (y sus intervalos de confianza al 95%) en el último modelo son:

  lower est. upper
(Intercept) 0.2506 0.5042 0.7577
periodo_ola -0.04972 -0.02125 0.007218
seasonsummer 0.3919 0.5011 0.6104
seasonautumn -0.09493 0.04086 0.1767
seasonwinter -0.09101 0.05284 0.1967

 

Para validar el modelo comprobamos homoscedasticidad (test de Levene) y normalidad (test de Shapiro-Wilk) en los residuos:

Levene.test.P Shapiro.test.P
0.9 0.9605

 

Ahora bien, hemos de tener en cuenta que nos encontramos ante un diseño no equilibrado (el número de datos es distinto en cada estación) y por tanto el resultado puede cambiar si los factores se incluyen el el modelo en otro orden. En el análisis anterior primero introdujimos el efecto del periodo de ola y luego el de la estación; que ambos hayan resultado significativos quiere decir que las diferencias entre estaciones subsisten aún cuando se haya ajustado por periodo de ola (esto es, a idéntico periodo de ola, la llegada de microplásticos varía de una estación a otra). Sin embargo, si el modelo se ajusta en orden inverso (primero la estación y luego el periodo de ola) al análisis de la varianza resultante es:

  numDF denDF F-value p-value
(Intercept) 1 56 590.4 0
season 3 23 43.74 1.149e-09
periodo_ola 1 23 2.242 0.1479

En otras palabras, el efecto de la estación se “come” al efecto del periodo de ola, y si se tiene en cuenta en primer lugar la variabilidad dependiente de la estación, esta es mayor que la variabilidad debida al periodo de ola. Si ajustáramos un modelo que tuviese en cuenta sólo el efecto del periodo de ola obtendríamos:

  numDF denDF F-value p-value
(Intercept) 1 56 126 6.661e-16
periodo_ola 1 26 6.575 0.01647 
##                  lower        est.       upper
## (Intercept)  0.6549290  1.14122465  1.62752029
## periodo_ola -0.1217325 -0.06792901 -0.01412552
## attr(,"label")
## [1] "Fixed effects:"

¿Añade algo el efecto del viento?

  numDF denDF F-value p-value
(Intercept) 1 56 590.6 0
periodo_ola 1 22 30.29 1.57e-05
viento 1 22 0.6697 0.4219
season 3 22 34.37 1.761e-08

El p-valor asociado al viento no resulta significativo una vez que se tiene en cuenta el efecto del periodo de ola.

 

 

 

 

Análisis con modelos GAM

Ajustamos un modelo GAM con las variables anteriores, teniendo en cuenta que es un modelo de efectos mixtos (el día y las réplicas son efectos aleatorios, el resto de variables son efectos fijos).

## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## tmic ~ season + s(ola)
## 
## Parametric coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.319756   0.043609   7.332 1.72e-10 ***
## seasonsummer 0.509610   0.056421   9.032 8.41e-14 ***
## seasonautumn 0.004754   0.067180   0.071    0.944    
## seasonwinter 0.019824   0.070960   0.279    0.781    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##        edf Ref.df     F p-value
## s(ola)   1      1 0.017   0.897
## 
## R-sq.(adj) =  -0.0695   
##   Scale est. = 0.0037128  n = 84

## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## tmic ~ season + s(periodo_ola, bs = "cr", k = 6)
## 
## Parametric coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   0.30880    0.04196   7.360 1.53e-10 ***
## seasonsummer  0.50381    0.05443   9.255 3.13e-14 ***
## seasonautumn  0.04273    0.06726   0.635    0.527    
## seasonwinter  0.05550    0.07142   0.777    0.439    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                  edf Ref.df     F p-value
## s(periodo_ola) 1.128  1.128 1.799   0.153
## 
## R-sq.(adj) =  0.176   
##   Scale est. = 0.003878  n = 84

El ajuste anterior nos indica que existen diferencias significativas entre el verano y el resto de las estaciones, así como que no hay efecto de la altura ni el periodo de ola. Los valores edf mostrados en la tabla anterior indican que la relación de ambas variables con el total de microplásticos sigue unas curvas cuyos grados son 1.65 y 1.25; téngase en cuenta que el grado 1 corresponde a una recta y el grado dos a una parábola. Los efectos de la altura y periodo de ola (durante el invierno) se recogen en las gráficas siguientes:

Ambos gráficos muestran que en el tramo en que se realiza el ajuste es muy similar a una recta, por lo que en realidad aplicar aquí un modelo GAMM no aporta gran novedad con respecto al modelo lineal mixto ajustado más arriba.

 

Tampoco se detecta efecto del viento con el modelo GAMM:

## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## total ~ s(viento)
## 
## Approximate significance of smooth terms:
##           edf Ref.df     F p-value
## s(viento)   1      1 1.493   0.225