diff Moyenne_geom.r @ 0:985f8839aebd draft default tip

planemo upload for repository https://github.com/galaxyecology/tools-ecology/tree/master/tools/Geom_mean_workflow commit 3f11e193fd9ba5bf0c706cd5d65d6398166776cb

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/Moyenne_geom.r	Sat Nov 25 15:18:01 2023 +0000
@@ -0,0 +1,156 @@
+#### loading required R libraries
+#### chargement des packages R utilisés
+library(gdata)
+library(XLConnect)
+library(rms)
+
+###### overall parameters and settings
+###### paramètres globaux utilisés
+
+args = commandArgs(trailingOnly=TRUE) 
+if (length(args)==0)
+{
+    stop("This tool needs at least one argument")
+}else{
+    data <- args[1]
+    sep <- args[2]
+    HR <- args[3]
+     
+}
+
+if (HR =="false"){HR<-FALSE} else {HR<-TRUE}
+
+###nrep: number of samples used to calculate geometric means
+###nrep: nombre d'échantillons utilisés pour calculer les moyennes géométriques
+nrep<-10000
+
+#______________________________________________________________________________________________________________________________________________________________________________________________
+###### common functions
+###### fonction utiles pour la suite
+
+		convert.to.numeric<-function(x){
+		t(apply(x,1,function(x){as.double(sub(" ","",as.character(x)))}))}
+
+		
+		### calculus of the logarithm of nrep geometric means, sampling based on a lognormal distribution with the same moments as the empirical ones (means & Ics)
+			#to prevent negative values
+		### calcul du logarithme de nrep moyennes géométriques, l'échantillonnage étant fait avec la distribution lognormale de mêmes moments que les momenst empriques (means et ICs)
+			#pour éviter d'avoir des valeurs négatives
+
+		lgeomean<-function(means,ICs,nrep)
+		{#means: vector: mean estimates for the different categories 
+		#ICs: vector: in proportion to the mean, difference between the extremum of the 95% confidence interval and the mean
+		require(mvtnorm)
+		#calculation of the parameters of the log normal distribution (on the log scale)
+		#cf. http://127.0.0.1:26338/library/stats/html/Lognormal.html
+		logsigma<-sqrt(log((ICs/qnorm(0.975)/means)^2+1))
+		logmean<-log(means)-1/2*logsigma^2
+
+		#gaussian sampling on the log scale then taking exponential
+		temp<-exp(rmvnorm(nrep,mean=logmean,sigma=diag(logsigma*logsigma)))
+
+		#taking geometric mean over categories, but kept on the log scale
+		geomm.rep<-apply(temp,1,function(x){(mean(log(x),na.rm=TRUE))})
+		#c(mean(geomm.rep),sd(geomm.rep))
+		geomm.rep}
+#_______________________________________________________________________________________________________________________________________________________________________________________________
+
+###### importation des données
+###### importation of data
+temp<-read.csv(file=data,sep=sep,header=HR,encoding="UTF-8")
+
+data2008_2012<-temp[4:14,]
+data2013_2017<-temp[21:31,]
+
+meandata2008_2012<-convert.to.numeric(data2008_2012[,c(3,6,9)])
+ICdata2008_2012<-convert.to.numeric(data2008_2012[,c(5,8,11)])
+meandata2013_2017<-convert.to.numeric(data2013_2017[,c(3,6,9)])
+ICdata2013_2017<-convert.to.numeric(data2013_2017[,c(5,8,11)])
+
+####### code to calculate (nrep) logarithms of geometric means by region (Greco)
+####### code pour calculer les nrep logarithmes de moyennes géométriques par région (GRECO)
+
+set.seed(1)
+#first period
+#première période
+rest2008_2012<-sapply(1:dim(data2008_2012)[1],function(region){lgeomean(meandata2008_2012[region,],ICdata2008_2012[region,],nrep)})
+
+set.seed(3)
+#first period but with different seed
+#première période mais avec une graine différente
+rest2008_2012_s3<-sapply(1:dim(data2008_2012)[1],function(region){lgeomean(meandata2008_2012[region,],ICdata2008_2012[region,],nrep)})
+
+set.seed(2)
+#second period
+#seconde période
+rest2013_2017<-sapply(1:dim(data2013_2017)[1],function(region){lgeomean(meandata2013_2017[region,],ICdata2013_2017[region,],nrep)})
+
+
+####### code to summarize the above nrep logarithms of geometric means by region into the statistics of an overall geometric mean across regions, taking the first period as reference
+###### code pour passer des nrep logarithmes de moyenne géométrique par région aux statistiques de la moyenne géométrique globale, en prennat la première période comme référence
+
+#for the first period
+#pour la première période
+Mean_2008_2012_scaled<-{temp<-apply(rest2008_2012_s3,1,function(x){mean(x)})-apply(rest2008_2012,1,function(x){mean(x)});c(mean(exp(temp)),sd(exp(temp)),quantile(exp(temp),prob=c(0.025,0.975)))}
+
+#for the second period
+#pour la seconde période
+Mean_2013_2017_scaled<-{temp<-apply(rest2013_2017,1,function(x){mean(x)})-apply(rest2008_2012,1,function(x){mean(x)});c(mean(exp(temp)),sd(exp(temp)),quantile(exp(temp),prob=c(0.025,0.975)))}
+
+
+
+############### NATIONAL OUPUTS:
+############### SORTIES NATIONALES:
+
+res2008_2012_scaled_df = data.frame(Mean_2008_2012_scaled)
+res2008_2012_scaled_df=`rownames<-`(res2008_2012_scaled_df,c("mean","sd","2,5%","97,5%"))
+
+res2013_2017_scaled_df = data.frame(Mean_2013_2017_scaled)
+res2013_2017_scaled_df=`rownames<-`(res2013_2017_scaled_df,c("mean","sd","2,5%","97,5%"))
+
+
+write.csv(res2008_2012_scaled_df, file = "res2008_2012_scaled.csv")
+write.csv(res2013_2017_scaled_df,file= "res2013_2017_scaled.csv")
+
+############### REGIONAL OUPUTS:
+############### SORTIES REGIONALES (GRECO):
+
+regres2008_2012_scaled<-apply(rest2008_2012_s3-rest2008_2012,2,function(x){temp<-x;c(mean=mean(exp(temp)),sd=sd(exp(temp)),quantile(exp(temp),prob=c(0.025,0.975)))})
+regres2013_2017_scaled<-apply(rest2013_2017-rest2008_2012,2,function(x){temp<-x;c(mean=mean(exp(temp)),sd=sd(exp(temp)),quantile(exp(temp),prob=c(0.025,0.975)))})
+dimnames(regres2008_2012_scaled)[[2]]<-as.character(data2008_2012[,2])
+dimnames(regres2013_2017_scaled)[[2]]<-as.character(data2013_2017[,2])
+
+write.csv(regres2008_2012_scaled, file = "regres2008_2012_scaled.csv")
+write.csv(regres2013_2017_scaled, file = "regres2013_2017_scaled.csv")
+
+############### data to make a bar plot of the national evolution rate 
+histo_data = data.frame(
+  variable_name = c(names(res2008_2012_scaled_df),names(res2013_2017_scaled_df)), 
+  variable = c(round(Mean_2008_2012_scaled[1]*100),round(Mean_2013_2017_scaled[1]*100)),
+  standard_deviation = c(Mean_2008_2012_scaled[2]*100,Mean_2013_2017_scaled[2]*100)
+)
+
+write.table(histo_data, file = "histo_data.tsv",row.names = F, col.names = T ,sep ="\t")
+
+############### data to make a map of the GRECO evolution rate
+
+rate2008_2012 = data.frame(round(regres2008_2012_scaled[1,1:11]*100))
+rate2013_2017 = data.frame(round(regres2013_2017_scaled[1,1:11]*100))
+
+evol_rate = rate2013_2017-rate2008_2012
+evol_rate = cbind(data2013_2017[,2],evol_rate)
+colnames(evol_rate)<-c("Regions","Evolution_rate")
+
+
+write.table(evol_rate,"evolution_rate.tsv",sep="\t",quote=F,row.names=F,col.names=T)
+
+
+
+
+
+
+
+
+
+
+