Mercurial > repos > ecology > stoc_mainglm

comparison FunctTrendSTOCGalaxy.r @ 0:df3ce23d0d23 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
"planemo upload for repository https://github.com/Alanamosse/Galaxy-E/tree/stoctool/tools/stoc commit f82f897ab22464de40c878e17616333855814e25"
author ecology
date Thu, 02 Apr 2020 03:32:56 -0400
parents
children a961be0a964f
comparison
equal deleted inserted replaced
-1:000000000000 0:df3ce23d0d23
1 #!/usr/bin/env Rscript
2
3
4 ##################################################################################################################################
5 ############## FUNCTION TO CALCULATE AND PLOT EVOLUTION OF SPECIES POPULATION function:main.glm ##############################
6 ##################################################################################################################################
7
8 #### Based on Romain Lorrillière R script
9 #### Modified by Alan Amosse and Benjamin Yguel for integrating within Galaxy-E
10
11 ##### workes with the R version 3.5.1 (2018-07-02)
12 ##### Package used with the version:
13 #suppressMessages(library(lme4)) version 1.1.18.1
14 #suppressMessages(library(ggplot2)) version 3.0.0
15 #suppressMessages(library(speedglm)) version 0.3.2
16 #suppressMessages(library(arm)) version 1.10.1
17 #suppressMessages(library(reshape)) version 0.8.8
18 #suppressMessages(library(data.table)) version 1.12.0
19 #suppressMessages(library(reshape2)) version 1.4.3
20
21
22
23 ######################################### debut de la fonction makeTableAnalyse / stard of the function makeTableAnalyse
24 ## mise en colonne des especes et rajout de zero mais sur la base des carrés selectionné sans l'import / Species are placed in separated columns and addition of zero on plots where at least one selected species is present
25
26 makeTableAnalyse <- function(data) {
27 tab <- reshape(data
28 ,v.names="abond"
29 ,idvar=c("carre","annee")
30 ,timevar="espece"
31 ,direction="wide")
32 tab[is.na(tab)] <- 0 ###### remplace les na par des 0 / replace NAs by 0
33
34 colnames(tab) <- sub("abond.","",colnames(tab))### remplace le premier pattern "abond." par le second "" / replace the column names "abond." by ""
35 return(tab)
36 }
37
38 ######################################### fin de la fonction makeTableAnalyse / end of the function makeTableAnalyse
39
40
41
42
43
44 ############################################# les fonctions qui filtrent les données pas suffisantes pour analyses fiables / The filtering functions removing species with not enough data to perform accurate analyses
45
46 filter_absent_species<-function(tab){
47 ##################### Filtre les espèces jamais présentes (abondance=0) / Filter of species with 0 abundance
48 ################################################################################# PARTIE POTENTIELLEMENT ISOLABLE ET INSERABLE AVANT LA BOUCLE = permet de gagner du temps sur la boucle car supprime sps pas vu, donc pas repris par la boucle
49
50 ## Fait la somme des abondances totales par espèce / calculate the sum of all abundance per species
51 if(ncol(tab)==3) {
52 tabSum <- sum(tab[,3])## cas d'une seule especes (problème de format et manip un peu differente) / when selecting only one species, use a different method
53 names(tabSum) <- colnames(tab)[3]
54 } else { ## cas de plusieurs espèce/ when selecting more than one species
55 tabSum <- colSums(tab[,-(1:2)])
56 }
57 ## colNull= espece(s) toujours absente /species with 0 total abundance
58 colNull <- names(which(tabSum==0))
59 ## colconserve= espece(s) au moins presente 1 fois/ species at least with 1 presence
60 colConserve <- names(which(tabSum>0))
61 ## Affichage des espèces rejetees / show species eliminated for the analyses
62 if(length(colNull)>0){
63 cat("\n",length(colNull)," Species removed from the analysis, abundance is always 0.\n\n",sep="") #Espèces enlevées de l'analyse car abondance toujours égale a 0\n\n",sep="")
64 #tabNull <- data.frame(Code_espece = colNull, nom_espece = tabsp[colNull,"nom"])
65 #cat("\n\n",sep="")
66 tab <- tab[,c("carre","annee",colConserve)]
67 }
68 ################################################################################ FIN DE LA PARTIE ISOLABLE
69 return(tab)
70 }
71
72
73
74
75 ###################### Filtre les especes trop rare pour avoir des analyses robustes i.e. espèce non presente la 1ère année, avec plus de 3 ans consecutif sans données et moins de 3 ans consécutif avec données
76 ###################### Filter too rare species for accurate analysis i.e. species absent the first year, with more than 3 consecutive years with 0 abundance, or with less than 3 consecutive years with presence
77
78 ###
79 filter_rare_species<-function(tab){
80 exclude_threshold <- NULL
81 ## calcul et filtre pour chaque (colonne) espece / measure and filter for each species
82 for(i in 3:ncol(tab)) {
83 ## v =abondance par annee / v= abundance per year
84 v <- tapply(tab[,i],tab$annee,sum) ####################
85 ## v0 =presence(1) abscence(0) per year
86 v0 <- ifelse(v>0,1,0) #####
87 tx <- paste(v0,collapse="") #### colle les 0 et 1 / stick the 0 and 1
88
89 p <- unlist(strsplit(tx,"0"))#### Enleve les 0, ce qui séparent les sequences de "1", les sequences de "1" = nbre d'années consécutives avec data / remove 0, splitting sequences of "1" which correspond to consecutve year with data (e.g. 111 = 3 years consecutive years with data)
90 p <- p[p!=""] #### ne garde pas les partie sans 1 ou 0 dans les sequences
91 ## gsSup0 = plus grande serie temporelle de presence =calcul du nbre de 1 consécutif max / calcul of the biggest temporal series which corresponds to the maximum number of consecutive "1"
92 gsSup0 <- max(nchar(p))####
93 ## gsInf0 plus grande serie temporelle d'absccence ou sans données = enlève les 1 séparant sequence de 0 qui correspondent au nbre d'année consecutive sans données / calcul of the biggest temporal series without data which corresponds to max numbzer fo consecutive "0"
94 gsInf0 <- max(nchar(unlist(strsplit(tx,"1")))) ####
95 ## y0is0 absence la premiere annee
96 y0is0 <- v0[1]==0 #### True ou false pour presence de "0"(=pas de données) dans la 1ère année / look if the first year of the time sequence analyzed has no data
97 ## seuil d'exclusion / exclusion threshold
98 exclude_threshold <- c(exclude_threshold,as.vector(ifelse( y0is0 | gsInf0 > 3 | gsSup0 < 3 ,"exclu","bon"))) ############## exclu sps absente la 1ère année, avec plus de 3 ans consécutifs sans données, et avec moins de 3 années consécutives sans données / indicate if the max consecutive year with data and without data, as well as whether the first year of the time sequence analyzed has data
99 }
100 names(exclude_threshold) <- colnames(tab)[3:ncol(tab)]
101
102 ## colonnes conservees avec assez de données / Column with enough data
103 colConserve <- names(exclude_threshold)[exclude_threshold=="bon"]
104
105
106 ## colonnes supprimees / Column that will conserved
107 colSupr <- names(exclude_threshold)[exclude_threshold=="exclu"]
108 tabCLEAN <- tab[,c("carre","annee",colConserve)] #### Garde les sps à conserver / select only species with enough data
109 lfiltre <- list(tabCLEAN=tabCLEAN,colConserve=colConserve,colSupr=colSupr)
110
111 #################################################################################
112
113 ## colConserve espece conservees / extract species that will be kept to print them
114 colConserve <- lfiltre$colConserve
115 ## colsupr espece trop rare et donc supprimée de l'analyse / extract species that will be deleted to print them
116 colSupr <- lfiltre$colSupr
117 ## affichage des especes retirer de l'analyse / print species that will be deleted
118 if(length(colSupr)>0){
119 cat("\n",length(colSupr)," Rare species removed from the analysis.\n\n",sep="")
120 #tabSupr <- subset(tabsp,espece %in% colSupr ,select=c("espece","nom"))
121 #tabSupr <- tabSupr[order(tabSupr$espece),]
122 #cat("\n\n",sep="")
123
124 }
125 if(length(colConserve)==0) {
126 mess <- "No species available to calculate abundance variation in this dataset."
127 stop(mess)
128 }
129
130 tabCLEAN <- lfiltre$tabCLEAN
131
132 #### MARCHE PAS NE SAIT PAS PQUOI
133 tabCLEAN <- melt(tabCLEAN, id.vars=c("carre", "annee")) #### remet le format de base :le nom d'espèce et abondance dans des colonnes séparées / back to the first format of the file: species name and abundance in separated column
134
135 colnames(tabCLEAN)[3:4] <- c("espece","abond")
136 tabCLEAN$annee <- as.numeric(as.character(tabCLEAN$annee))
137 ################################################################################
138 return(tabCLEAN)
139 }
140
141 ####################################################################################################################### fin des 2 fonctions de filtre des données / end of the two function to filter the data
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168 ############################################################################################ debut de la Function main.glm / start of the function main.glm
169
170 main.glm <- function(id="france",donneesAll=dataCLEAN,assessIC= TRUE,tabsp=tabsp,annees=annees,figure=TRUE,description=TRUE,tendanceSurFigure=TRUE, ###### declaration des arguments listSp=sp était avant declaré avant la fonction mais il me semble que ca marche aussi comme cela
171 seuilOccu=14,seuilAbond=NA) {
172
173
174
175 filesaveAn <- paste("Output/",id,"/variationsAnnuellesEspece_",id,".tabular", ##### Nom du dossier ET fichier de sortie des resultats par année / name of the output file with results for each years
176 sep = "")
177 filesaveTrend <- paste("Output/",id,"/tendanceGlobalEspece_",id,".tabular", ##### Nom du dossier ET fichier de sortie des resultats pour la période "annee" complete / name of the output file with the results for the period
178 sep = "")
179 fileSaveGLMs <- paste("Output/",id,"/listGLM_",id,sep = "") ##### Nom du dossier ET fichier de sortie des modèles lineaire generalisés / name of the output file of the generlized linear models
180
181
182
183
184 seuilSignif <- 0.05 ## seuil de significativite / significancy threshold
185
186
187 rownames(tabsp) <- tabsp$espece ## change nom des lignes de tabsp (table de reference des especes)
188
189
190 ##vpan vecteur des panels de la figure ###### POUR FAIRE LES GRAPHIQUES
191 vpan <- c("Variation abondance")
192 if(description) vpan <- c(vpan,"Occurrences","Abondances brutes")
193
194
195 ## specifications des variables temporelles necesaires pour les analyses / specification of temporal variable necessary for the analyses
196 annee <- sort(unique(donneesAll$annee))
197 nbans <- length(annee)
198 pasdetemps <- nbans-1
199 firstY <- min(annee)
200 lastY <- max(annee)
201
202
203
204
205
206 ## Ordre de traitement des especes ### order of species to be analyzed
207 spOrdre <- aggregate(abond~espece,data=donneesAll,sum) #### calcul les sommes des abondances pour ordonner / calculate the sum for the ordination
208 spOrdre <- merge(spOrdre,tabsp,by="espece") #### rajoute la colonne avec les abondances totales par espece / add a new column with the sum
209
210 spOrdre <- spOrdre[order(as.numeric(spOrdre$indicateur),spOrdre$abond,decreasing = TRUE),] #### mets les especes plus abondantes en premiers (plus long pour faire tourner le modèle) / order the species by abundance, the most abundant species being the less fast analysis
211
212
213 listSp <- spOrdre$espece
214 i <- 0
215 nbSp <- length(listSp)
216 # browser()
217 ## analyse par espece
218 ### browser()
219 ## affichage des especes conservees pour l'analyse ### PAS SUR QUE CE SOIT ENCORE UTILE
220 cat("\n",nbSp," Espèces conservées pour l'analyse\n\n",sep="")
221 rownames(tabsp) <- tabsp$espece
222 print(tabsp[,1:2])
223 #tabCons <- data.frame(Code_espece = listSp, nom_espece = tabsp[as.character(listSp),"nom"])
224 #print(tabCons)
225 cat("\n\n",sep="")
226 flush.console()
227
228
229 ## initialisation de la liste de sauvegarde
230
231
232 ##browser()
233
234 for (sp in listSp) { ######## Boucle pour analyse par espèce / loop for the analysis by species
235
236
237 i <- i + 1
238
239 d <- subset(donneesAll,espece==sp) ## d data pour l'espece en court / cut the data keeping only the i species
240
241 #nomSp <- as.character(tabsp[sp,"nom"]) ## info sp
242 nomSp <- tabsp$nom[which(tabsp$espece==sp)] ## info sp
243 cat("\n(",i,"/",nbSp,") ",sp," | ", nomSp,"\n",sep="")
244 flush.console()
245
246 #indic <- tabsp[sp,"indicateur"] ## indic :espece utilisee pour le calcul des indicateurs par groupe de specialisation / list the species used as species indicators by trophic specialization
247 indic <- tabsp$indicateur[which(tabsp$espece==sp)] ## indic :espece utilisee pour le calcul des indicateurs par groupe de specialisation / list the species used as species indicators by trophic specialization
248 nb_carre = tapply(rep(1,nrow(d)),d$annee,sum) ## nb_carre nombre de carre suivie par annee / number of plots per year
249
250 nb_carre_presence = tapply(ifelse(d$abond>0,1,0),d$annee,sum) ## nb_carre_presence nombre de carre de presence par annee / number the plots where the species were observed
251
252 tab2 <- data.frame(annee=rep(annee,2),val=c(nb_carre,nb_carre_presence),LL = NA,UL=NA, ## tab2 table de resultat d'analyse / data.frame of the analyses results
253 catPoint=NA,pval=NA,
254 courbe=rep(c("carre","presence"),each=length(annee)),panel=vpan[2])
255 tab2$catPoint <- ifelse(tab2$val == 0,"0",ifelse(tab2$val < seuilOccu,"infSeuil",NA))
256
257 abond <- tapply(d$abond,d$annee,sum) ## abond abondance par annee / abundance per year
258
259 tab3 <- data.frame(annee=annee,val=abond,LL = NA,UL=NA,catPoint=NA,pval=NA,courbe=vpan[3],panel=vpan[3]) ## table pour la figure / data.frame made to realize the graphical outputs
260 tab3$catPoint <- ifelse(tab3$val == 0,"0",ifelse(tab3$val < seuilAbond,"infSeuil",NA))
261
262 ## GLM pour calcul des tendances annuelles de l'evolution des populations / GLM to measure annual tendency of population evolution
263 formule <- as.formula("abond~as.factor(carre)+as.factor(annee)") #### specification du modèle = log lineaire / specifying the model = log linear
264 if(assessIC) {##### OPTION A RENTRER AU DEBUT PEUT ËTRE A METTRE DANS LES ARGUMENTS SI LAISSE LE CHOIX SINON L ARG PAR DEFAUT LORS DE LA DECLARATION DE LA FONCTION
265 glm1 <- glm(formule,data=d,family=quasipoisson) ##### fit model lineaire general avec intervalle de confiance disponible / fit linear and generalized model with confidence intervalle available
266 } else {
267 glm1 <- try(speedglm(formule,data=d,family=quasipoisson())) ##### fit modele lineaire et generaux pour les gros jeux de données / fit of linear and generalized model for large-medium dataset
268 if(class(glm1)[1]=="try-error")
269 glm1 <- glm(formule,data=d,family=quasipoisson) ##### comprends pas mais je pense que c'est speedglm qui marche pas avec toutes les données
270 }
271 sglm1 <- summary(glm1) #### sortie du modele / output of the model
272 sglm1 <- coefficients(sglm1) ### coefficient regression de chaque variable avec les résultats des tests statistiques / regression coefficient of each predictive variables with results of the statistical tests
273 sglm1 <- tail(sglm1,pasdetemps) #### recupére les derniers elements du modèle avec la taille de l'objet "pasdetemps" car le nombre de coef = nbre d'année et pas les coefficient de regression de la variable carre / retrieve only the coefficient regression of the variable year
274 coefan <- as.numeric(as.character(sglm1[,1]))#### coefficient de regression de la variable année (1 pour chaque année)
275
276 coefannee <- c(1,exp(coefan))## coefannee vecteur des variation d'abondance par annee avec transformation inverse du log :exp() / regression coefficient of the year back transformed from log(abundance) : exp()
277
278 erreuran <- as.numeric(as.character(sglm1[,2])) #### erreur standard sur le coefficient de regression de la variable annee / standard error on the regression coefficient of the year
279 erreurannee1 <- c(0,erreuran*exp(coefan))## erreur standard par année / the standard error per year ###### LA J AI UN DOUTE NORMALEMENT INTERVAL DE CONF C CI_lower <- coefficients(lin_mod)[2] - 1.96*summary(lin_mod)$coefficients[2,2]
280 ####CI_upper <- coefficients(lin_mod)[2] + 1.96*summary(lin_mod)$coefficients[2,2]
281
282 pval <- c(1,as.numeric(as.character(sglm1[,4])))###### p value
283
284 ## calcul des intervalle de confiance avec methode de bootstrap pour simuler des coef de regress sur lequel intervalle de conf sont mesurés/ calcul of the confidence interval using bootstrap method to simulate set regression coefficients and s.e.with uncertainty POURQUOI PAS UTILISE confint.glm() ou boot() ou ci.boot()
285
286 if(assessIC) {
287 glm1.sim <- sim(glm1)
288 ic_inf_sim <- c(1,exp(tail(apply(coef(glm1.sim), 2, quantile,.025),pasdetemps)))
289 ic_sup_sim <- c(1,exp(tail(apply(coef(glm1.sim), 2, quantile,.975),pasdetemps)))
290 } else {
291 ic_inf_sim <- NA
292 ic_sup_sim <- NA
293
294 }
295
296
297
298 tab1 <- data.frame(annee,val=coefannee, ## tab1 table pour la realisation des figures / table for the graphical outputs ### 2EME POUR GRAPH ici ce sont le coef de regress annee en fonction des annéés alors que tab3 c'est les abondance en fct des années et tab2 nombre de carré total et avec presence
299 LL=ic_inf_sim,UL=ic_sup_sim,
300 catPoint=ifelse(pval<seuilSignif,"significatif",NA),pval,
301 courbe=vpan[1],
302 panel=vpan[1])
303 ## netoyage des intervalle de confiance mal estimés et qd donnees pas suffisantes pour calcul d'IC /cleaning of wrong or biaised measures of the confidence interval
304 if(assessIC) {
305 tab1$UL <- ifelse( nb_carre_presence==0,NA,tab1$UL)
306 tab1$UL <- ifelse(tab1$UL == Inf, NA,tab1$UL)
307 tab1$UL <- ifelse(tab1$UL > 1.000000e+20, NA,tab1$UL)
308 tab1$UL[1] <- 1
309 tab1$val <- ifelse(tab1$val > 1.000000e+20,1.000000e+20,tab1$val)
310 }
311 ## indice de surdispersion / overdispersion index
312 ## browser()
313 if(assessIC) dispAn <- glm1$deviance/glm1$null.deviance else dispAn <- glm1$deviance/glm1$nulldev
314
315
316 ## tabAn table de sauvegarde des resultats par année / table of the results per year ###### reprends bcp de tabl DIFFERENCE AVEC tab2 c les abondances relatives, alors que nb de carre, nb de carre presnce, p val sont aussi ds tab2
317 tabAn <- data.frame(id,code_espece=sp, nom_espece = nomSp,indicateur = indic,annee = tab1$annee,
318 abondance_relative=round(tab1$val,3),
319 IC_inferieur = round(tab1$LL,3), IC_superieur = round(tab1$UL,3),
320 erreur_standard = round(erreurannee1,4),
321 p_value = round(tab1$pval,3),significatif = !is.na(tab1$catPoint),
322 nb_carre,nb_carre_presence,abondance=abond)
323
324 ## GLM pour calcul des tendance generale sur la periode avec modele log lineaire / GLM to measure the tendency of population evolution on the studied period with log linear model
325 formule <- as.formula(paste("abond~ as.factor(carre) + annee",sep="")) ###
326 # browser()
327
328
329 if(assessIC) {
330 md2 <- glm(formule,data=d,family=quasipoisson) }
331 else {
332 md2 <- try(speedglm(formule,data=d,family=quasipoisson()),silent=TRUE)
333
334 if(class(md2)[1]=="try-error")
335 md2 <- glm(formule,data=d,family=quasipoisson)
336 }
337
338
339 smd2 <- summary(md2) #### sortie du modele / output of the model
340 smd2 <- coefficients(smd2) ### coefficient regression de chaque variable avec les résultats des tests statistiques / regression coefficient of each predictive variables with results of the statistical tests
341 smd2 <- tail(smd2,1) ### coefficient regression de variable annee avec les résultats des tests statistiques / regression coefficient of the variable year with results of the statistical tests
342
343
344 coefan <- as.numeric(as.character(smd2[,1])) ## tendences sur la periode = coefficient regression de variable annee / tendency of population evolution on the studied period = regression coefficient of the variable year
345 trend <- round(exp(coefan),3)
346
347 pourcentage <- round((exp(coefan*pasdetemps)-1)*100,2) ## pourcentage de variation sur la periode / percentage of population variation on the studied period
348 pval <- as.numeric(as.character(smd2[,4]))
349
350 erreuran <- as.numeric(as.character(smd2[,2])) #### récuperer l'erreur standard / retrieve the error
351 ## erreur standard
352 erreurannee2 <- erreuran*exp(coefan)
353
354
355 ## calcul des intervalle de confiance avec methode de bootstrap pour simuler des coef de regress sur lequel intervalle de conf sont mesurés/ calculating the confidence interval based on bootstrap method to simulate set regression coefficients and s.e.with uncertainty
356 LL <- NA
357 UL <- NA
358 if(assessIC) {
359 md2.sim <- sim(md2)
360 LL <- round(exp(tail(apply(coef(md2.sim), 2, quantile,.025),1)),3)
361 UL <- round(exp(tail(apply(coef(md2.sim), 2, quantile,.975),1)),3)
362 } else {
363 LL <- NA
364 UL <- NA
365 }
366
367 ## tab1t table utile pour la realisation des figures / table used for the figures
368 tab1t <- data.frame(Est=trend,
369 LL , UL,
370 pourcent=pourcentage,signif=pval<seuilSignif,pval)
371
372
373 trendsignif <- tab1t$signif
374 pourcent <- round((exp(coefan*pasdetemps)-1)*100,3)
375 ## mesure de la surdispersion / overdispersion measurment
376
377 if(assessIC) dispTrend <- md2$deviance/md2$null.deviance else dispTrend <- md2$deviance/md2$nulldev
378
379
380
381 ## classement en categorie incertain /classifying wrong or not reliable results
382 # browser()
383 if(assessIC) {
384 if(dispTrend > 2 | dispAn > 2 | median( nb_carre_presence)<seuilOccu) catIncert <- "Incertain" else catIncert <-"bon" ##### en fonction de l'indice de surdispersion et presence < à seuil occurence / based on the overdispersion index and the presence on a minimum number of plots
385 vecLib <- NULL
386 if(dispTrend > 2 | dispAn > 2 | median( nb_carre_presence)<seuilOccu) {
387 if(median( nb_carre_presence)<seuilOccu) {
388 vecLib <- c(vecLib,"espece trop rare")
389 }
390 if(dispTrend > 2 | dispAn > 2) {
391 vecLib <- c(vecLib,"deviance")
392 }
393 }
394 raisonIncert <- paste(vecLib,collapse=" et ")
395 } else {
396 catIncert <- NA
397 raisonIncert <- NA
398 }
399
400
401
402 ## affectation des tendence EBCC / retrieve the trend of population evolution on the studied period
403 catEBCC <- NA
404 if(assessIC) catEBCC <- affectCatEBCC(trend = as.vector(trend),pVal = pval,ICinf=as.vector(LL),ICsup=as.vector(UL)) else catEBCC <- NA
405 ## table complete de resultats pour la periode etudiée / complete table with results for the studied period
406 # browser()
407 tabTrend <- data.frame(
408 id,code_espece=sp,nom_espece = nomSp,indicateur = indic,
409 nombre_annees = pasdetemps,premiere_annee = firstY,derniere_annee = lastY,
410 tendance = as.vector(trend) , IC_inferieur=as.vector(LL) , IC_superieur = as.vector(UL),pourcentage_variation=as.vector(pourcent),
411 erreur_standard = as.vector(round(erreurannee2,4)), p_value = round(pval,3),
412 significatif = trendsignif,categorie_tendance_EBCC=catEBCC,mediane_occurrence=median( nb_carre_presence) ,
413 valide = catIncert,raison_incertitude = raisonIncert)
414
415
416 if(assessIC) listGLMsp <- list(list(glm1,glm1.sim,md2,md2.sim)) else listGLMsp <- list(list(glm1,md2))
417 names(listGLMsp)[[1]] <-sp
418 fileSaveGLMsp <- paste(fileSaveGLMs,"_",sp,".Rdata",sep="")
419
420 save(listGLMsp,file=fileSaveGLMsp)
421 cat("--->",fileSaveGLMsp,"\n")
422 flush.console()
423
424 if(sp==listSp[1]) {
425 glmAn <- tabAn
426 glmTrend <- tabTrend
427 } else {
428 glmAn <- rbind(glmAn,tabAn)
429 glmTrend <- rbind(glmTrend,tabTrend)
430 }
431 ## les figures
432 if(figure) {
433 ## table complete pour la figure en panel par ggplot2
434 ## table pour graphe en panel par ggplot2
435 if(description) dgg <- rbind(tab1,tab2,tab3) else dgg <- tab1
436 ## les figures
437
438 ggplot.espece(dgg,tab1t,id,serie=NULL,sp,valide=catIncert,nomSp,description,tendanceSurFigure,seuilOccu=14,vpan = vpan,assessIC=assessIC)
439
440 }
441
442
443
444
445 }
446
447 write.table(glmAn,filesaveAn,row.names=FALSE,quote=FALSE,sep="\t",dec=".",fileEncoding="UTF-8")
448 cat("--->",filesaveAn,"\n")
449 write.table(glmTrend,filesaveTrend,row.names=FALSE,quote=FALSE,sep="\t",dec=".",fileEncoding="UTF-8")
450 cat("--->",filesaveTrend,"\n")
451
452
453 flush.console()
454
455
456
457 }
458 ########################################################################################################## Fin de la fonction main.glm / end of the function main.glm
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474 ########################################################################################################### fonction appelée par main.glm renvoyant la categorie European Bird Census Council en fonction des resultats des modèles / function called by main.glm to classify results depending on the quality of the data and analyses
475 ## renvoie la categorie EBCC de la tendance en fonction
476 ## trend l'estimateur de la tendance / estimation of the trends
477 ## pVal la p value
478 ## ICinf ICsup l intervalle de confiance a 95 pourcent
479 affectCatEBCC <- function(trend,pVal,ICinf,ICsup){
480 catEBCC <- ifelse(pVal>0.05,
481 ifelse(ICinf < 0.95 | ICsup > 1.05,"Incertain","Stable"),
482 ifelse(trend<1,
483 ifelse(ICsup<0.95,"Fort declin","Declin moderee"),
484 ifelse(ICinf>1.05,"Forte augmentation","Augmentation modere")))
485 return(catEBCC)
486 }
487
488 ############################################################################################################ fin de la fonction renvoyant la categorie EBCC / end of the function main.glm
489
490
491
492
493
494
495
496
497 ############################################################################################################ fonction graphique appelée par main.glm / function called by main.glm for graphical output
498 ggplot.espece <- function(dgg,tab1t,id,serie=NULL,sp,valide,nomSp=NULL,description=TRUE,
499 tendanceSurFigure=TRUE,seuilOccu=14, vpan,assessIC=TRUE) {
500
501 # serie=NULL;nomSp=NULL;description=TRUE;valide=catIncert
502 # tendanceSurFigure=TRUE;seuilOccu=14
503 require(ggplot2)
504
505 figname<- paste("Output/",id,"/",ifelse(valide=="Incertain","Incertain/",""),
506 sp,"_",id,serie, ".png",
507 sep = "")
508 ## coordonnee des ligne horizontal de seuil pour les abondances et les occurences
509 hline.data1 <- data.frame(z = c(1), panel = c(vpan[1]),couleur = "variation abondance",type="variation abondance")
510 hline.data2 <- data.frame(z = c(0,seuilOccu), panel = c(vpan[2],vpan[2]),couleur = "seuil",type="seuil")
511 hline.data3 <- data.frame(z = 0, panel = vpan[3] ,couleur = "seuil",type="seuil")
512 hline.data <- rbind(hline.data1,hline.data2,hline.data3)
513 titre <- paste(nomSp)#,"\n",min(annee)," - ",max(annee),sep="")
514
515 ## texte de la tendance / text for the population evolution trend
516 tab1 <- subset(dgg,panel =="Variation abondance")
517 pasdetemps <- max(dgg$annee) - min(dgg$annee) + 1
518 if(assessIC){
519 txtPente1 <- paste(tab1t$Est,
520 ifelse(tab1t$signif," *","")," [",tab1t$LL," , ",tab1t$UL,"]",
521 ifelse(tab1t$signif,paste("\n",ifelse(tab1t$pourcent>0,"+ ","- "),
522 abs(tab1t$pourcent)," % en ",pasdetemps," ans",sep=""),""),sep="")
523 }else{
524 txtPente1 <- ifelse(tab1t$signif,paste("\n",ifelse(tab1t$pourcent>0,"+ ","- "),
525 abs(tab1t$pourcent)," % en ",pasdetemps," ans",sep=""),"")
526
527 }
528 ## table du texte de la tendance / table of the text for the population evolution trend
529 tabTextPent <- data.frame(y=c(max(c(tab1$val,tab1$UL),na.rm=TRUE)*.9),
530 x=median(tab1$annee),
531 txt=ifelse(tendanceSurFigure,c(txtPente1),""),
532 courbe=c(vpan[1]),panel=c(vpan[1]))
533 ## les couleurs / the colors
534 vecColPoint <- c("#ffffff","#eeb40f","#ee0f59")
535 names(vecColPoint) <- c("significatif","infSeuil","0")
536 vecColCourbe <- c("#3c47e0","#5b754d","#55bb1d","#973ce0")
537 names(vecColCourbe) <- c(vpan[1],"carre","presence",vpan[3])
538 vecColHline <- c("#ffffff","#e76060")
539 names(vecColHline) <- c("variation abondance","seuil")
540
541 col <- c(vecColPoint,vecColCourbe,vecColHline)
542 names(col) <- c(names(vecColPoint),names(vecColCourbe),names(vecColHline))
543
544 ## si description graphique en 3 panels
545 if(description) {
546 p <- ggplot(data = dgg, mapping = aes(x = annee, y = val))
547 ## Titre, axes ...
548 p <- p + facet_grid(panel ~ ., scale = "free") +
549 theme(legend.position="none",
550 panel.grid.minor=element_blank(),
551 panel.grid.major.y=element_blank()) +
552 ylab("") + xlab("Annee")+ ggtitle(titre) +
553 scale_colour_manual(values=col, name = "" ,
554 breaks = names(col))+
555 scale_x_continuous(breaks=min(dgg$annee):max(dgg$annee))
556 p <- p + geom_hline(data =hline.data,mapping = aes(yintercept=z, colour = couleur,linetype=type ),
557 alpha=1,size=1.2)
558 if(assessIC){ ############# ONLY FOR THE CONFIDENCE INTERVAL
559 p <- p + geom_ribbon(mapping=aes(ymin=LL,ymax=UL),fill=col[vpan[1]],alpha=.2)
560 p <- p + geom_pointrange(mapping= aes(y=val,ymin=LL,ymax=UL),fill=col[vpan[1]],alpha=.2)
561 }
562 p <- p + geom_line(mapping=aes(colour=courbe),size = 1.5)
563 p <- p + geom_point(mapping=aes(colour=courbe),size = 3)
564 p <- p + geom_point(mapping=aes(colour=catPoint,alpha=ifelse(!is.na(catPoint),1,0)),size = 2)
565 p <- p + geom_text(data=tabTextPent, mapping=aes(x,y,label=txt),parse=FALSE,color=col[vpan[1]],fontface=2, size=4)
566 ggsave(figname, p,width=16,height=21, units="cm")
567 print (figname) ##### CAN BE REMOVED IF YOU DO NOT WANT THE GRAPH TO BE PLOTTED
568 } else {
569
570 p <- ggplot(data = subset(dgg,panel=="Variation abondance"), mapping = aes(x = annee, y = val))
571 ## Titre, axes ...
572 p <- p + facet_grid(panel ~ ., scale = "free") +
573 theme(legend.position="none",
574 panel.grid.minor=element_blank(),
575 panel.grid.major.y=element_blank()) +
576 ylab("") + xlab("Annee")+ ggtitle(titre) +
577 scale_colour_manual(values=col, name = "" ,
578 breaks = names(col))+
579 scale_x_continuous(breaks=min(dgg$annee):max(dgg$annee))
580 p <- p + geom_hline(data =subset(hline.data,panel=="Variation abondance"),mapping = aes(yintercept=z, colour = couleur,linetype=type ),
581 alpha=1,size=1.2)
582
583 if(assessIC){ ############# ONLY FOR THE CONFIDENCE INTERVAL
584 p <- p + geom_ribbon(mapping=aes(ymin=LL,ymax=UL),fill=col[vpan[1]],alpha=.2)
585 p <- p + geom_pointrange(mapping= aes(y=val,ymin=LL,ymax=UL),fill=col[vpan[1]],alpha=.2)
586 }
587 p <- p + geom_line(mapping=aes(colour=courbe),size = 1.5)
588 p <- p + geom_point(mapping=aes(colour=courbe),size = 3)
589 p <- p + geom_point(mapping=aes(colour=catPoint,alpha=ifelse(!is.na(catPoint),1,0)),size = 2)
590 p <- p + geom_text(data=tabTextPent, mapping=aes(x,y,label=txt),parse=FALSE,color=col[vpan[1]],fontface=2, size=4)
591 ggsave(figname, p,width=15,height=9,units="cm")
592 print (figname) ##### CAN BE REMOVED IF YOU DO NOT WANT THE GRAPH TO BE PLOTTED
593 }
594 }
595 ############################################################################################################ fin fonction graphique / end of function for graphical output
596
597
598
599
600 #################################################################################################################### debut de la fonction de moyenne geometrique pondere / start of the geometric weighted mean function
601 geometriqueWeighted <- function(x,w=1) exp(sum(w*log(x))/sum(w))
602 #################################################################################################################### fin de la fonction de moyenne geometrique pondere / end of the geometric weighted mean function
603
604
605
606 ##################################################################################################################### debut de la fonction analyseGroupe / start of the function analyseGroupe
607 ## Analyse par groupe de specialisation a partir des resulats de variation d'abondance par especes / analysis by specialization group based on results of the analysis of population evolution trend
608 #
609
610
611 analyseGroupe <- function(id="france",tabsp=tabsp,donnees=donnees,donneesTrend=donneesTrend,ICfigureGroupeSp=TRUE,powerWeight=2,
612 correctionAbondanceNull = 0.000001,
613 groupeNom = c("generaliste","milieux batis","milieux forestiers","milieux agricoles"),
614 groupeCouleur = c("black","firebrick3","chartreuse4","orange")) {
615
616
617
618
619
620 ## donnees tendances globales / results of the global trends
621 donneesTrend <- subset(donneesTrend, select = c(code_espece,valide,mediane_occurrence))
622
623 ## table de reference espece / reference table for species
624 tabsp <- subset(tabsp, select= c(sp,nom,indicateur, specialisation))
625 donnees <- merge(donnees,donneesTrend,by="code_espece")
626 donnees <- merge(donnees,tabsp,by.x="code_espece",by.y="sp")
627 ## table de correspondance de biais en fonction des medianes des occuerences
628
629
630 nameFileSpe <- paste("Output/",id,"/variationsAnnuellesGroupes_",id, ############# Declare le fichier de sortie des variations annuelles par groupe / declare the name of the outputfile for annual population evolution trend by group
631 ".tabular",sep="" )
632 nameFileSpepng <- paste("Output/",id,"/variationsAnnuellesGroupes_",id, ############# Declare le fichier de sortie graphique des variations annuelles par groupe / declare the name of the graphical output file for annual population evolution trend by group
633 ".png",sep="" )
634
635 grpe <- donnees$specialisation
636
637 ####### valeur seuil sont obtenues à partir de simulations / threshold values are obtained from simulations
638 ff <- function(x,y) max(which(y<=x)) ## fonction pour recherche le poid associé à valeur max parmi valeur seuil d'occurence inferieur ou egale à occurence mediane obs / function to retrieve the weight associated with the max occurence threshold equal or smaller than the occurence mediane observed
639
640 IncertW <- ifelse(donnees$valide=="Incertain",tBiais$biais[sapply(as.vector(donnees$mediane_occurrence),ff,y=tBiais$occurrenceMed)],1) ## pr verifier poids de l'espèce dans analyse, récupére seuil occurence minimum pour lequel tendance pas bonne, et compare avec mediane occurence des données / to check the weight of species in the analysis, this retrieve occurence threshold with wich real occurence measured on data are compared in order to verify the accuracy of the trend measurment
641 ## poids du a la qualite de l'estimation
642 # erreur_stW <- 1/((donnees$erreur_st+1)^powerWeight)
643 # erreur_stW <- ifelse( is.na(donnees$IC_superieur),0,erreur_stW)
644 erreur_stW <- ifelse(is.na(donnees$IC_superieur),0,1)##### si pas d'interval de confiance met 0 et donne un poid de 0 à l'esps (voir ci dessous) / if no confidence interval calculated give a weight of 0 for the sps
645 ## calcul du poids total de chaque espèce / calcul of the weight of each species
646 W <- IncertW * erreur_stW
647
648 ## variable de regroupement pour les calculs par groupe de specialisation et par an / variables gathered to identify group for the calculation (per specialization and per year)
649 grAn <- paste(donnees$specialisation,donnees$annee,sep="_")
650 ## data frame pour le calcul / dataframe made for the calcul
651 dd <- data.frame(grAn,annee = donnees$annee, grpe,W,ab=donnees$abondance_relative,ICinf= donnees$IC_inferieur, ICsup= ifelse(is.na(donnees$IC_superieur),10000,donnees$IC_superieur))
652 ## table resumer de tous les poids / table to sum up the weights of each species depending on the incertainty in the calcul of the poulation evolution trends
653 ddd <- data.frame(code_espece = donnees$code_espece,nom_espece = donnees$nom_espece,annee = donnees$annee,
654 groupe_indicateur = grpe,
655 poids_erreur_standard = round(erreur_stW,3), poids_incertitude = round(IncertW,3),poids_final = round(W,3),
656 abondance_relative=donnees$abondance_relative,
657 IC_inferieur= donnees$IC_inferieur,
658 IC_superieur= ifelse(is.na(donnees$IC_superieur),10000,donnees$IC_superieur),
659 valide = donnees$valide, mediane_occurrence = donnees$mediane_occurrence)
660
661 nomFileResum <- paste("Output/",id,"/donneesGroupes_",id, ###### declaration du nom du repertoire et des fichiers de sortie / declaring the name of the output folder and files
662 ".tabular",sep="" )
663 write.table(ddd,nomFileResum,row.names=FALSE,sep="\t",dec=".",fileEncoding="UTF-8")
664 cat("-->",nomFileResum,"\n")
665
666 ## calcul des moyennes ponderees par groupe par an et pour les abondance et les IC / calcul of weighted means per specialization group and per year for the abundance and confidence interval
667 for(j in 5:7) dd[,j] <- ifelse(dd[,j]==0,correctionAbondanceNull,dd[,j])
668 ag <- apply(dd[,5:7], 2, ######## sur les abondances relatives, les ICinf et ICsup
669 function(x) {
670 sapply(split(data.frame(dd[,1:4], x), dd$grAn), ###### fait les moyennes pondérés par groupe grAn / calculate the weighted mean by group grAn
671 function(y) round(geometriqueWeighted(y[,5], w = y$W),3))
672 })
673 ## gg <- subset(dd,as.character(dd$grAn)=="milieux forestier_2014") #############################################################
674
675 ag <- ifelse(is.na(ag),1,ag)
676 ag <- as.data.frame(ag)
677 ag$grAn <- rownames(ag)
678 dbon <- subset(donnees,valide=="bon")
679 dIncert <- subset(donnees,valide=="Incertain")
680 ## calcul nombre d'espece "bonne" pour le calcul / calculating the number of species with low level of incertainty, "good" species
681 bon <- tapply(dbon$nom,dbon$specialisation,FUN=function(X)length(unique(X)) )
682 bon <- ifelse(is.na(bon),0,bon)
683 tbon <- data.frame(groupe=names(bon),bon)
684 ## calcul nombre d'especes "incertaines" pour le calcul / calculating the number of species with high level of incertainty, "bad" species
685 Incert <- tapply(dIncert$nom,dIncert$specialisation,FUN=function(X)length(unique(X)) )
686 Incert <- ifelse(is.na(Incert),0,Incert)
687 tIncert <- data.frame(groupe=names(Incert),Incertain=Incert)
688
689 tIncert <- merge(tIncert,tbon,by="groupe")
690
691 ## table de données avec les moyennes ponderees par groupe / table of the data with the weighted mean by group
692 da <- merge(unique(dd[,1:3]),ag,by="grAn")[,-1]
693 colnames(da) <- c("annee","groupe","abondance_relative","IC_inferieur","IC_superieur")
694
695 da$annee <- as.numeric(da$annee)
696 da <- merge(da,tIncert,by="groupe") #### ajoute le nombre d'espece "incertaines" et "bonne" aux resultats / add the number of "good" and "bad" species to the overall resutls
697 da <- subset(da, groupe != "non")
698 colnames(da)[6:7] <- c("nombre_especes_incertaines","nombre_espece_bonnes")
699 a <- data.frame(id,da)
700 write.table(da,file=nameFileSpe,row.names=FALSE,quote=FALSE,sep="\t",dec=".",fileEncoding="UTF-8")
701
702 cat("-->",nameFileSpe,"\n")
703 yearsrange <- c(min(da$annee),max(da$annee))
704
705 ## figure par ggplot2 / plots with ggplot2
706 titre <- paste("Variation de l'indicateur groupe de specialisation",sep="")
707
708 vecCouleur <- setNames(groupeCouleur,groupeNom)
709 #browser()
710 p <- ggplot(data = da, mapping = aes(x = annee, y = abondance_relative, colour=groupe,fill=groupe))
711 p <- p + geom_hline(aes(yintercept = 1), colour="white", alpha=1,size=1.2)
712 if(ICfigureGroupeSp)
713 p <- p + geom_ribbon(mapping=aes(ymin=IC_inferieur,ymax=IC_superieur),linetype=2,alpha=.1,size=0.1)
714 p <- p + geom_line(size=1.5)
715 p <- p + ylab("") + xlab("Annee")+ ggtitle(titre)
716 if(!is.null(groupeNom)) p <- p + scale_colour_manual(values=vecCouleur, name = "" )+
717 scale_x_continuous(breaks=unique(da$annee))
718 if(!is.null(groupeNom)) p <- p + scale_fill_manual(values=vecCouleur, name="")
719 p <- p + theme(panel.grid.minor=element_blank(), panel.grid.major.y=element_blank())
720 ggsave(nameFileSpepng, p,width=17,height=10,units="cm")
721
722 # cat(" <==",nameFileSpepng,"\n")
723
724 ## calul pour chaque groupe une pente de regression d'evolution des abondances sur la periode étudiée / calculating for each group the regression slope for the abundance evolution on the studied period
725 vecSpe <- unique(da$groupe)
726 datasum <- data.frame(groupe=NULL,tendance=NULL,pourcentage_variation=NULL)
727 for(spe in 1:4){
728 # print(spe)
729 subtab <- subset(da,groupe==vecSpe[spe])
730 if(nrow(subtab)>1) {
731 sumlm <- summary(lm(abondance_relative~annee,data=subtab)) ##### recupère les resultats du modèle linéaire / retrieve the results of the linear model
732 subdatasum <- data.frame(groupe=vecSpe[spe],
733 tendance=round(sumlm$coefficients[2,1],3),
734 pourcentage_variation=round(sumlm$coefficients[2,1]*(nrow(subtab)-1)*100,3)) #### assemble les resultats pour en faire une sortie / bring together the results for an output file
735 datasum <- rbind(datasum,subdatasum)
736
737 }
738
739 }
740 datasum <- merge(datasum,tIncert,by="groupe") ####
741 datasum <- data.frame(id,datasum)
742 #datasum$cat_tendance_EBCC <- affectCatEBCC(trend,pVal,ICinf,ICsup
743 namefilesum <- paste("Output/",id,"/tendancesGlobalesGroupes_",id,
744 ".tabular",sep="" )
745 write.table(datasum,file=namefilesum,row.names=FALSE,quote=FALSE,sep="\t",dec=".",fileEncoding="UTF-8")
746 cat("-->",namefilesum,"\n")
747 }
748
749 ################################################################################################################## fin de la fonction analyseGroupe / end of the function analyseGroupe
750
751
752
753
754
755
756
757 ################################################################################################################### debut de la fonction check_file / start of the function check_file
758 # Fonction pour verifier les données d'entrée / General function to check integrity of input file. Will check numbers and contents of variables(colnames).
759 #return an error message and exit if mismatch detected
760 #Faut rentrer le nom du jeu de données, le nbre et le nom des variables / Enter dataset name, expected number and names of variables. + an exit error message to guide user.
761
762 check_file<-function(dataset,err_msg,vars,nb_vars){
763 if(ncol(dataset)!=nb_vars){ #Verifiction de la présence du bon nb de colonnes, si c'est pas le cas= message d'erreur / checking for right number of columns in the file if not = error message
764 cat("\nerr nb var\n")
765 stop(err_msg, call.=FALSE)
766 }
767
768 for(i in vars){
769 if(!(i %in% names(dataset))){
770 stop(err_msg,call.=FALSE)
771 }
772 }
773 }
774
775 #####################################################################################################################
776