Mercurial > repos > ecology > cb_ivr

comparison div_index.r @ 1:b67730406f1b draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload for repository https://github.com/Marie59/champ_blocs commit 0d86db7d42b608c386a54500064f5f9c9d7019a4
author ecology
date Wed, 04 Jan 2023 13:21:30 +0000
parents
children bcbad4f83dec
comparison
equal deleted inserted replaced
0:8c6142630659 1:b67730406f1b
1 # author: "Jonathan Richir"
2 # date: "01 October 2022"
3
4
5 #Rscript
6
7 ###############################
8 ## ##
9 ###############################
10
11 #####Packages : dplyr
12 # tidyr
13 # readr
14 # writexl
15 # stringr
16 # readxl
17 # tibble
18 # lubridate
19 # cowplot
20 # magrittr
21 # rmarkdown
22 library(magrittr)
23 library(dplyr)
24 #####Load arguments
25
26 args <- commandArgs(trailingOnly = TRUE)
27
28 #####Import data
29
30 if (length(args) < 1) {
31 stop("This tool needs at least 1 argument")
32 }else {
33 qecnato0 <- args[1]
34
35 }
36
37 qecnato0 <- readRDS(qecnato0)
38
39
40 # first, create vector (4) for qecb and fishing by region (same as above)
41
42 bret_egmp_basq_qecb <- c(
43 "X..algues.brunes",
44 "X..algues.rouges",
45 "X..algues.vertes",
46 "X..Cladophora",
47 "X..Lithophyllum",
48 "Nb.Littorina.obtusata",
49 "Nb.Gibbula.cineraria",
50 "Nb.Gibbula.pennanti",
51 "Nb.Gibbula.umbilicalis",
52 "Nb.Phallusia.mamillata",
53 "Nb.Tethya.aurantium",
54 "Nb.Spirobranchus.lamarckii.total",
55 "Nb.spirorbis.total",
56 "X..Eponges",
57 "X..Ascidies.Coloniales",
58 "X..Ascidies.Solitaires",
59 "X..Bryozoaires.Dresses",
60 "X..Balanes.Vivantes"
61 #, "X..Recouvrement.Sediment"
62 #, "X..Roche.Nue"
63 #, "X..Surface.Accolement"
64 )
65
66 egmp_basq_qecb <- c("Nb.Crassostrea.gigas", "Nb.Ostrea.edulis", "X..Mytilus.sp.", "X..Hermelles", "X..Hydraires")
67
68 bret_egmp_basq_fishing <- c("Nb.Cancer.pagurus..Tourteau.",
69 "Nb.Necora.puber..Etrille.",
70 "Nb.Carcinus.maenas..Crabe.vert.",
71 "Nb.Nucella.lapilus..Pourpre.",
72 "Nb.Galathea..Galathées.",
73 "Nb.Lophozozymus.incisus..ancien.Xantho.incisus.",
74 "Nb.Palaemon.sp..Crevette.bouquet.ou.crevette.rose.",
75 "Nb.Haliotis.tuberculata..Ormeau.",
76 "Nb.Littorina.littorea..Bigorneau.",
77 "Nb.Xantho.pilipes..Xanthe.poilu.",
78 "Nb.Mimachlamys.varia..Pétoncle.noir.")
79
80 egmp_basq_fishing <- c("Nb.Eriphia.verrucosa..Crabe.verruqueux.", "Nb.Octopus.vulgaris..Poulpe.", "Nb.Paracentrotus.lividus..Oursin.", "Nb.Stramonita.haemastoma..Pourpre.bouche.de.sang.")
81
82 # here I can choose to either replace spirorbis and/or spirobranchus by their log10 transformation in bret_egmp_basq_qecb vector
83 bret_egmp_basq_qecb <- replace(bret_egmp_basq_qecb, bret_egmp_basq_qecb == "Nb.spirorbis.total", "log10.Nb.spirorbis.total")
84
85 ## Diversity index
86
87
88 # adiv contains two main functions for species diversity indices: speciesdiv, which includes widely used indices such as species richness and the Shannon index, and divparam, which includes indices that have a parameter to control the importance given to rare versus abundant species in diversity measurements (Pavoine (2020) - adiv: An r package to analyse biodiversity in ecology).
89
90 # NB: just like for dissimilarity distance matrices, no sense to use the "fishing" variable lists, because either they are present for the bloc mobile and not for the bloc fixe (therefore false higher diversity for bloc mobile), either they are repeated between face supérieure and face inférieure of bloc mobile.
91
92 # function in a loop
93
94 row.names(qecnato0) <- c(paste0(qecnato0$region.site_year_month_day, "_", qecnato0$Quadrat.bis, "_", qecnato0$Type.Bloc, "_", qecnato0$Numéro.Bloc.échantillon, "_", qecnato0$Face))
95
96 # later on I can copy-paste above code to recreate variable names vector
97 #bret_egmp_basq_qecb
98 #egmp_basq_qecb
99 #Bret_EGMP.BASQ_fishing
100 #EGMP.BASQ_fishing
101
102 # remove boulder variables
103 bret_egmp_basq_qecb <- bret_egmp_basq_qecb[! bret_egmp_basq_qecb %in% c("X..Recouvrement.Sediment", "X..Roche.Nue", "X..Surface.Accolement")]
104
105 qecnato0$period <- as.character(qecnato0$period)
106 qecnato0$Face <- as.character(qecnato0$Face)
107
108 div_list <- vector("list", length(unique(qecnato0$site_year_month_day)))
109
110 for (i in c(1:nrow(qecnato0))) {
111 div_i <- dplyr::filter(qecnato0, site_year_month_day == qecnato0$site_year_month_day[i])
112
113 ifelse(unique(div_i$region) == "Bretagne", var. <- c(bret_egmp_basq_qecb), var. <- c(bret_egmp_basq_qecb, egmp_basq_qecb)) # Qu. : Why can't R's ifelse statements return vectors? => you can circumvent the problem if you assign the result inside the ifelse.
114
115 #8 remove empty row cfr: In speciesdiv(div_i[, var.]) & divparam(div_i[, var.]) : empty communities should be discarded
116 div_i <- dplyr::filter(div_i, rowSums(div_i[, var.]) > 0)
117
118 div_i_speciesdiv <- adiv::speciesdiv(div_i[, var.])
119 adiv_i_df <- data.frame(div_i_speciesdiv)
120
121 div_i_divparam <- adiv::divparam(div_i[, var.], q = c(0, 0.25, 0.5, 1, 2, 4, 8)) # When q increases, abundant species are overweighted compared to rare species, we thus expect that the evenness in species weights decreases.
122
123
124 par(mfrow = c (1, 1))
125 plot(adiv::divparam(div_i[, var.], q = 0), main = unique(div_i$site_year_month_day))
126 plot(adiv::divparam(div_i[, var.], q = 0:10), legend = FALSE, main = unique(div_i$site_year_month_day))
127
128 adiv_i_df$x <- div_i_divparam$div$`1`
129 colnames(adiv_i_df)[which(colnames(adiv_i_df) == "x")] <- paste0("Para. ISD, q = ", div_i_divparam$q[1], " (equi. richness)")
130 adiv_i_df$x <- div_i_divparam$div$`2`
131 colnames(adiv_i_df)[which(colnames(adiv_i_df) == "x")] <- paste0("Para. ISD, q = ", div_i_divparam$q[2])
132 adiv_i_df$x <- div_i_divparam$div$`3`
133 colnames(adiv_i_df)[which(colnames(adiv_i_df) == "x")] <- paste0("Para. ISD, q = ", div_i_divparam$q[3])
134 adiv_i_df$x <- div_i_divparam$div$`4`
135 colnames(adiv_i_df)[which(colnames(adiv_i_df) == "x")] <- paste0("Para. ISD, q = ", div_i_divparam$q[4])
136 adiv_i_df$x <- div_i_divparam$div$`5`
137 colnames(adiv_i_df)[which(colnames(adiv_i_df) == "x")] <- paste0("Para. ISD, q = ", div_i_divparam$q[5], " (equi. Simpson)")
138 adiv_i_df$x <- div_i_divparam$div$`6`
139 colnames(adiv_i_df)[which(colnames(adiv_i_df) == "x")] <- paste0("Para. ISD, q = ", div_i_divparam$q[6])
140 adiv_i_df$x <- div_i_divparam$div$`7`
141 colnames(adiv_i_df)[which(colnames(adiv_i_df) == "x")] <- paste0("Para. ISD, q = ", div_i_divparam$q[7])
142
143 # plot
144 par(mfrow = c(3, 2))
145 sapply(names(adiv_i_df[, c(1, 8, 2, 3, 12, 4:7, 9:11, 13:ncol(adiv_i_df))]),
146 function(cname) {
147 png(paste0(cname, "_histo.png"))
148 hist(adiv_i_df[, c(1, 8, 2, 3, 12, 4:7, 9:11, 13:ncol(adiv_i_df))][[cname]], main = "", xlab = cname, breaks = length(unique(adiv_i_df[, c(1, 8, 2, 3, 12, 4:7, 9:11, 13:ncol(adiv_i_df))][[cname]])))
149 dev.off()
150 }
151
152 )
153 par(mfrow = c(1,1))
154
155 div_list[[i]] <- adiv_i_df
156
157 rm(div_i, adiv_i_df, div_i_speciesdiv, div_i_divparam)
158
159 }
160
161 # for the error message due to richness NA data => 35 observations in 21 surveys; no reason to remove these data "...remo", was checked in the complete script.
162
163
164 div_df <- do.call("rbind", div_list)
165
166
167 # There is an issue with region.terri that are merged with no "_" ...
168
169 div_df <- tibble::add_column(div_df, rownames. = rownames(div_df), .before = "richness")
170 div_df <- tidyr::separate(div_df, rownames., into = c("region.terri.", "site_year_month_day", "Quadrat.bis", "Type.Bloc", "Numéro.Bloc.échantillon", "Face"), sep = "_")
171
172 # I therefore add these lines to solve that issue
173
174 div_df <- tibble::add_column(div_df, terri. = substring(div_df$region.terri., nchar(div_df$region.terri.)-3), .after = "region.terri.")
175
176 div_df$region.terri. <- substring(div_df$region.terri., 1, nchar(div_df$region.terri)-4)
177 div_df <- dplyr::rename(div_df, region = region.terri.)
178
179 div_df$site_year_month_day <- paste0(div_df$terri., "_", div_df$site_year_month_day)
180 div_df <- subset(div_df, select = -c(terri.))
181
182 div_df$Type.Bloc <- as.factor(div_df$Type.Bloc)
183 div_df$Face <- as.factor(div_df$Face)
184 div_df$Numéro.Bloc.échantillon <- as.integer(div_df$Numéro.Bloc.échantillon)
185
186 saveRDS(div_df, "div_df.RDS")
187 write.table(div_df, "Valeurs_stat.tabular", row.names = FALSE, quote = FALSE, sep = "\t", dec = ".", fileEncoding = "UTF-8")