Sargassum arribazon

Some graphs for Nancy

Authors
Affiliations

Diego Lizcano

WCS - WildMon

Nancy Cabanillas-Terán

ECOSUR

Published

June 4, 2024

Question

Which is the relationship between DO (oxigen) and distance to coast line on Sargassum plumes?

Sargassum plume

Sargassum plume

Set up analysis

Load libraries and set some options.

Code 
library(readxl) # Read Excel Files # Read Excel Files
library(gt) # Easily Create Presentation-Ready Display Tables
library(lubridate) # Make Dealing with Dates a Little Easier
library(stringr) # Simple, Consistent Wrappers for Common String Operations
library(readxl) # Read Excel Files # Read Excel Files
library(sf) # Simple Features for R
# library(MuMIn) # multimodel inference
library(visreg) # see trend
library(MASS) # stepAIC
# library(terra)
library(lme4) # mix model
library(sjPlot) # Data Visualization for Statistics in Social Science
library(see) # Model Visualisation Toolbox for 'easystats' and 'ggplot2'
library(performance) # Assessment of Regression Models Performance
library(mapview) # Interactive Viewing of Spatial Data in R
library(corrplot) # Visualization of a Correlation Matrix
library(DT) # A Wrapper of the JavaScript Library 'DataTables'
library(measurements) # transform DMS to lat lon
library(ggeffects) # Create Tidy Data Frames of Marginal Effects for 'ggplot' from Model Outputs
library(grateful)

library(tidyverse) # Easily Install and Load the 'Tidyverse'



options(scipen=99999)
options(max.print=99999)
options(stringsAsFactors=F)

Load data

No Map, the data has not latitude longitude

Code 
sargasso <- read_excel("C:/CodigoR/Nancy2/data/PARAMETROS_FQ6.xlsx", 
    sheet = "ALTA_BAJA", col_types = c("text", 
         "text", "text", "text", "numeric", 
         "text", "text", "text", "text", "numeric", 
         "date", "numeric", "numeric", "numeric", 
         "text", "numeric"))


# sargasso$lat <- as.numeric(conv_unit(sargasso$Lattitude, from = "deg_min_sec", to = "dec_deg"))
# sargasso$lon <- conv_unit(sargasso$Longitude , from = "deg_min_sec", to = "dec_deg") 
# sargasso$lon <- as.numeric(sargasso$lon)*(-1)

######################
# Make map #
######################
# convierte covs a puntos terra
# puntos <- vect(sargasso, geom=c("lon", "lat"), crs="EPSG:4326")
# convierte a sf

# sargasso_sf <- sf::st_as_sf(sargasso, 
#                             coords = c("lon", "lat"),
#                             crs="EPSG:4326")
# 
# mapview(sargasso_sf, zcol="DO")

Relation DO[mgL] and distance to coastline, Plume, Season, and locality

ggplot2 graph (Loess method)

Code 
# Add the regression line Loess method
ggplot(sargasso, aes(x=Distance_coast, y=DO, color=Plume)) + 
  geom_point()+
  geom_smooth() + # method=lm
  scale_color_manual (values=c("#33ccf3", "#f35a33" )) + # change color
  facet_grid(Season ~ Locality) + #, scales = "free") # no free scale
  ylab("DO mg/L")

Code 

# # Marginal density plot of x (top panel)
# xdensity <- ggplot(sargasso, aes(Distance_coast, fill=Plume)) + 
#   geom_density(alpha=.5) + 
#   scale_fill_manual(values = c('#999999','#E69F00')) + # + 
#   # theme(legend.position = "none")
#    facet_grid(Season ~ Locality, scales = "free")
# xdensity

Regresion models

A simple model (lm1)

\[ DO \sim DistanceCoast + Plume + Season + Locality \]

Code 
sargasso$Temp <- as.numeric(sargasso$Temp)

lm1 <- lm(DO ~ Distance_coast + Plume + Season + Locality, data = sargasso)
# lm1t <- lm(DO ~ Distance_coast + Temp + Plume + Season + Locality, data = sargasso)


dat <- predict_response(lm1, terms = c( "Distance_coast", "Plume", "Season", "Locality" ))
a <- plot(dat, facets = TRUE)
a # + scale_color_manual (values=c("#33ccf3", "#f35a33" ))

Code 


out1 <- check_model(lm1)
plot(out1)

Code 

summary(lm1)
#> 
#> Call:
#> lm(formula = DO ~ Distance_coast + Plume + Season + Locality, 
#>     data = sargasso)
#> 
#> Residuals:
#>     Min      1Q  Median      3Q     Max 
#> -4.3958 -0.8859 -0.0941  0.7744  3.8080 
#> 
#> Coefficients:
#>                             Estimate Std. Error t value Pr(>|t|)    
#> (Intercept)                2.4060491  0.1052386  22.863  < 2e-16 ***
#> Distance_coast             0.0026642  0.0002512  10.606  < 2e-16 ***
#> Plumewithout               1.6802324  0.1094203  15.356  < 2e-16 ***
#> SeasonNo massive arrivals  4.1045588  0.1048126  39.161  < 2e-16 ***
#> LocalityXahuayxol         -1.6919933  0.1346895 -12.562  < 2e-16 ***
#> LocalityXcalak            -0.9164167  0.1244600  -7.363 5.08e-13 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 1.345 on 697 degrees of freedom
#> Multiple R-squared:  0.7468, Adjusted R-squared:  0.745 
#> F-statistic: 411.2 on 5 and 697 DF,  p-value: < 2.2e-16
tab_model(summary(lm1))
  DO
Predictors Estimates CI p
(Intercept) 2.41 2.20 – 2.61 <0.001
Distance coast 0.00 0.00 – 0.00 <0.001
Plume [without] 1.68 1.47 – 1.90 <0.001
Season [No massive
arrivals]		 4.10 3.90 – 4.31 <0.001
Locality [Xahuayxol] -1.69 -1.96 – -1.43 <0.001
Locality [Xcalak] -0.92 -1.16 – -0.67 <0.001
Observations 703
R2 / R2 adjusted 0.747 / 0.745
Code 
#visreg(glm2, xvar = c("Distance_coast"), 
#       yvar="DO",
#       by="Plume")

A model with 2 interaction between distance and Season (lm2)

\[ DO \sim DistanceCoast * Season + Plume + Locality \]

Code 
lm2 <- lm(formula = DO ~ Distance_coast * Season +  Plume + Locality,   data = sargasso)

dat2 <- predict_response(lm2, terms = c("Distance_coast", "Plume", "Season", "Locality"))
plot(dat2, facets = TRUE)

Code 

plot_model(lm2, type = "est")

Code 
# plot_model(lm2, type = "re") # in case random effect
plot_model(lm2, type = "pred")
#> $Distance_coast

#> 
#> $Season

#> 
#> $Plume

#> 
#> $Locality

Code 


out2 <- check_model(lm2, panel=FALSE)
plot(out2)
#> $PP_CHECK

#> 
#> $NCV

#> 
#> $HOMOGENEITY

#> 
#> $OUTLIERS

#> 
#> $VIF

#> 
#> $QQ

Code 

tab_model(summary(lm2))
  DO
Predictors Estimates CI p
(Intercept) 1.68 1.46 – 1.90 <0.001
Distance coast 0.01 0.01 – 0.01 <0.001
Season [No massive
arrivals]		 5.20 4.95 – 5.46 <0.001
Plume [without] 1.61 1.42 – 1.81 <0.001
Locality [Xahuayxol] -1.52 -1.76 – -1.28 <0.001
Locality [Xcalak] -0.61 -0.84 – -0.39 <0.001
Distance coast × Season
[No massive arrivals]		 -0.01 -0.01 – -0.00 <0.001
Observations 703
R2 / R2 adjusted 0.794 / 0.792

complex to interpret

Three-Way-Interaction (lm3)

\[ DO \sim DistanceCoast * Season * Plume \]

Code 
lm3 <- lm(formula = DO ~ Distance_coast * Season * Plume ,   data = sargasso)

plot_model(lm3, type = "int")
#> [[1]]

#> 
#> [[2]]

#> 
#> [[3]]

#> 
#> [[4]]

Code 

out3 <- check_model(lm3)
plot(out3)

Code 

tab_model(summary(lm3))
  DO
Predictors Estimates CI p
(Intercept) 1.35 1.09 – 1.61 <0.001
Distance coast 0.01 0.01 – 0.01 <0.001
Season [No massive
arrivals]		 4.90 4.54 – 5.26 <0.001
Plume [without] 1.26 0.85 – 1.67 <0.001
Distance coast × Season
[No massive arrivals]		 -0.01 -0.01 – -0.01 <0.001
Distance coast × Plume
[without]		 -0.00 -0.00 – -0.00 0.003
Season [No massive
arrivals] × Plume
[without]		 0.69 0.13 – 1.25 0.015
(Distance coast × Season
[No massive arrivals]) ×
Plume [without]		 0.00 0.00 – 0.00 0.032
Observations 703
R2 / R2 adjusted 0.760 / 0.758

Mix effect model. Random effect in Locality (lm4)

\[ DO \sim DistanceCoast + Season + Plume + (1|Locality) \]

Code 
lm4 <- lmer(formula = DO ~ Distance_coast + Season + Plume + (1|Locality),   data = sargasso)


dat4 <- predict_response(lm4, terms = c("Distance_coast", "Plume", "Season", "Locality"))
plot(dat4, facets = TRUE)

Code 

plot_model(lm4, type = "est")

Code 
plot_model(lm4, type = "re")

Code 
plot_model(lm4, type = "pred")
#> $Distance_coast

#> 
#> $Season

#> 
#> $Plume

Code 


out4 <- check_model(lm4)#, panel=FALSE)
plot(out4)

Code 

tab_model((lm4))
  DO
Predictors Estimates CI p
(Intercept) 1.54 0.56 – 2.52 0.002
Distance coast 0.00 0.00 – 0.00 <0.001
Season [No massive
arrivals]		 4.10 3.90 – 4.31 <0.001
Plume [without] 1.68 1.46 – 1.89 <0.001
Random Effects
σ2 1.81
τ00 Locality 0.71
ICC 0.28
N Locality 3
Observations 703
Marginal R2 / Conditional R2 0.677 / 0.768

Mix effect model. Interaction and random effect in Locality (lm5)

\[ DO \sim DistanceCoast * Season + Plume + (1|Locality) \]

Code 
lm5 <- lmer(formula = DO ~ Distance_coast * Season + Plume + (1|Locality),   data = sargasso)


dat5 <- predict_response(lm5, terms = c("Distance_coast", "Plume", "Season", "Locality"))
plot(dat5, facets = TRUE)

Code 


plot_model(lm5, type = "est")

Code 
plot_model(lm5, type = "re")

Code 
plot_model(lm5, type = "pred")
#> $Distance_coast

#> 
#> $Season

#> 
#> $Plume

Code 

out5 <- check_model(lm5,  panel=FALSE)
plot(out5)
#> $PP_CHECK

#> 
#> $NCV

#> 
#> $HOMOGENEITY

#> 
#> $OUTLIERS

#> 
#> $VIF

#> 
#> $QQ

#> 
#> [[7]]

Code 

tab_model((lm5))
  DO
Predictors Estimates CI p
(Intercept) 0.97 0.09 – 1.86 0.032
Distance coast 0.01 0.01 – 0.01 <0.001
Season [No massive
arrivals]		 5.21 4.95 – 5.46 <0.001
Plume [without] 1.61 1.42 – 1.81 <0.001
Distance coast × Season
[No massive arrivals]		 -0.01 -0.01 – -0.00 <0.001
Random Effects
σ2 1.47
τ00 Locality 0.58
ICC 0.28
N Locality 3
Observations 703
Marginal R2 / Conditional R2 0.728 / 0.805

Which model can be the best?

Lets compare the model performance

Code 

result <- compare_performance(lm1, lm2, lm3, lm4, lm5, rank=TRUE)

DT::datatable(result)
Code 

plot(result)

Code 

test_performance(lm2, lm3, lm4, lm5)
#> Name |   Model |      BF
#> ------------------------
#> lm2  |      lm |        
#> lm3  |      lm | < 0.001
#> lm4  | lmerMod | < 0.001
#> lm5  | lmerMod |   0.012
#> Each model is compared to lm2.

Larger values indicate better model performance. Hence, points closer to the center indicate worse fit indices.

Package Citation

Code 
pkgs <- cite_packages(output = "paragraph", out.dir = ".")
knitr::kable(pkgs)
x
We used R version 4.3.2 (R Core Team 2023) and the following R packages: corrplot v. 0.92 (Wei and Simko 2021), DT v. 0.32 (Xie, Cheng, and Tan 2024), ggeffects v. 1.6.0 (Lüdecke 2018), gt v. 0.10.1 (Iannone et al. 2024), knitr v. 1.46 (Xie 2014, 2015, 2024), lme4 v. 1.1.35.3 (Bates et al. 2015), mapview v. 2.11.2 (Appelhans et al. 2023), MASS v. 7.3.60 (Venables and Ripley 2002), measurements v. 1.5.1 (Birk 2023), performance v. 0.11.0.9 (Lüdecke, Ben-Shachar, et al. 2021), rmarkdown v. 2.27 (Xie, Allaire, and Grolemund 2018; Xie, Dervieux, and Riederer 2020; Allaire et al. 2024), see v. 0.8.4.1 (Lüdecke, Patil, et al. 2021), sf v. 1.0.15 (Pebesma 2018; Pebesma and Bivand 2023), sjPlot v. 2.8.16 (Lüdecke 2024), tidyverse v. 2.0.0 (Wickham et al. 2019), visreg v. 2.7.0 (Breheny and Burchett 2017).
Code 
# pkgs

Información de la sesión en R.

Code 
sessionInfo()
#> R version 4.3.2 (2023-10-31 ucrt)
#> Platform: x86_64-w64-mingw32/x64 (64-bit)
#> Running under: Windows 10 x64 (build 19042)
#> 
#> Matrix products: default
#> 
#> 
#> locale:
#> [1] LC_COLLATE=Spanish_Colombia.utf8  LC_CTYPE=Spanish_Colombia.utf8   
#> [3] LC_MONETARY=Spanish_Colombia.utf8 LC_NUMERIC=C                     
#> [5] LC_TIME=Spanish_Colombia.utf8    
#> 
#> time zone: America/Bogota
#> tzcode source: internal
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#>  [1] forcats_1.0.0        dplyr_1.1.4          purrr_1.0.2         
#>  [4] readr_2.1.5          tidyr_1.3.1          tibble_3.2.1        
#>  [7] ggplot2_3.5.1        tidyverse_2.0.0      ggeffects_1.6.0     
#> [10] measurements_1.5.1   DT_0.32              corrplot_0.92       
#> [13] mapview_2.11.2       performance_0.11.0.9 see_0.8.4.1         
#> [16] sjPlot_2.8.16        lme4_1.1-35.3        Matrix_1.6-1.1      
#> [19] MASS_7.3-60          visreg_2.7.0         sf_1.0-15           
#> [22] stringr_1.5.1        lubridate_1.9.3      gt_0.10.1           
#> [25] readxl_1.4.3        
#> 
#> loaded via a namespace (and not attached):
#>   [1] RColorBrewer_1.1-3      rstudioapi_0.16.0       jsonlite_1.8.8         
#>   [4] datawizard_0.10.0.4     magrittr_2.0.3          TH.data_1.1-2          
#>   [7] estimability_1.5        farver_2.1.1            nloptr_2.0.3           
#>  [10] rmarkdown_2.27          vctrs_0.6.5             minqa_1.2.6            
#>  [13] effectsize_0.8.8        base64enc_0.1-3         terra_1.7-71           
#>  [16] htmltools_0.5.7         haven_2.5.4             raster_3.6-26          
#>  [19] cellranger_1.1.0        sjmisc_2.8.10           pracma_2.4.4           
#>  [22] KernSmooth_2.23-22      htmlwidgets_1.6.4       sandwich_3.1-0         
#>  [25] zoo_1.8-12              emmeans_1.10.1          uuid_1.2-0             
#>  [28] lifecycle_1.0.4         iterators_1.0.14        pkgconfig_2.0.3        
#>  [31] sjlabelled_1.2.0        R6_2.5.1                fastmap_1.1.1          
#>  [34] digest_0.6.34           colorspace_2.1-0        patchwork_1.2.0        
#>  [37] leafem_0.2.3            crosstalk_1.2.1         labeling_0.4.3         
#>  [40] fansi_1.0.6             timechange_0.3.0        mgcv_1.9-0             
#>  [43] compiler_4.3.2          proxy_0.4-27            withr_3.0.0            
#>  [46] doParallel_1.0.17       brew_1.0-10             DBI_1.2.2              
#>  [49] sjstats_0.19.0          leaflet_2.2.1           classInt_0.4-10        
#>  [52] caTools_1.18.2          tools_4.3.2             units_0.8-5            
#>  [55] qqconf_1.3.2            glue_1.7.0              satellite_1.0.5        
#>  [58] nlme_3.1-163            grid_4.3.2              generics_0.1.3         
#>  [61] qqplotr_0.0.6           gtable_0.3.4            leaflet.providers_2.0.0
#>  [64] tzdb_0.4.0              class_7.3-22            hms_1.1.3              
#>  [67] sp_2.1-3                xml2_1.3.6              utf8_1.2.4             
#>  [70] ggrepel_0.9.5           foreach_1.5.2           pillar_1.9.0           
#>  [73] robustbase_0.99-2       splines_4.3.2           lattice_0.22-5         
#>  [76] survival_3.5-7          tidyselect_1.2.1        knitr_1.46             
#>  [79] svglite_2.1.3           stats4_4.3.2            xfun_0.44              
#>  [82] leafpop_0.1.0           DEoptimR_1.1-3          stringi_1.8.3          
#>  [85] yaml_2.3.8              boot_1.3-28.1           evaluate_0.23          
#>  [88] codetools_0.2-19        twosamples_2.0.1        cli_3.6.2              
#>  [91] xtable_1.8-4            parameters_0.21.7.1     systemfonts_1.0.5      
#>  [94] pbmcapply_1.5.1         munsell_0.5.0           jquerylib_0.1.4        
#>  [97] Rcpp_1.0.12             coda_0.19-4.1           png_0.1-8              
#> [100] parallel_4.3.2          ellipsis_0.3.2          bayestestR_0.13.2.1    
#> [103] opdisDownsampling_1.0.1 bitops_1.0-7            mvtnorm_1.2-4          
#> [106] scales_1.3.0            e1071_1.7-14            insight_0.19.11.2      
#> [109] rlang_1.1.3             multcomp_1.4-25

References

Allaire, JJ, Yihui Xie, Christophe Dervieux, Jonathan McPherson, Javier Luraschi, Kevin Ushey, Aron Atkins, et al. 2024. rmarkdown: Dynamic Documents for r. https://github.com/rstudio/rmarkdown.
Appelhans, Tim, Florian Detsch, Christoph Reudenbach, and Stefan Woellauer. 2023. mapview: Interactive Viewing of Spatial Data in r. https://CRAN.R-project.org/package=mapview.
Bates, Douglas, Martin Mächler, Ben Bolker, and Steve Walker. 2015. “Fitting Linear Mixed-Effects Models Using lme4.” Journal of Statistical Software 67 (1): 1–48. https://doi.org/10.18637/jss.v067.i01.
Birk, Matthew A. 2023. measurements: Tools for Units of Measurement. https://CRAN.R-project.org/package=measurements.
Breheny, Patrick, and Woodrow Burchett. 2017. “Visualization of Regression Models Using Visreg.” The R Journal 9 (2): 56–71.
Iannone, Richard, Joe Cheng, Barret Schloerke, Ellis Hughes, Alexandra Lauer, and JooYoung Seo. 2024. gt: Easily Create Presentation-Ready Display Tables. https://CRAN.R-project.org/package=gt.
Lüdecke, Daniel. 2018. “ggeffects: Tidy Data Frames of Marginal Effects from Regression Models.” Journal of Open Source Software 3 (26): 772. https://doi.org/10.21105/joss.00772.
———. 2024. sjPlot: Data Visualization for Statistics in Social Science. https://CRAN.R-project.org/package=sjPlot.
Lüdecke, Daniel, Mattan S. Ben-Shachar, Indrajeet Patil, Philip Waggoner, and Dominique Makowski. 2021. “performance: An R Package for Assessment, Comparison and Testing of Statistical Models.” Journal of Open Source Software 6 (60): 3139. https://doi.org/10.21105/joss.03139.
Lüdecke, Daniel, Indrajeet Patil, Mattan S. Ben-Shachar, Brenton M. Wiernik, Philip Waggoner, and Dominique Makowski. 2021. “see: An R Package for Visualizing Statistical Models.” Journal of Open Source Software 6 (64): 3393. https://doi.org/10.21105/joss.03393.
Pebesma, Edzer. 2018. “Simple Features for R: Standardized Support for Spatial Vector Data.” The R Journal 10 (1): 439–46. https://doi.org/10.32614/RJ-2018-009.
Pebesma, Edzer, and Roger Bivand. 2023. Spatial Data Science: With applications in R. Chapman and Hall/CRC. https://doi.org/10.1201/9780429459016.
R Core Team. 2023. R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing. https://www.R-project.org/.
Venables, W. N., and B. D. Ripley. 2002. Modern Applied Statistics with s. Fourth. New York: Springer. https://www.stats.ox.ac.uk/pub/MASS4/.
Wei, Taiyun, and Viliam Simko. 2021. R Package “corrplot”: Visualization of a Correlation Matrix. https://github.com/taiyun/corrplot.
Wickham, Hadley, Mara Averick, Jennifer Bryan, Winston Chang, Lucy D’Agostino McGowan, Romain François, Garrett Grolemund, et al. 2019. “Welcome to the tidyverse.” Journal of Open Source Software 4 (43): 1686. https://doi.org/10.21105/joss.01686.
Xie, Yihui. 2014. “knitr: A Comprehensive Tool for Reproducible Research in R.” In Implementing Reproducible Computational Research, edited by Victoria Stodden, Friedrich Leisch, and Roger D. Peng. Chapman; Hall/CRC.
———. 2015. Dynamic Documents with R and Knitr. 2nd ed. Boca Raton, Florida: Chapman; Hall/CRC. https://yihui.org/knitr/.
———. 2024. knitr: A General-Purpose Package for Dynamic Report Generation in r. https://yihui.org/knitr/.
Xie, Yihui, J. J. Allaire, and Garrett Grolemund. 2018. R Markdown: The Definitive Guide. Boca Raton, Florida: Chapman; Hall/CRC. https://bookdown.org/yihui/rmarkdown.
Xie, Yihui, Joe Cheng, and Xianying Tan. 2024. DT: A Wrapper of the JavaScript Library “DataTables”. https://CRAN.R-project.org/package=DT.
Xie, Yihui, Christophe Dervieux, and Emily Riederer. 2020. R Markdown Cookbook. Boca Raton, Florida: Chapman; Hall/CRC. https://bookdown.org/yihui/rmarkdown-cookbook.

Reuse

Citation

BibTeX citation:
@online{lizcano,
  author = {Lizcano, Diego and Cabanillas-Terán, Nancy},
  title = {*Sargassum* Arribazon},
  date = {},
  url = {https://dlizcano.github.io/Sargasso/Sargasso.html},
  langid = {en}
}
For attribution, please cite this work as:
Lizcano, Diego, and Nancy Cabanillas-Terán. n.d. “*Sargassum* Arribazon.” https://dlizcano.github.io/Sargasso/Sargasso.html.