Mercurial > repos > ecology > pampa_glmcomm

comparison FunctExeCalcGLMGalaxy.r @ 4:44b9069775ca draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
"planemo upload for repository https://github.com/ColineRoyaux/PAMPA-Galaxy commit 3df1978827a91be30e815dee2ed83a92862d1b1c"
author ecology
date Sun, 22 Nov 2020 18:40:23 +0000
parents f0dc3958e65d
children 8db02073b671
comparison
equal deleted inserted replaced
3:97967fbcf6da 4:44b9069775ca
1 #Rscript 1 #Rscript
2 2
3 ##################################################################################################################### 3 #####################################################################################################################
4 ##################################################################################################################### 4 #####################################################################################################################
5 ################################# Compute a Generalized Linear Model from your data ################################# 5 ################################# Compute a Generalized Linear Model from your data #################################
6 ##################################################################################################################### 6 #####################################################################################################################
7 ##################################################################################################################### 7 #####################################################################################################################
8 8
9 ###################### Packages 9 ###################### Packages
10 suppressMessages(library(multcomp)) 10 suppressMessages(library(multcomp))
11 suppressMessages(library(DHARMa))
11 suppressMessages(library(glmmTMB)) ###Version: 0.2.3 12 suppressMessages(library(glmmTMB)) ###Version: 0.2.3
12 suppressMessages(library(gap)) 13 suppressMessages(library(gap))
14
13 15
14 ###################### Load arguments and declaring variables 16 ###################### Load arguments and declaring variables
15 17
16 args = commandArgs(trailingOnly=TRUE) 18 args <- commandArgs(trailingOnly = TRUE)
17 #options(encoding = "UTF-8")
18 19
19 if (length(args) < 10) { 20 if (length(args) < 10) {
20 stop("At least 4 arguments must be supplied : \n- two input dataset files (.tabular) : metrics table and unitobs table \n- Interest variable field from metrics table \n- Response variable from unitobs table.", call.=FALSE) #si pas d'arguments -> affiche erreur et quitte / if no args -> error and exit1 21 stop("At least 4 arguments must be supplied : \n- two input dataset files (.tabular) : metrics table and unitobs table \n- Interest variable field from metrics table \n- Response variable from unitobs table.", call. = FALSE) # if no args -> error and exit1
21 22
22 } else { 23 } else {
23 Importdata <- args[1] ###### file name : metrics table 24 import_data <- args[1] ###### file name : metrics table
24 ImportUnitobs <- args[2] ###### file name : unitobs informations 25 import_unitobs <- args[2] ###### file name : unitobs informations
25 colmetric <- as.numeric(args[3]) ###### Selected interest metric for GLM 26 colmetric <- as.numeric(args[3]) ###### Selected interest metric for GLM
26 listFact <- strsplit(args [4],",")[[1]] ###### Selected response factors for GLM 27 list_fact <- strsplit(args [4], ",")[[1]] ###### Selected response factors for GLM
27 listRand <- strsplit(args [5],",")[[1]] ###### Selected randomized response factors for GLM 28 list_rand <- strsplit(args [5], ",")[[1]] ###### Selected randomized response factors for GLM
28 colFactAna <- args[6] ####### (optional) Selected splitting factors for GLMs 29 col_fact_ana <- args[6] ####### (optional) Selected splitting factors for GLMs
29 Distrib <- args[7] ###### (optional) Selected distribution for GLM 30 distrib <- args[7] ###### (optional) Selected distribution for GLM
30 log <- args[8] ###### (Optional) Log on interest metric ? 31 glm_out <- args[8] ###### (Optional) GLM object as Rdata output ?
31 aggreg <- args[9] ###### Aggregation level of the data table 32 aggreg <- args[9] ###### Aggregation level of the data table
32 source(args[10]) ###### Import functions 33 source(args[10]) ###### Import functions
33 34
34 } 35 }
35 #### Data must be a dataframe with at least 3 variables : unitobs representing location and year ("observation.unit"), species code ("species.code") and abundance ("number") 36
36 37
37 38 #### d_ata must be a dataframe with at least 3 variables : unitobs representing location and year ("observation.unit"), species code ("species.code") and abundance ("number")
38 #Import des données / Import data 39
39 obs<- read.table(Importdata,sep="\t",dec=".",header=TRUE,encoding="UTF-8") # 40
40 obs[obs == -999] <- NA 41 #Import des données / Import data
42 obs <- read.table(import_data, sep = "\t", dec = ".", header = TRUE, encoding = "UTF-8") #
43 obs[obs == -999] <- NA
41 metric <- colnames(obs)[colmetric] 44 metric <- colnames(obs)[colmetric]
42 tabUnitobs <- read.table(ImportUnitobs,sep="\t",dec=".",header=TRUE,encoding="UTF-8") 45 tab_unitobs <- read.table(import_unitobs, sep = "\t", dec = ".", header = TRUE, encoding = "UTF-8")
43 tabUnitobs[tabUnitobs == -999] <- NA 46 tab_unitobs[tab_unitobs == -999] <- NA
44 47
45 if (colFactAna != "None") 48 if (col_fact_ana != "None") {
46 { 49 fact_ana <- colnames(tab_unitobs)[as.numeric(col_fact_ana)]
47 FactAna <- colnames(tabUnitobs)[as.numeric(colFactAna)] 50 if (class(tab_unitobs[fact_ana]) == "numeric" || fact_ana == "observation.unit") {
48 if (class(tabUnitobs[FactAna]) == "numeric" || FactAna == "observation.unit"){stop("Wrong chosen separation factor : Analysis can't be separated by observation unit or numeric factor")} 51 stop("Wrong chosen separation factor : Analysis can't be separated by observation unit or numeric factor")
52 }
49 }else{ 53 }else{
50 FactAna <- colFactAna 54 fact_ana <- col_fact_ana
51 } 55 }
52 56
53 57 vars_data1 <- NULL
54 #factors <- fact.det.f(Obs=obs) 58 err_msg_data1 <- "The input metrics dataset doesn't have the right format. It needs to have at least the following 2 variables :\n- observation.unit (or year and site)\n- numeric or integer metric\n"
55 59 check_file(obs, err_msg_data1, vars_data1, 2)
56 vars_data1<- NULL 60
57 err_msg_data1<-"The input metrics dataset doesn't have the right format. It needs to have at least the following 2 variables :\n- observation.unit (or year and site)\n- numeric or integer metric\n" 61 vars_data2 <- c("observation.unit", list_fact, list_rand)
58 check_file(obs,err_msg_data1,vars_data1,2) 62 err_msg_data2 <- "The input unitobs dataset doesn't have the right format. It needs to have at least the following 2 variables :\n- observation.unit (or year and site)\n- factors used in GLM (habitat, year and/or site)\n"
59 63 check_file(tab_unitobs, err_msg_data2, vars_data2[vars_data2 != "None"], 2)
60 vars_data2 <- c(listFact,listRand) 64
61 err_msg_data2<-"The input unitobs dataset doesn't have the right format. It needs to have at least the following 2 variables :\n- observation.unit (or year and site)\n- factors used in GLM (habitat, year and/or site)\n" 65 if (all(c(list_fact, list_rand) == "None")) {
62 check_file(tabUnitobs,err_msg_data2,vars_data2[vars_data2 != "None"],2) 66 stop("GLM needs to have at least one response variable.")
67 }
68
69 if (list_fact[1] == "None" || all(is.element(list_fact, list_rand))) {
70 stop("GLM can't have only random effects.")
71 }
72
63 73
64 #################################################################################################### 74 ####################################################################################################
65 ########## Computing Generalized Linear Model ## Function : modeleLineaireWP2.unitobs.f ############ 75 ######### Computing Generalized Linear Model ## Function : glm_community ############
66 #################################################################################################### 76 ####################################################################################################
67 77
68 modeleLineaireWP2.unitobs.f <- function(metrique, listFact, listRand, FactAna, Distrib, log=FALSE, tabMetrics, tableMetrique, tabUnitobs, unitobs="observation.unit", nbName="number") 78 glm_community <- function(metrique, list_fact, list_rand, fact_ana, distrib, tab_metrics, tab_metrique, tab_unitobs, unitobs = "observation.unit", nb_name = "number") {
69 { 79 ## Purpose: Monitoring steps for GLM on community data
70 ## Purpose: Monitoring steps for GLM on unitobs
71 ## ---------------------------------------------------------------------- 80 ## ----------------------------------------------------------------------
72 ## Arguments: metrique : selected metric 81 ## Arguments: metrique : selected metric
73 ## listFact : Factors for GLM 82 ## list_fact : Factors for GLM
74 ## listRand : Random factors for GLM 83 ## list_rand : Random factors for GLM
75 ## factAna : Separation factor for GLMs 84 ## fact_ana : Separation factor for GLMs
76 ## Distrib : selected distribution for model 85 ## distrib : selected distribution for model
77 ## log : log transformation on data ? boolean 86 ## tab_metrics : data table metrics
78 ## tabMetrics : data table metrics 87 ## tab_metrique : data table's name
79 ## tableMetrique : data table's name 88 ## tab_unitobs : data table unitobs
80 ## tabUnitobs : data table unitobs
81 ## ---------------------------------------------------------------------- 89 ## ----------------------------------------------------------------------
82 ## Author: Yves Reecht, Date: 18 août 2010, 15:59 modified by Coline ROYAUX 04 june 2020 90 ## Author: Yves Reecht, Date: 18 août 2010, 15:59 modified by Coline ROYAUX 04 june 2020
83 91
84 tmpData <- tabMetrics 92 tmpd_ata <- tab_metrics
85 93
86 if (listRand[1] != "None") 94 out_fact <- .GlobalEnv$organise_fact(list_rand = list_rand, list_fact = list_fact)
87 { 95 resp_fact <- out_fact[[1]]
88 if (all(is.element(listFact,listRand)) || listFact[1] == "None") 96 list_f <- out_fact[[2]]
89 { 97 list_fact <- out_fact[[3]]
90 RespFact <- paste("(1|",paste(listRand,collapse=") + (1|"),")")
91 listF <- NULL
92 listFact <- listRand
93 }else{
94 listF <- listFact[!is.element(listFact,listRand)]
95 RespFact <- paste(paste(listF, collapse=" + ")," + (1|",paste(listRand,collapse=") + (1|"),")")
96 listFact <- c(listF,listRand)
97 }
98 }else{
99 listF <- listFact
100 RespFact <- paste(listFact, collapse=" + ")
101 }
102 98
103 ##Creating model's expression : 99 ##Creating model's expression :
104 100 expr_lm <- eval(parse(text = paste(metrique, "~", resp_fact)))
105 #if (log == FALSE) {
106 exprML <- eval(parse(text=paste(metrique, "~", RespFact)))
107 #}else{
108 # exprML <- eval(parse(text=paste("log(",metrique,")", "~", RespFact)))
109 #}
110 101
111 ##Creating analysis table : 102 ##Creating analysis table :
112 103
113 listFactTab <- c(listFact, FactAna) 104 list_fact_tab <- c(list_fact, fact_ana)
114 listFactTab <- listFactTab[listFactTab != "None"] 105 list_fact_tab <- list_fact_tab[list_fact_tab != "None"]
115 106
116 if (all(is.na(match(tmpData[,unitobs],tabUnitobs[,unitobs])))) {stop("Observation units doesn't match in the two input tables")} 107 if (all(is.na(match(tmpd_ata[, unitobs], tab_unitobs[, unitobs])))) {
117 108 stop("Observation units doesn't match in the two input tables")
118 if(! is.element("species.code",colnames(tmpData))) 109 }
119 { 110
120 col <- c(unitobs,metrique) 111 if (! is.element("species.code", colnames(tmpd_ata))) {
121 tmpData <- cbind(tmpData[,col], tabUnitobs[match(tmpData[,unitobs],tabUnitobs[,unitobs]),listFactTab]) 112 col <- c(unitobs, metrique)
122 colnames(tmpData) <- c(col,listFactTab) 113 tmpd_ata <- cbind(tmpd_ata[, col], tab_unitobs[match(tmpd_ata[, unitobs], tab_unitobs[, unitobs]), list_fact_tab])
123 114 colnames(tmpd_ata) <- c(col, list_fact_tab)
124 for (i in listFactTab) { 115
116 for (i in list_fact_tab) {
125 switch(i, 117 switch(i,
126 tmpData[,i] <- as.factor(tmpData[,i])) 118 tmpd_ata[, i] <- as.factor(tmpd_ata[, i]))
127 } 119 }
128 }else{ 120 }else{
129 stop("Warning : wrong data frame, data frame should be aggregated by observation unit (year and site)") 121 stop("Warning : wrong data frame, data frame should be aggregated by observation unit (year and site)")
130 } 122 }
131 123
132 ## Suppress unsed 'levels' : 124 ## Suppress unsed 'levels' :
133 tmpData <- dropLevels.f(tmpData) 125 tmpd_ata <- .GlobalEnv$drop_levels_f(tmpd_ata)
134 126
135 ## Automatic choice of distribution if none is selected by user : 127 ## Automatic choice of distribution if none is selected by user :
136 if (Distrib == "None") 128
137 { 129 chose_distrib <- .GlobalEnv$distrib_choice(distrib = distrib, metrique = metrique, data = tmpd_ata)
138 switch(class(tmpData[,metrique]), 130
139 "integer"={loiChoisie <- "poisson"}, 131 if (fact_ana != "None" && nlevels(tmpd_ata[, fact_ana]) > 1) {
140 "numeric"={loiChoisie <- "gaussian"}, 132 ana_cut <- levels(tmpd_ata[, fact_ana])
141 stop("Selected metric class doesn't fit, you should select an integer or a numeric variable")) 133 }else{
142 }else{ 134 ana_cut <- NULL
143 loiChoisie <- Distrib 135 }
144 } 136
145 137 ##Create results table :
146 138 lev <- unlist(lapply(list_f, FUN = function(x) {
147 ## Compute Model(s) : 139 levels(tmpd_ata[, x])
148 140 }))
149 if (FactAna != "None" && nlevels(tmpData[,FactAna]) > 1) 141 row <- c("global", ana_cut)
150 { 142
151 Anacut <- levels(tmpData[,FactAna]) 143 if (is.element("year", list_f) && ! is.element("year", list_rand)) {
152 }else{ 144 tab_sum <- .GlobalEnv$create_res_table(list_rand = list_rand, list_fact = list_fact, row = row, lev = unlist(c("year", lev)), distrib = chose_distrib)
153 Anacut <- NULL 145 }else{
154 } 146 tab_sum <- .GlobalEnv$create_res_table(list_rand = list_rand, list_fact = list_fact, row = row, lev = lev, distrib = chose_distrib)
155 147 }
156 ##Create results table : 148 ### creating rate table
157 lev <- unlist(lapply(listF,FUN=function(x){levels(tmpData[,x])})) 149 tab_rate <- data.frame(analysis = row, complete_plan = NA, balanced_plan = NA, NA_proportion_OK = NA, no_residual_dispersion = NA, uniform_residuals = NA, outliers_proportion_OK = NA, no_zero_inflation = NA, observation_factor_ratio_OK = NA, enough_levels_random_effect = NA, rate = NA)
158 150
159 if (listRand[1] != "None") ## if random effects 151 ##plural analysis
160 { 152 for (cut in ana_cut) {
161 TabSum <- data.frame(analysis=c("global", Anacut),AIC=NA,BIC=NA,logLik=NA, deviance=NA,df.resid=NA) 153 cutd_ata <- tmpd_ata[grep(cut, tmpd_ata[, fact_ana]), ]
162 colrand <- unlist(lapply(listRand, 154 cutd_ata <- .GlobalEnv$drop_levels_f(cutd_ata)
163 FUN=function(x){lapply(c("Std.Dev","NbObservation","NbLevels"), 155
164 FUN=function(y){paste(x,y,collapse = ":") 156 res <- ""
165 }) 157 resy <- ""
166 })) 158
167 TabSum[,colrand] <- NA 159 if (list_rand[1] != "None") {
168 160 res <- tryCatch(glmmTMB::glmmTMB(expr_lm, family = chose_distrib, data = cutd_ata), error = function(e) {
169 if (! is.null(lev)) ## if fixed effects + random effects 161 })
170 { 162 if (is.element("year", list_f) && ! is.element("year", list_rand)) { #Model with year as continuous
171 colcoef <- unlist(lapply(c("(Intercept)",lev), 163 cutd_ata$year <- as.numeric(cutd_ata$year)
172 FUN=function(x){lapply(c("Estimate","Std.Err","Zvalue","Pvalue","signif"), 164 resy <- tryCatch(glmmTMB::glmmTMB(expr_lm, family = chose_distrib, data = cutd_ata), error = function(e) {
173 FUN=function(y){paste(x,y,collapse = ":") 165 })
174 }) 166 cutd_ata$year <- as.factor(cutd_ata$year)
175 })) 167 }else{
176 }else{ ## if no fixed effects 168 resy <- ""
177 colcoef <- NULL 169 }
178 } 170
179 171 }else{
180 }else{ ## if no random effects 172 res <- tryCatch(glm(expr_lm, data = cutd_ata, family = chose_distrib), error = function(e) {
181 TabSum <- data.frame(analysis=c("global", Anacut),AIC=NA,Resid.deviance=NA,df.resid=NA,Null.deviance=NA,df.null=NA) 173 })
182 174 if (is.element("year", list_f)) { #Model with year as continuous
183 switch(loiChoisie, 175 cutd_ata$year <- as.numeric(cutd_ata$year)
184 "gaussian"={colcoef <- unlist(lapply(c("(Intercept)",lev), 176 resy <- tryCatch(glm(expr_lm, family = chose_distrib, data = cutd_ata), error = function(e) {
185 FUN=function(x){lapply(c("Estimate","Std.Err","Tvalue","Pvalue","signif"), 177 })
186 FUN=function(y){paste(x,y,collapse = ":") 178 cutd_ata$year <- as.factor(cutd_ata$year)
187 }) 179 }else{
188 }))}, 180 resy <- ""
189 "quasipoisson"={colcoef <- unlist(lapply(c("(Intercept)",lev), 181 }
190 FUN=function(x){lapply(c("Estimate","Std.Err","Tvalue","Pvalue","signif"), 182
191 FUN=function(y){paste(x,y,collapse = ":")
192 })
193 }))},
194 colcoef <- unlist(lapply(c("(Intercept)",lev),
195 FUN=function(x){lapply(c("Estimate","Std.Err","Zvalue","Pvalue","signif"),
196 FUN=function(y){paste(x,y,collapse = ":")
197 })
198 })))
199
200 }
201
202 TabSum[,colcoef] <- NA
203
204 ### creating rate table
205 TabRate <- data.frame(analysis=c("global", Anacut), complete_plan=NA, balanced_plan=NA, NA_proportion_OK=NA, no_residual_dispersion=NA, uniform_residuals=NA, outliers_proportion_OK=NA, no_zero_inflation=NA, observation_factor_ratio_OK=NA, enough_levels_random_effect=NA, rate=NA)
206
207 for (cut in Anacut)
208 {
209 cutData <- tmpData[grep(cut,tmpData[,FactAna]),]
210 cutData <- dropLevels.f(cutData)
211
212 res <-""
213
214 if (listRand[1] != "None")
215 {
216 res <- tryCatch(glmmTMB(exprML,family=loiChoisie, data=cutData), error=function(e){})
217 }else{
218 res <- tryCatch(glm(exprML,data=cutData,family=loiChoisie), error=function(e){})
219 } 183 }
220 184
221 ## Write results : 185 ## Write results :
222 if (! is.null(res)) 186 if (! is.null(res)) {
223 { 187 file_save_glm_cut <- paste("GLM_", cut, ".Rdata", sep = "")
224 TabSum <- sortiesLM.f(objLM=res, TabSum=TabSum, metrique=metrique, 188 save(res, file = file_save_glm_cut)
225 factAna=factAna, cut=cut, colAna="analysis", lev=lev, #modSel=iFactGraphSel, listFactSel=listFactSel, 189
226 listFact=listFact, 190 tab_sum <- sorties_lm_f(obj_lm = res, obj_lmy = resy, tab_sum = tab_sum, metrique = metrique,
227 Data=cutData, #Log=Log, 191 fact_ana = fact_ana, cut = cut, col_ana = "analysis", lev = lev, #modSel = iFactGraphSel, list_fact_sel = list_fact_sel,
228 type=ifelse(tableMetrique == "unitSpSz" && factAna != "size.class", 192 list_fact = list_fact, list_rand = list_rand,
229 "CL_unitobs", 193 d_ata = cutd_ata)
230 "unitobs")) 194
231 195 tab_rate[tab_rate[, "analysis"] == cut, c(2:11)] <- .GlobalEnv$note_glm_f(data = cutd_ata, obj_lm = res, metric = metrique, list_fact = list_fact, details = TRUE)
232 TabRate[TabRate[,"analysis"]==cut,c(2:11)] <- noteGLM.f(data=cutData, objLM=res, metric=metrique, listFact=listFact, details=TRUE) 196
233 197 }else{
234 }else{ 198 cat("\nCannot compute GLM for level", cut, "Check if one or more factor(s) have only one level, or try with another distribution for the model in advanced settings \n\n")
235 cat("\nCannot compute GLM for level",cut,"Check if one or more factor(s) have only one level, or try with another distribution for the model in advanced settings \n\n") 199 }
236 } 200
237 201 }
238 } 202
239 203 ## Global analysis :
240 ## Global analysis : 204
241 205 res_g <- ""
242 if (listRand[1] != "None") 206 res_gy <- ""
243 { 207
244 resG <- glmmTMB(exprML,family=loiChoisie, data=tmpData) 208 if (list_rand[1] != "None") {
245 }else{ 209 res_g <- glmmTMB::glmmTMB(expr_lm, family = chose_distrib, data = tmpd_ata)
246 resG <- glm(exprML,data=tmpData,family=loiChoisie) 210 if (is.element("year", list_fact) && ! is.element("year", list_rand)) { #Model with year as continuous
211 yr <- tmpd_ata$year
212 tmpd_ata$year <- as.numeric(tmpd_ata$year)
213 res_gy <- glmmTMB::glmmTMB(expr_lm, family = chose_distrib, data = tmpd_ata)
214 tmpd_ata$year <- as.factor(tmpd_ata$year)
215 tmpd_ata$year <- yr
216 }else{
217 res_gy <- ""
218 }
219
220 }else{
221 res_g <- glm(expr_lm, data = tmpd_ata, family = chose_distrib)
222 if (is.element("year", list_fact)) { #Model with year as continuous
223 yr <- tmpd_ata$year
224 tmpd_ata$year <- as.numeric(tmpd_ata$year)
225 res_gy <- glm(expr_lm, family = chose_distrib, data = tmpd_ata)
226 tmpd_ata$year <- yr
227 }else{
228 res_gy <- ""
229 }
247 } 230 }
248 231
249 ## write results : 232 ## write results :
250 TabSum <- sortiesLM.f(objLM=resG, TabSum=TabSum, metrique=metrique, 233
251 factAna=factAna, cut="global", colAna="analysis", lev=lev, #modSel=iFactGraphSel, listFactSel=listFactSel, 234 save(res_g, file = "global_GLM.Rdata")
252 listFact=listFact, 235
253 Data=tmpData, #Log=Log, 236 tab_sum <- sorties_lm_f(obj_lm = res_g, obj_lmy = res_gy, tab_sum = tab_sum, metrique = metrique,
254 type=ifelse(tableMetrique == "unitSpSz" && factAna != "size.class", 237 fact_ana = fact_ana, cut = "global", col_ana = "analysis", lev = lev, #modSel = iFactGraphSel, list_fact_sel = list_fact_sel,
255 "CL_unitobs", 238 list_fact = list_fact, list_rand = list_rand,
256 "unitobs")) 239 d_ata = tmpd_ata)
257 240
258 TabRate[TabRate[,"analysis"]=="global",c(2:11)] <- noteGLM.f(data=tmpData, objLM=resG, metric=metrique, listFact=listFact, details=TRUE) 241 tab_rate[tab_rate[, "analysis"] == "global", c(2:11)] <- .GlobalEnv$note_glm_f(data = tmpd_ata, obj_lm = res_g, metric = metrique, list_fact = list_fact, details = TRUE)
259 noteGLMs.f(tabRate=TabRate,exprML=exprML,objLM=resG, file_out=TRUE) 242 .GlobalEnv$note_glms_f(tab_rate = tab_rate, expr_lm = expr_lm, obj_lm = res_g, file_out = TRUE)
243
260 ## simple statistics and infos : 244 ## simple statistics and infos :
261 filename <- "GLMSummaryFull.txt" 245 filename <- "GLMSummaryFull.txt"
262 246 info_stats_f(filename = filename, d_ata = tmpd_ata, agreg_level = aggreg, type = "stat",
263 ## Save data on model : 247 metrique = metrique, fact_graph = fact_ana, #fact_graph_sel = modSel,
264 248 list_fact = list_fact)#, list_fact_sel = list_fact_sel)
265 infoStats.f(filename=filename, Data=tmpData, agregLevel=aggreg, type="stat", 249
266 metrique=metrique, factGraph=factAna, #factGraphSel=modSel, 250 tab_sum$separation <- fact_ana
267 listFact=listFact)#, listFactSel=listFactSel) 251
268 252 return(tab_sum)
269 return(TabSum)
270 253
271 } 254 }
272 255
273 ################# Analysis 256 ################# Analysis
274 257
275 Tab <- modeleLineaireWP2.unitobs.f(metrique=metric, listFact=listFact, listRand=listRand, FactAna=FactAna, Distrib=Distrib, log=log, tabMetrics=obs, tableMetrique=aggreg, tabUnitobs=tabUnitobs, nbName="number") 258 tab <- glm_community(metrique = metric, list_fact = list_fact, list_rand = list_rand, fact_ana = fact_ana, distrib = distrib, tab_metrics = obs, tab_metrique = aggreg, tab_unitobs = tab_unitobs, nb_name = "number")
276 259
277 write.table(Tab,"GLMSummary.tabular", row.names=FALSE, sep="\t", dec=".",fileEncoding="UTF-8") 260 write.table(tab, "GLMSummary.tabular", row.names = FALSE, sep = "\t", dec = ".", fileEncoding = "UTF-8")