Mercurial > repos > ecology > pampa_glmcomm
comparison FunctExeCalcGLMGalaxy.r @ 4:44b9069775ca draft
"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") |