Mercurial > repos > ecology > stoc_trend_indic

comparison FunctExeTemporalAnalysisofIndicatorTrait.r @ 0:916b49d725ba 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:35:29 -0400
parents
children 5c244408661f
comparison
equal deleted inserted replaced
-1:000000000000 0:916b49d725ba
1 #!/usr/bin/env Rscript
2
3 ########################################################################################################################################################################
4 ############## FUNCTION TO CALCULATE AND PLOT VARIATION IN TIME OF INDICATOR OR MEAN TRAIT VALUE OF COMMUNITIES function:csi_cti_ctri ##############################
5 ########################################################################################################################################################################
6
7 #### Based on Romain Lorrillière R script
8 #### Modified by Alan Amosse and Benjamin Yguel for integrating within Galaxy-E
9 ### made with R version 3.5.1
10
11
12 suppressMessages(library(RODBC)) ##Version: 1.3-15
13 suppressMessages(library(reshape)) ##Version: 0.8.8
14 suppressMessages(library(data.table)) ##Version: 1.12.0
15 suppressMessages(library(rgdal)) ##Version: 1.3-4
16 suppressMessages(library(lubridate)) ##Version: 1.7.4
17 #suppressMessages(library(RPostgreSQL)) ##Version: 0.6-2
18 suppressMessages(library(doBy)) ##Version: 4.6-2
19 suppressMessages(library(arm)) ##Version: 1.10-1
20 suppressMessages(library(ggplot2)) ##Version: 3.1.0
21 suppressMessages(library(scales)) ##Version: 1.0.0
22 suppressMessages(library(mgcv)) ##Version: 1.8-24
23 #suppressMessages(library(visreg)) ##Version: 2.5-0
24 suppressMessages(library(plyr)) ##Version: 1.8.4
25 #suppressMessages(library(lme4)) ##Version: 1.1-18-1
26 #suppressMessages(library(lmerTest)) ##Version: 3.1-0
27 suppressMessages(library(glmmTMB)) ###Version: 0.2.3
28
29
30 ###########
31 #delcaration des arguments et variables/ declaring some variables and load arguments
32
33 args = commandArgs(trailingOnly=TRUE)
34 #for( i in 1:length(args)){print(args[i])}
35 if (length(args)<12) {
36 stop("At least 12 arguments must be supplied :\n- An input dataset filtered (.tabular). May come from the filter rare species tool.\n- A species detail table (.tabular)\n- A species ssi/sti table.\n- A table with plots coordinates.\n- table with csi calculated before 2001.\n\n", call.=FALSE) #si pas d'arguments -> affiche erreur et quitte / if no args -> error and exit1
37 } else {
38 Datafiltered<-args[1] ###### Nom du fichier avec extension ".typedefichier", peut provenir de la fonction "FiltreEspeceRare" / file name without the file type ".filetype", may result from the function "FiltreEspeceRare"
39 tabSpecies<-args[2] ###### Nom du fichier avec extension ".typedefichier", fichier mis à disposition dans Galaxy-E avec specialisation à l'habitat des especes et si espece considérée comme indicatrice / file name without the file type ".filetype", file available in Galaxy-E containing habitat specialization for each species and whether or not they are considered as indicator
40 tabtrait<-args[3] ##### Nom du fichier avec extension ".typedefichier", fichier mis à disposition dans Galaxy-E avec degre de specialisation de l espece et affinite thermique /file name without the file type ".filetype", file available in Galaxy-E containing specilalization degree as well as thermic preferences
41 coordCarre<-args[4] #### Nom du fichier avec extension ".typedefichier", fichier mis à disposition dans Galaxy-E avec les coordonnees gps des carres /file name without the file type ".filetype", file available in Galaxy-E containing gps coordinates of the plots
42 Var <- args[5] #### Nom du trait dans fichier de traits "nomdutrait" exemple: "ssi" pour l'indice de specialisation par sps / Name of the trait in the file containing trait data
43 indicator <- args[6] #### Nom de l'indicateur ou du trait par communauté ex pour ssi c'est csi calculé au niveau communauté / Name of the indicator or the trait per community ex: for the ssi, it is the csi measured at the community level
44 methode <- args[7] #### Methode d'analyse de l'evolution du trait ou de l'indicateur, lmer pour modèloe mixte seul ou gam pour generalized additive model / name of the models used to analyze evolution of mean trait or indicator
45 dd <- args[8] ##### Nom du fichier si déjà un fichier avec trait moyen par communauté, avec une colonne annee appelé "year" et une colonne plot appelé "carre" correspondant à l'echelle des communautés etudiées / name of the file if a file with the mean trait value per community is already prepared with one column named "year" for the year, one column named "carre" for the plots (the scale of the community measurment)
46 id<-args[9]#Id name for output res repo
47 plot_smooth<-args[10]#TRUE or FALSE
48 ic<-args[11]#TRUE or FALSE
49 source(args[12])
50 }
51
52 #Import des données / Import data
53 tabCLEAN <- read.table(Datafiltered,sep="\t",dec=".",header=TRUE) #### charge le fichier de données d abondance / load abundance of species
54 tabsp <- read.table(tabSpecies,sep="\t",dec=".",header=TRUE) #### charge le fichier de donnees sur nom latin, vernaculaire et abbreviation, espece indicatrice ou non / load the file with information on species specialization and if species are indicators
55
56 vars_tabCLEAN<-c("carre","annee","espece","abond")
57 err_msg_tabCLEAN<-"The input dataset filtered doesn't have the right format. It need to have the following 4 variables :\n- carre\n- annee\n- espece\n- abond\n"
58
59 vars_tabsp<-c("espece","nom","nomscientific","indicateur","specialisation")
60 err_msg_tabsp<-"\nThe species dataset filtered doesn't have the right format. It need to have the following 4 variables :\n- espece\n- nom\n- nomscientific\n- indicateur\n- specialisation\n"
61
62 check_file(tabCLEAN,err_msg_tabCLEAN,vars_tabCLEAN,4)
63 check_file(tabsp,err_msg_tabsp,vars_tabsp,5)
64
65
66 if(!dd==""){
67 vars_dd<-c("carre","year","longitude_grid_wgs84","latitude_grid_wgs84","indic") #### si vous avez déjà votre tableau d'analyse indic correspond au trait moyen par communauté ou au calcul de l'indicateur
68 err_msg_dd<-"\nThe dataset for analysis doesn't have the right format. It need to have the following 5 variables :\n- carre\n- year\n- longitude_grid_wgs84\n- latitude_grid_wgs84\n- indic\n"
69 check_file(dd,err_msg_dd,vars_dd,5)
70 dd <- read.table(dd,sep="\t",dec=".",header=TRUE) #### charge le fichier pour analyse si déjà construit (voir ci dessus pour les détails ) / load the required file for the analysis if already prepared (see above for details)
71 }else{
72 dd<-NULL
73 }
74
75 spTrait=read.table(tabtrait,sep="\t",dec=".",header=TRUE) ############# species_indicateur_fonctionnel.csv pour le STOC sinon fichier avec traits pour calcul du trait moyen par communauté / file with the trait for the community weighted mean calculation
76 coordCarre=read.table(coordCarre,sep="\t",dec=".",header=TRUE) ######## carre.csv charge les coordonnées des carrés qui sont utilisés comme covariable / load the gps coordinates of the plots, is used as covariable in the models
77
78 dir.create(paste("Output/",sep=""),recursive=TRUE,showWarnings=FALSE)##### Creation du dossier de sortie
79 #cat(paste("Create Output/","\n",sep=""))
80
81 ############################# The function
82
83
84 csi_cti_ctri <- function(tabCLEAN=tabCLEAN,coordCarre=coordCarre,spTrait=spTrait,dd=NULL,ic=TRUE, Var="ssi",indicator="csi", methode="gam", ####### Var= nom du trait dans le fichier de trait (pour le calcul du csi, le trait est ssi par exemple); Indicator= nom de l'indicateur ou du trait moyen par communauté ; methode: choisir modele "gam" ou "glmmtmb" ; ic pour calcul de l'interval de confiance plus rapide sans mais moins fiable / Var=name of the trait in the trait file; Indicator= name of the indicator or of the community weighted mean trait used in graphical output and output files; methode is the statistical model use for the analysis "gam" or "glmmtmb" ; ic is for the calculation of confidence interval faster without but less reliable
85 firstYear=NULL,lastYear=NULL,altitude=800,departement=NULL,onf=TRUE,distance_contact=NULL, #### altitude, departement onf, distance de contact = Argument non utilise, se trouvait dans requete sql / altitude, departement onf, distance de contact = not use anymore was in a postgres request
86 spExcluPassage1=c("MOTFLA","SAXRUB","ANTPRA","OENOEN","PHYTRO"),# (Prince et al. 2013 Env. Sc. and Pol.) + "OENOEN","PHYTRO" avis d'expert F. Jiguet, #### Argument non utilise, se trouvait dans requete sql / not use anymore was in a postgres request
87 seuilAbondance=.99,plot_smooth=TRUE, ###### init_1989 si TRUE, option que pour csi et besoin du fichier des csi calculés sur les données avant 2001 (pas forcement fiable car protocole un peu different) / init_1989 if TRUE, only working for csi, and use calculation of csi based on data before 2001 (protocol was bit different, not totally reliable)
88 champSp = "code_sp", sp=NULL,champsHabitat=FALSE, #### Argument non utilise, se trouvait dans requete sql / not use anymore was in a postgres request
89 anglais=FALSE,seuilSignif=0.05,##### #### anglais=FALSE Argument non utilise, se trouvait dans requete sql / not use anymore was in a postgres request
90 couleur="#4444c3",
91 titreY=indicator,titreX="Années",titre=indicator,
92 savePostgres=FALSE,output=FALSE, ##### OPTION "output" pour afficher le resultat dans R / OPTION "output" is only to show the result in the R window
93 operateur=c("Lorrilliere Romain","lorrilliere@mnhn.fr"), encodingSave="ISO-8859-1",fileName="dataCSI",id="France"){ ####### nom des fichiers de sorties et de l'operateur / name of the output files and of the operator
94
95
96 start <- Sys.time()
97 dateExport <- format(start,"%Y-%m-%d")
98
99 if(is.null(firstYear)) firstYear <- 2001
100 if(is.null(lastYear)) lastYear <- 9999
101 if(is.null(altitude)) altitude <- 10000
102
103 ############################################# calcul de l'indicateur ou du trait moyen pondéré par communauté / Calculation of the weighted mean trait or of the indicator per community i.e. calculation of the csi, cti, ctri (corresponds to Community weighted mean or CWM for traits)
104 if (is.null(dd)){
105 colnames(spTrait)[colnames(spTrait) == Var] <- "trait"
106 spTrait$trait <- spTrait$trait
107 ###browser()
108 tabCLEAN$trait <- spTrait$trait[match(tabCLEAN$espece,spTrait$pk_species)] ### recupere donnee du trait par espece calcule / retrieve trait data for each species
109 tabCLEAN=na.omit(tabCLEAN) ##### pour faire les moyennes pondérées sur les espèces avec des données de trait (donc pas de prise en compte des sps sans traits dans l'abondance totale par carré)
110 traitcarre <- aggregate(trait*abond~annee+carre,tabCLEAN,sum) ### somme des traits par annee et carre pondere par les abondances / sum of the trait per year and per plots weighted by abundances
111 abcarre <- aggregate(abond~annee+carre,tabCLEAN,sum) ### somme des abondances totales par annee et carre / sum of total abundance per year and plots
112 indic <- traitcarre[,3]/abcarre[,3] #### le trait moyen par carre = indicateur par carre et annee / mean trait per plots = indicator per year and plots
113 dd <- data.frame(indic,traitcarre$carre,traitcarre$annee)
114 names(dd)[2] <- "carre"
115 names(dd)[3] <- "year"
116 dd$longitude_grid_wgs84 <- coordCarre$longitude_grid_wgs84[match(dd$carre,coordCarre$pk_carre)] #### recupere coordonnées gps / retrieve gps coordinates
117 dd$latitude_grid_wgs84 <- coordCarre$latitude_grid_wgs84[match(dd$carre,coordCarre$pk_carre)] #### recupere coordonnées gps / retrieve gps coordinates
118 dd$id_plot <- dd$carre ### id_plot nom données aux carrés dans le script /id_plot is use as the name of the plot in the following script
119 }else{
120 colnames(dd)[colnames == "indicator"] <- "indic"
121 }
122
123 ############################################ fin du calcul de l'indicateur ou du trait moyen pondéré par communauté / end of the calculation of the indicator or the mean trait value per community
124 annee <- sort(unique(dd$year))
125 nban <- length(annee)
126 pasdetemps <-nban-1
127
128 if(methode == "gam") {
129 cat("Methode: gam\n")
130 ## Utilisation des modèles GAMM pour obtenir les tendances d evolution par an du csi cti ou ctri !!!! Marche pas si peu de données !!!! / Use of GAMM model for the estimation of the annual variations of the csi cti or ctri !!! does not work with few data !!!
131 cat("\nEstimation de la variation annuelle ",indicator,"~ factor(year)+s(longitude_grid_wgs84,latitude_grid_wgs84,bs='sos'),random=reStruct(object = ~ 1| id_plot,correlation=corAR1(form=~year)\n",sep="")
132 gammf <- gamm(indic ~ factor(year)+s(longitude_grid_wgs84,latitude_grid_wgs84,bs="sos"), data=dd,random=reStruct(object = ~ 1| id_plot, pdClass="pdDiag"),correlation=corAR1(form=~year)) #### spline sur les coordonnées, effet aleatoire sur les carres, methode autoregressive sur l'année N-1 / spline on the gps coordinates, random effect on the plots, autoregressive method on the year-1
133 # gammf <- gamm(indic ~ factor(year)+s(longitude_grid_wgs84,latitude_grid_wgs84,bs="sos"), data=dd,random=reStruct(object = ~ 1| id_plot, pdClass="pdDiag"),correlation=corAR1(form=~year)) #### spline sur les coordonnées, effet aleatoire sur les carres, methode autoregressive sur l'année N-1 / spline on the gps coordinates, random effect on the plots, autoregressive method on the year-1
134 sgammf<-summary(gammf$gam)
135 coefdata=coefficients(gammf$gam) ### recupere les coefficient de regression de la variable "annee" / retrieve the regression coefficient of the variable "year"
136 coefannee <- c(0,sgammf$p.coeff[2:nban]) ### meme chose que au dessus / same as before
137 erreuran <- c(0,sgammf$se[2:nban])### recupere les erreurs standard des coefficient de regression de la variable "annee" / retrieve the standard errors of the regression coefficient of the variable "year"
138 pval <- c(1,sgammf$p.pv[2:nban])### recupere les p value de la variable "annee" / retrieve the p value of the variable year
139
140
141 ## calcul des intervalles de confiance / confidence interval calculation
142 if(ic) {
143 # gammf.sim <- sim(gammf) ###################### VERSION ROMAIN mais fct sim() ne marche pas avec GAMM / old version using function sim() but did not work with Gamm models
144 # ic_inf_sim <- c(0,tail(apply(coef(gammf.sim), 2, quantile,.025),pasdetemps))
145 # ic_sup_sim <- c(0,tail(apply(coef(gammf.sim), 2, quantile,.975),pasdetemps))
146 icalpha05 <- as.data.frame(confint(gammf$gam))[2:nban,1:2] ########## VERSION BENJ
147 ic_inf_sim <- icalpha05[,1]
148 ic_inf_sim <- c("NA",ic_inf_sim[1:nban-1])
149 ic_sup_sim <- icalpha05[,2]
150 ic_sup_sim <- c("NA",ic_sup_sim[1:nban-1])
151 } else{
152 ic_inf_sim <- "not assessed"
153 ic_sup_sim <- "not assessed"
154 }
155
156 tabfgamm <- data.frame(model = "gamm factor(year) plot",annee,coef=coefannee,se = erreuran,pval,signif=pval<seuilSignif,Lower_ci=ic_inf_sim,upper_ci=ic_sup_sim,indicator=indicator) #### recupère les resultats des modèles avec interval de confiance / retrieve results of the models used with confidence interval
157 write.table(tabfgamm,paste("Output/",indicator,"_gammParannee_",id,".tabular",sep=""),row.names=FALSE,sep="\t")
158 gg <- ggplot(data=tabfgamm,aes(x=annee,y=coef))
159 gg= gg + geom_line(size=1.5,colour=couleur)+ geom_point(size=3,colour=couleur) + geom_point(size=1.5,colour="white")
160
161 if (ic){
162 gg <- gg + geom_errorbar(aes(ymin=coef-se, ymax=coef+se), width=0,colour=couleur,alpha=0.5)
163 gg <- gg + geom_ribbon(aes(ymin=coef-se, ymax=coef+se),fill = couleur,alpha=.2)
164 }
165 gg <- gg + labs(y=indicator,x="annee")+scale_x_continuous(breaks=pretty_breaks())
166 ggsave(paste("Output/fig",indicator,"_carre_",id,".png",sep=""),gg)
167
168
169 if(plot_smooth) { #### Representation graphique de l'evolution annuelle des indicateurs / Graphical representation of the annual evolution of the indicators
170 cat("\nGam pour la figure ",indicator,"~s(year),random=reStruct(object = ~ 1| id_plot,correlation=corAR1(form=~year)\n",sep="")
171 ## create a sequence of temperature that spans your temperature #####not use anymore
172 ## http://zevross.com/blog/2014/09/15/recreate-the-gam-partial-regression-smooth-plots-from-r-package-mgcv-with-a-little-style/ #### method for the plot
173 ####dd$yearf=factor(dd$year) PAS BON j'ai modifié la ligne suivante en mettant s(year) et plus s(yearsf)
174 gammgg <- gamm(indic ~ s(year), data=dd,random=reStruct(object = ~ 1| id_plot, pdClass="pdDiag"),correlation=corAR1(form=~year)) #### spline sur l'année, effet aleatoire des carres sur ordonnée à l'origine, methode autoregressive sur l'année N-1 / spline on the year, random effect of the plots on the intercept, autoregressive method on the year-1
175 maxyear<-max(dd$year)
176 minyear<-min(dd$year)
177 year.seq<-sort(unique(c(minyear:maxyear,(seq(minyear, maxyear,length=1000)))))
178 year.seq<-data.frame(year=year.seq)
179 # predict only the temperature term (the sum of the ########### ???? not use anymore
180 # term predictions and the intercept gives you the overall########### ???? not use anymore
181 # prediction)########### ???? not use anymore
182 preds<-predict(gammgg$gam, newdata=year.seq, type="terms", se.fit=TRUE) #### Utilise model pour predire les valeurs de indic sur sequence d'années defini au dessus \ Use of the model to predict value of the indicator in the year sequence define above
183 # set up the temperature, the fit and the upper and lower########### ???? not use anymore
184 # confidence interval########### ???? not use anymore
185 year <-year.seq$year
186 realYear <- sort(unique(dd$year))
187 fit<-as.vector(preds$fit)
188 init <- fit[1]
189 fit.up95 <- fit-1.96*as.vector(preds$se.fit)
190 fit.low95 <- fit+1.96*as.vector(preds$se.fit)
191 # ggGamData <- data.frame(year=year, csi=fit,ic_low95 = fit.low95, ic_up95 = fit.up95)
192 fit <- fit - init ### Réechelonne les predictions du modèle sur la 1ère valeure de la prediction ? ne sait pas pourquoi
193 fit.up95 <- fit.up95 - init ### Réechelonne IC superieur sur la 1ère valeure de la prediction ? ne sait pas pourquoi
194 fit.low95 <- fit.low95 - init ### Réechelonne IC inferieur sur la 1ère valeure de la prediction ? ne sait pas pourquoi
195 ggGamData <- data.frame(year=year, indic=fit,ic_low95 = fit.low95, ic_up95 = fit.up95) ####### Recupère les resultats des modèles / retrieve the results of the models
196
197
198 ## The ggplot:
199 gg <- ggplot(data=ggGamData,aes(x=year,y=indic))+ geom_line(size=1,colour=couleur)
200 if (ic) {
201 gg <- gg + geom_ribbon(aes(ymin=ic_low95, ymax=ic_up95),fill = couleur,alpha=.2)
202 }
203 gg <- gg + geom_point(data = subset(ggGamData,year %in% realYear),size=3,colour=couleur) + geom_point(data = subset(ggGamData,year %in% realYear),size=1.5,colour="white")
204 gg <- gg + labs(y=titreY,x=titreX,title=titre)+scale_x_continuous(breaks=pretty_breaks())
205 ggsave(paste("Output/fig",indicator,"_plot",id,".png",sep=""),gg)
206 #cat("\n--> Output/fig",indicator,"_plot",id,".png\n",sep="")
207
208 tabPredict <- subset(ggGamData,year %in% realYear)########### Tableau des resultats pour ne prendre que les valeurs d'IC pour l'année pas entre les années (spline sur annee) !!!plus utilisé!! / Table of the results not taking confidence interval between year but at each year (because of the spline of year)
209 colnames(tabPredict)[1:2] <- c("annee",paste(indicator,"_predict",sep=""))
210 tabgamm <- merge(tabfgamm,tabPredict,by="annee") #### Desactivation car merge sortie de modèles différents (le modèle dont on tire les coef de regression pour année, avec spli ne sur les coordonnées geo vs celui pour faire la figure avec splin sur année uniquement) / not use anymore (as before) because use the results of the restricted model with the spline on the year while the better analysis is on full model with the spline on gps coordinates
211 ##tabgamm <- tabfgamm #### remplace la ligne au dessus / replace the line above
212 write.table(tabPredict,paste("Output/",indicator,"_gammsmooth",id,".tabular",sep=""),row.names=FALSE,sep="\t")
213 }else {
214 tabgamm = tabfgamm
215 }
216 ###### CROCHET deplacer apres le dernier ggplot ici mais dans version Romain CROCHET placé AVANT le dernier ggplot qui n'était pas dans la condition if(init_1989)
217 cat("\nEstimation de la tendence ",indicator,"~ year+s(longitude_grid_wgs84,latitude_grid_wgs84,bs='sos'),random=reStruct(object = ~ 1| id_plot,correlation=corAR1(form=~year)\n")
218 gammc <- gamm(indic~year+s(longitude_grid_wgs84,latitude_grid_wgs84,bs="sos"), data=dd,random=reStruct(object = ~ 1| id_plot, pdClass="pdDiag"),correlation=corAR1(form=~year))### spline sur les coordonnées, effet aleatoire des carres sur ordonnée à l'origine, methode autoregressive sur l'année N-1 / spline on the gps coordinates, random effect of the plots on intercept, autoregressive method on the year-1
219 sgammc=summary(gammc$gam)
220 coefannee <- sgammc$p.coeff[2] #### coefficient de regression de la variable année / regression coefficient of the variable "year"
221 ## erreur standard / standard error
222 erreuran <- sgammc$se[2]
223 ## p value
224 pval <- sgammc$p.pv[2]
225
226 #### Calcul des intervalles de confiances / calculation of the confidence intervals
227 if(ic) {
228 # gammc.sim <- sim(gammc)###################### VERSION ROMAIN mais fct sim() marche pas avec Gamm / old version using function sim() but did not work with Gamm models
229 # ic_inf_sim <- c(0,tail(apply(coef(gammc.sim), 2, quantile,.025),pasdetemps))
230 # ic_sup_sim <- c(0,tail(apply(coef(gammc.sim), 2, quantile,.975),pasdetemps))
231 icalpha052 <- as.data.frame(confint(gammc$gam))[2,1:2] ########## VERSION BENJ
232 ic_inf_sim2 <- icalpha052[,1]
233 ic_sup_sim2 <- icalpha052[,2]
234 } else{
235 ic_inf_sim2 <- "not assessed"
236 ic_sup_sim2 <- "not assessed"
237 }
238 tabcgamm <- data.frame(model = "gamm numeric(year) plot",annee = NA,coef = coefannee,se = erreuran,pval,signif = pval<seuilSignif, indicator= indicator , Lower_ci = as.factor(ic_inf_sim2), upper_ci = as.factor(ic_sup_sim2), csi_predict =NA ,ic_low95 =NA,ic_up95=NA)#### recupère les resultats des modèles avec interval de confiance / retrieve results of the models used with confidence interval
239 ########### MODIF tabcgamm en remplacant ic_low95 ic_up95 par Lower_ci et upper_ci pour coller avec les sorties des modèles "pour les stats" et non celui utilisé pour le graphe uniquement, et rajout des colonnes spécifique a model ggsmooth pour les garder
240 #### ai rajouté aussi as.factor(ic) car ne savait pas pourquoi mais tabfgamm sont en facteur et besoin de la meme class pour rbind() ci dessous
241 tabcgamm <- tabcgamm[,colnames(tabgamm)] ##
242 tabgamm <- tabgamm[,colnames(tabcgamm)] #### recupère que les colonnes de tabcgamm donc perds les infos du modèle du ggsmooth si pas declarer dans le tableau tabcgamm (maintenant c fait)
243 tabgamm <- rbind(tabgamm,tabcgamm)
244 write.table(tabgamm,paste("Output/",indicator,"_gammCOMPLET_",id,".tabular",sep=""),row.names=FALSE,sep="\t")
245 #cat("\n --> Output/",indicator,"_gammPlot_",id,".tabular\n",sep="")
246 }
247
248
249 if (methode == "glmmtmb") {
250 #cat("Method : lmer \n")
251 cat("Method : glmmTMB \n")
252 ###################
253 ### Utilisation des modèles mixtes pour obtenir les tendances d evolution par an du csi cti ou ctri / Use of mixte model for the estimation of the annual variations of the csi cti or ctri
254 #cat("\nEstimation de la variation annuelle lmer(",indicator,"~ factor(year)+(1|id_plot)\n",sep="")
255 cat("\nEstimation de la variation annuelle glmmTMB(",indicator,"~ factor(year)+(1|id_plot)\n",sep="")
256 #md.f <- lmer(indic~ factor(year)+(1|id_plot),data=dd) ##### effet aleatoire liés aux carrés sur l'ordonnée à l'origine / random effects of plots on intercept
257 md.f <- glmmTMB(indic~ factor(year)+(1|id_plot),data=dd)
258 smd.f <- summary(md.f)
259 # coefdata.f <- as.data.frame(smd.f$coefficients) ### version pour sortie lmer()
260 coefdata.f <- as.data.frame(smd.f$coefficients$cond[-1,])
261 coefdata.f <- data.frame(model="Annual fluctuation", variable = rownames(coefdata.f),coefdata.f)
262 # ggdata <<- data.frame(year=c(as.numeric(substr(coefdata.f$variable[-1],13,16))),##### version pour sortie lmer()
263 ggdata <- data.frame(year=c(as.numeric(substr(coefdata.f$variable,13,17))),
264 estimate=c(coefdata.f$Estimate),
265 se=c(coefdata.f$Std..Error)) ##################### resultat du modèle / results of the models
266 #ggdata$estimate <- ggdata$estimate
267 #ggdata$se.supR <- ggdata$estimate + ggdata$se ############################################################################## METHODE ROMAIN
268 #ggdata$se.infR <- ggdata$estimate - ggdata$se
269 # ggdata$estimate2 <- c(coefdata.f$Estimate[1],coefdata.f$Estimate[1] + coefdata.f$Estimate[-1])
270 # ggdata$se.sup2 <- ggdata$estimate2 + ggdata$se
271 # ggdata$se.inf2 <- ggdata$estimate2 - ggdata$se
272 #prof <- profile(md.f) #### Nouvel interval de confiance avec utilisation du logarithme des ecarts types / logarithms of standard deviations are used, while varianceProf converts from the standard-deviation to the variance scale
273 MODconfint <- confint (md.f) #### plus rapide de ne pas passer par la fonction profile et pas indispensable fonctionne aussi directement sur modele mixte md.f / more rapid using both function profile and confint but works also directly on output of the model
274 se.sup <- MODconfint[2:nban,2]#### [2:nban+2,2] version pour sortie lmer()
275 se.inf <- MODconfint[2:nban,1]#### [2:nban+2,2] version pour sortie lmer()
276 if (ic) {
277 ggdata$se.sup <- se.sup
278 ggdata$se.inf <- se.inf
279 } else{
280 ggdata$se.sup <- "not assessed"
281 ggdata$se.inf <- "not assessed"
282 }
283 coefdata.f$se.inf <- ggdata$se.inf
284 coefdata.f$se.sup <- ggdata$se.sup
285 #gg <<- ggplot(ggdata,aes(x=year,y=estimate))+ geom_ribbon(ymin=ggdata$se.infR,ymax=ggdata$se.supR,alpha=.25)+geom_errorbar(ymin=ggdata$se.infR,ymax=ggdata$se.supR,width=0,alpha=.25)+ geom_point() + geom_line() + ylim(min(ggdata$se.infR),max(ggdata$se.supR)) + labs(x="Years",y=paste(indic," variation",sep="")) ##### AVEC INTERVAL ROMAIN
286 gg <- ggplot(ggdata,aes(x=year,y=estimate))+ geom_point() + geom_line() + labs(x="Years",y=paste(indicator," variation",sep="")) #####+ ylim(min(ggdata$se.inf),max(ggdata$se.sup))
287 gg <- gg + geom_line(size=1.5,colour=couleur)+ geom_point(size=3,colour=couleur) + geom_point(size=1.5,colour="white")
288 if (ic) {
289 gg <- gg + geom_ribbon(aes(ymin=ggdata$se.inf,ymax=ggdata$se.sup),alpha=.25,fill = couleur) + geom_errorbar(ymin=ggdata$se.inf,ymax=ggdata$se.sup,width=0,alpha=.25)
290 }
291 ggfile <- paste("Output/",indicator,"_glmmTMB_",id,".png",sep="")
292 ggsave(ggfile,gg)
293
294 ############ Estimation de la tendance sur la periode étudiée / Trends estimation on the time period studied
295 #cat("\nEstimation de la tendance lmer(",indicator,"~ year+(1|id_plot)\n",sep="")
296 cat("\nEstimation de la tendance glmmTMB(",indicator,"~ year+(1|id_plot)\n",sep="")
297 #md.c <- lmer(indic~ year+(1|id_plot),data=dd)##### effet aleatoire liés aux carrés sur l'ordonnée à l'origine / random effects of plots on intercept ### version lmer
298 md.c <- glmmTMB(indic~ year+(1|id_plot),data=dd)
299 smd.c<-summary(md.c)
300 # coefdata.c <- as.data.frame(smd.c$coefficients) #### pour la version lmer
301 coefdata.c <- as.data.frame(smd.c$coefficients$cond)[2,]
302 #profc=profile(md.c) ######### Ajout des intervalles de confiances / addition of the confidence intervals
303 if (ic) {
304 MODconfint=confint(md.c) ### plus rapide de ne pas passer par profile
305 se.inf=MODconfint[2,1]### [4,1] pour la version lmer
306 se.sup=MODconfint[2,2]### [4,2] pour la version lmer
307 } else{
308 se.inf <- "not assessed"
309 se.sup <- "not assessed"
310 }
311 coefdata.c <- data.frame(model = "Linear trend", variable = rownames(coefdata.c),coefdata.c,se.inf,se.sup)
312 coefdata <- rbind(coefdata.c,coefdata.f)
313 write.table(coefdata,paste("Output/GlmmTMB_coefficient_",indicator,id,".tabular",sep=""),row.names=FALSE,sep="\t")
314 write.table(ggdata,paste("Output/ggdata_",indicator,id,".tabular",sep=""),row.names=FALSE,sep="\t")
315 smd.file <- paste("Output/summary_lmer_",indicator,"_",id,".txt",sep="")
316 #####################
317 }
318 }
319
320
321
322
323
324 ##################
325 ### Do your analysis
326
327 csi_cti_ctri(tabCLEAN=tabCLEAN,coordCarre=coordCarre,spTrait=spTrait,dd=NULL,Var=Var,indicator=indicator,ic=ic,plot_smooth = plot_smooth,methode=methode)#,init_1989 = FALSE) ##### exemple pour l'indicateur csi sans csi déjà calculé donc à partir du ssi avec interval de confiance et utilisant modele mixte / example for the csi index which is not already calculated from the ssi with confidence interval using the mixte model
328