##############################################
#
# Code de
propagation de l'incertitude pour l'analyse du risque de niveau 2 de toutes les
CIE
# Élaboré aux
fins d'utilisation avec R par Gerald Singh
# 2 Mai 2014
#
-------------------------------------------
#
# Script
généralisé pour tous les types de CIE (espèce, habitat et propriétés des
communautés/des écosystèmes) et
# efficacité
améliorée, dans la mesure du possible. Il convient de noter que d'autres
ajustements devront être
# apportés à la
mise en forme des figures afin de refléter le nombre variable de CIE, etc.
#
# Lors de
l'application à une nouvelle évaluation du risque, il convient de noter que des
ajustements peuvent être
# requis pour
affiner la mise en forme des figures propres au CERE (après
# la ligne
370).
#
# For more
information: M. Thiess
#
mary.thiess@dfo-mpo.gc.ca
#
# Références:
# [1] Murray,
C. C., Mach, M.E. et O, M. (2016). Évaluation pilote des risques pour
l'écosystème en vue d'évaluer
# les
risques cumulatifs pour les espèces dans la zone de gestion intégrée de la côte
nord du Pacifique.
# Secr.
can. de consult. sci. du MPO. Doc. de rech. 2016/049. vii + 59 p.
#
##############################################
# Nouvel espace de travail
rm(list =
ls(all=TRUE))
# Ouvrir les trousses nécessaires
require(psych)
require(plyr)
require(dtplyr)
require(data.table)
###############################################
# PARAMÈTRES DÉFINIS PAR L'UTILISATEUR
###############################################
# Définir le nombre d'échantillons aléatoires à
utiliser pour produire des estimations d'incertitude:
nreps = 10000
# Sélectionner la distribution à utiliser pour obtenir
des échantillons aléatoires:
distn =
"normal" # Can be
"trunc.norm", "normal", or "multinomial"
# Pour reproduire exactement les résultats du document
de recherche, définir le point de départ
# du générateur de nombres aléatoires:
#set.seed(101)
# Créer des figures récapitulatives?
do.figures =
"Y" # Can be "Y"
or "N"
# Variables de CIE, d'activités et d'agents de stress
interaction.nm =
c("SEC.Type","SEC","Activity","Sub.Activity","Stressor")
# Variables
d'exposition
exp.nm =
c("Area","Depth","Temporal","Intensity","Intensity2")
# Variables de résilience
res.nm =
c("AcuteChange","ChronicChange")
# Variables de rétablissement
rec.nm =
c("Fecundity","NatMort","AgeMat","LifeStage","PopConn","Listed","MaxAge","MaxSize","vonBert")
###############################################
# FONCTIONS DÉFINIES PAR L'UTILISATEUR
###############################################
# Sélectionne un sous-ensemble de facteurs de risque
dans la trame de données saisie:
risk.subset =
function(x,df.type){
t.score =
matrix(unlist(df.type[,x]),ncol=length(x),byrow=FALSE)
dimnames(t.score)[[2]] = x
return(t.score)
}
# Sélectionne un échantillon aléatoire (n=nreps) à
l'aide d'une distribution normale tronquée.
# des valeurs aux éléments du système de cotation qui
se trouvent hors des limites supérieure
# et inférieure:
norm.fn = function(mean.mat, sd.mat,lwr.bd,upr.bd){
output.rand =
array(NA,dim=c(no.Int,ncol(mean.mat),nreps))
for (k in 1:nreps){
output.rand[,,k] = mean.mat
+
matrix(rnorm(no.Int*ncol(mean.mat),0,sd.mat),no.Int,ncol(mean.mat))
}
#
Réattribuer des valeurs aux éléments hors de la plage de cotes:
output.rand[output.rand<lwr.bd] = lwr.bd
output.rand[output.rand>upr.bd] =
upr.bd
return(output.rand)
}
# Génère une distribution normale tronquée délimitée
par a et b:
rtnorm <-
function(n, mean, sd, a = -Inf, b = Inf){
qnorm(runif(n, pnorm(a, mean, sd),
pnorm(b, mean, sd)), mean, sd)
}
# Sélectionne un échantillon aléatoire (n=nreps) à
l'aide d'une distribution normale tronquée:
trnorm.fn = function(mean.mat,sd.mat,lwr.bd,upr.bd){
output.rand = array(NA,dim=c(no.Int,ncol(mean.mat),nreps))
for (k in 1:nreps){
output.rand[,,k] =
matrix(rtnorm(no.Int*ncol(mean.mat),mean.mat,sd.mat,lwr.bd,upr.bd),
no.Int,ncol(mean.mat),byrow=FALSE)
}
return(output.rand)
}
#################################
# Préliminaires
#################################
# Ouvrir le fichier contenant les colonnes relatives
aux CIE, activités, agents de stress et cotes des facteurs de risque:
imp =
read.csv(file.choose(), stringsAsFactors = FALSE, strip.white = TRUE)
#Créer un dossier pour y conserver les fichiers de
résultats:
d =
paste("OUTPUT_",Sys.Date(),sep="")
dir.create(d)
# Vérifier si des variables sont manquantes (définies
dans les PARAMÈTRES DÉFINIS PAR L'UTILISATEUR) dans le fichier d'entrée.
# Interrompre le processus si des variables sont
manquantes:
allvars <-
c(exp.nm, res.nm, rec.nm, interaction.nm)
if( !all(allvars
%in% names(imp)) ){
notfound <- allvars[!allvars %in%
names(imp)]
stop( "ERROR: Column name(s) '",
paste( notfound, collapse="', '" ), "' not found in input
file." )
}
# Exclure toutes les combinaisons CIE-Activité-Agent
de stress dont la cote est 0 pour ChangementAigu
# et ChangementChronique (Resilience components):
imp1 = imp[imp[,"AcuteChange"]>0 |
imp[,"ChronicChange"]>0,]
row.del =
nrow(imp)-nrow(imp1)
cat(paste(nrow(imp),"
interactions in original input file.\n",sep=""))
cat(paste(nrow(imp1),"
interactions in file for analyses.\n",sep=""))
cat(paste(row.del,"
interactions excluded due to treatment of resilience scores equal to 0 (AC
& CC = 0).\n",sep=""))
# Séparer les colonnes Cotes de risque et Incertitude
des variables de catégories d'interaction.
# Matrice des types de CIE, noms de CIE, activités et
noms des agents de stress:
SEC.A.S = imp1[,interaction.nm]
# Matrice de toutes les cotes de risque et
d'incertitude:
SEC.allRisks =
imp1[,!(names(imp1) %in% interaction.nm)]
# Extraire les cotes de risque et les estimations de
l'incertitude. Compiler dans des matrices distinctes.
# Matrices avec cotes de risque uniquement (une pour
l'exposition, une autre pour la résilience et une autre pour le
rétablissement):
Exposure.score =
risk.subset(exp.nm,SEC.allRisks)
Resilience.score
= risk.subset(res.nm, SEC.allRisks)
Recovery.score =
risk.subset(rec.nm, SEC.allRisks)
# Trame de données des estimations de l'incertitude
relatives aux cotes de risque (sélectionne les colonnes dont les noms
commencent par « I_ »):
Uncertainty =
imp1[grepl("U_", names(imp1))]
# Matrice des cotes d'incertitude échelonnées (1 = 0,2
à 5= 1), comme dans l'étude de Murray et al. (2016) :
U.SD =
Uncertainty*0.2
dimnames(U.SD)[[2]]
= c(exp.nm,res.nm,rec.nm)
# Matrice des cotes d'incertitude uniquement (une pour
l'exposition, une autre pour la résilience et une autre pour le
rétablissement):
U.SD.exp =
risk.subset(exp.nm, U.SD)
U.SD.res =
risk.subset(res.nm, U.SD)
U.SD.rec =
risk.subset(rec.nm, U.SD)
# Quantifier les paramètres de base:
# Listes des CIE et des agents de stress uniques:
SEC = unique(SEC.A.S$SEC)
SEC.A.S$SEC.no =
as.numeric(factor(SEC.A.S$SEC,levels=c(SEC)))
Act =
unique(SEC.A.S$Activity)
SEC.A.S$Activity
= factor(SEC.A.S$Activity, levels = Act, ordered=TRUE)
SEC.A.S$A.S =
paste(SEC.A.S$Sub.Activity,"-", SEC.A.S$Stressor)
SA.S =
unique(SEC.A.S$A.S[order(SEC.A.S$Activity,SEC.A.S$A.S)])
SEC.A.S$A.S.no =
as.numeric(factor(SEC.A.S$A.S, levels=SA.S, ordered=TRUE))
# Nombre de CIE, activités, agents de stress et
composantes de risque uniques:
no.SEC =
length(SEC)
no.Int =
nrow(SEC.A.S)
no.RiskComp =
length(exp.nm)+length(res.nm)+length(rec.nm)
#######################################################
#
EXPOSITION
#######################################################
# Définir les limites supérieure et inférieure:
lwr.bd.exp = 1
upr.bd.exp = 4
# Générer un échantillon aléatoire (taille = nreps) à
l'aide de la distribution adéquate (voir la section PARAMÈTRES DÉFINIS PAR
L'UTILISATEUR):
if (distn == "normal"){
Exposure.rand
= norm.fn(Exposure.score, U.SD.exp, lwr.bd.exp, upr.bd.exp)
}
if (distn == "trunc.norm"){
Exposure.rand
= trnorm.fn(Exposure.score, U.SD.exp, lwr.bd.exp, upr.bd.exp)
}
if (distn == "multinom"){
# Pas
encore mis en œuvre
}
# Attribuer des noms à la deuxième marge:
dimnames(Exposure.rand)[[2]] = exp.nm
# Calculer le chevauchement et l'intensité:
Overlap =
aaply(Exposure.rand[,c("Area","Depth","Temporal"),],1,geometric.mean,na.rm=T)
Intensity =
aaply(Exposure.rand[,c("Intensity","Intensity2"),],1,geometric.mean,na.rm=T)
# Calculer l'exposition:
Exposure = Overlap*Intensity
cat("Exposure calculations complete.\n")
#######################################################
#
CONSÉQUENCE
#######################################################
# RÉSILIENCE
# Définir les limites supérieure et inférieure:
lwr.bd.res = 0
upr.bd.res = 3
# Générer un échantillon aléatoire (taille = nreps) à
l'aide de la distribution adéquate (voir la section PARAMÈTRES DÉFINIS PAR
L'UTILISATEUR):
if (distn == "normal"){
Resilience.rand =
norm.fn(Resilience.score,U.SD.res,lwr.bd.res,upr.bd.res)
}
if (distn == "trunc.norm"){
Resilience.rand =
trnorm.fn(Resilience.score,U.SD.res,lwr.bd.res,upr.bd.res)
}
if (distn == "multinom"){
# Pas
encore mis en œuvre
}
# Attribuer des noms à la deuxième marge:
dimnames(Resilience.rand)[[2]] = res.nm
# Calculer la résilience:
Resilience =
apply(Resilience.rand,c(1,3),sum,na.rm=TRUE)
#
RÉTABLISSEMENT
# Définir les limites supérieure et inférieure:
lwr.bd.rec = 1
upr.bd.rec = 3
# Générer un échantillon aléatoire (taille = nreps) à
l'aide de la distribution adéquate (voir la section PARAMÈTRES DÉFINIS PAR
L'UTILISATEUR):
# Il convient de noter que la commande
suppressWarnings() est utilisée en raison du grand nombre d'avertissements
découlant de tous les
# facteurs de rétablissement = SO.
if (distn == "normal"){
Recovery.rand
= suppressWarnings(norm.fn(Recovery.score,U.SD.rec,lwr.bd.rec,upr.bd.rec))
}
if (distn ==
"trunc.norm"){
Recovery.rand =
suppressWarnings(trnorm.fn(Recovery.score,U.SD.rec,lwr.bd.rec,upr.bd.rec))
}
if (distn == "multinom"){
# Pas
encore mis en œuvre
}
# Attribuer des noms à la deuxième marge:
dimnames(Recovery.rand)[[2]]=
rec.nm
# Calculer le
rétablissement:
Recovery =
apply(Recovery.rand,c(1,3),mean,na.rm=TRUE)
# Calculer les conséquences:
Consequence = Resilience * Recovery
cat("Consequence calculations complete.\n")
###############################################
#
Calculer et résumer le risque
###############################################
# Statistiques récapitulatives relatives à
l'exposition et aux conséquences de chaque interaction:
Exp.Con.Stat = matrix(NA,no.Int,6)
for(i in
1:no.Int) {
Exp.Con.Stat[i,1] =
mean(Exposure[i,],na.rm=T)
Exp.Con.Stat[i,2] =
quantile(Exposure[i,],probs=c(0.1),names=F,na.rm=T)
Exp.Con.Stat[i,3] =
quantile(Exposure[i,],probs=c(0.9),names=F,na.rm=T)
Exp.Con.Stat[i,4] =
mean(Consequence[i,],na.rm=T)
Exp.Con.Stat[i,5] =
quantile(Consequence[i,],probs=c(0.1),names=F,na.rm=T)
Exp.Con.Stat[i,6] = quantile(Consequence[i,],probs=c(0.9),names=F,na.rm=T)
}
colnames(Exp.Con.Stat) = c("Mean
Exposure","ME.10","ME.90","Mean
Consequence","MC.10","MC.90")
# Tableau récapitulatif des cotes attribuées à
l'exposition et aux conséquences de chaque interaction entre CIE et agents de
stress:
Exp.Con.Stat.table = cbind(SEC.A.S,Exp.Con.Stat)
cat("Exposure
and consequence statistics calculated.\n")
# Calculer le risque posé par chaque agent de stress
(Risk.sc) pour chaque CIE pour l'ensemble des échantillons aléatoires
(n=nreps):
Risk.sc = matrix(NA,no.Int,nreps)
Risk.sc = Exposure * Consequence
# Calculer les statistiques récapitulatives de risque
pour chaque combinaison CIE et activité/agent de stress
# (Risque MÉDIAN, quantiles de 10 à 90 %):
Risk.sc.Stats = matrix(NA,no.Int,3)
for(i in
1:no.Int) {
Risk.sc.Stats[i,1] =
median(Risk.sc[i,],na.rm=T)
Risk.sc.Stats[i,2] =
quantile(Risk.sc[i,],probs=c(0.1),names=F,na.rm=T)
Risk.sc.Stats[i,3] =
quantile(Risk.sc[i,],probs=c(0.9),names=F,na.rm=T)
}
colnames(Risk.sc.Stats)
= c("Median Risk","MR.10","MR.90")
# Tableau récapitulatif des cotes de risque médian :
Risk.sc.Stats.table = cbind(SEC.A.S,Risk.sc.Stats)
cat("Risk
statistics calculated.\n")
# Créer un tableau récapitulatif tiré (cf. doc. de recherche,
tableau 3.8) en le limitant aux six premiers
# agents de stress pour chaque CIE :
t3.8 = cbind(SEC.A.S,Risk.sc.Stats,"Mean
Exposure" = Exp.Con.Stat[,"Mean Exposure"],
"ME.10" = Exp.Con.Stat[,"ME.10"],"ME.90" =
Exp.Con.Stat[,"ME.90"],
"Mean Consequence" =
Exp.Con.Stat[,"Mean Consequence"],
"MC.10" =
Exp.Con.Stat[,"MC.10"],"MC.90" =
Exp.Con.Stat[,"MC.90"])
t3.8.order =
t3.8[order(t3.8$SEC.no,-t3.8$"Median Risk"),]
t3.8.order.dt =
as.data.table(t3.8.order)
t3.8.top6 =
t3.8.order.dt[,head(.SD,6), by=.(SEC)]
###################################
# RISQUE
CUMULATIF
###################################
CRisk = rowsum(Risk.sc,SEC.A.S$SEC)
strcount =
rowSums(table(SEC.A.S$SEC,SEC.A.S$Stressor))
# Statistiques récapitulatives pour les cotes de
risque cumulatif:
CRisk.Stats =
array(NA,dim=c(nrow(CRisk),4))
for(i in
1:nrow(CRisk)) {
CRisk.Stats[i,1] =
mean(CRisk[SEC[i],])
CRisk.Stats[i,2] =
quantile(CRisk[SEC[i],],probs=c(0.1),names=F)
CRisk.Stats[i,3] =
quantile(CRisk[SEC[i],],probs=c(0.9),names=F)
CRisk.Stats[i,4] = strcount[SEC[i]]
}
rownames(CRisk.Stats)
= SEC
colnames(CRisk.Stats)
= c("Mean CRisk","10%Q","90%Q","Stressor
Count")
# Tableau récapitulatif du risque cumulatif pour les
CIE découlant de tous les agents de stress connexes :
CRisk.order =
CRisk.Stats[order(CRisk.Stats[,"Mean CRisk"],decreasing=TRUE),]
cat("Cumulative
risk calculated.\n")
###################################
# PUISSANCE
###################################
Potency =
rowsum(Risk.sc,SEC.A.S$A.S)
sppcount =
colSums(table(SEC.A.S$SEC,SEC.A.S$A.S))
# Statistiques récapitulatives des cotes de puissance
:
Potency.Stats =
array(NA,dim=c(nrow(Potency),4))
for(i in
1:nrow(Potency)) {
Potency.Stats[i,1] =
mean(Potency[SA.S[i],])
Potency.Stats[i,2] =
quantile(Potency[SA.S[i],],probs=c(0.1),names=F)
Potency.Stats[i,3] =
quantile(Potency[SA.S[i],],probs=c(0.9),names=F)
Potency.Stats[i,4] = sppcount[SA.S[i]]
}
rownames(Potency.Stats)
= SA.S
colnames(Potency.Stats)
= c("Mean Potency","10%Q","90%Q","Species
Count")
# Tableau récapitulatif de la puissance des agents de
stress pour les CIE connexes:
Potency.order = Potency.Stats[order(Potency.Stats[,"Mean
Potency"],decreasing=TRUE),]
cat("Potency calculated.\n")
###########################################################
# ÉCRIRE LES RÉSULTATS
###########################################################
setwd(d)
write.table(Exp.Con.Stat.table,file="ExpCon.Stats.csv",col.names=NA,sep=",")
write.table(Risk.sc.Stats.table,file="MedianRisk.Stats.csv",col.names=NA,sep=",")
write.table(t3.8.order,file=paste("SKB_FullOrderedRiskResults_",distn,".csv",sep="")
,col.names=NA,sep=",")
write.table(t3.8.top6,file=paste("SKB_Top6RiskResults_",distn,".csv",sep=""),
col.names=NA,sep=",")
write.table(CRisk.order,file="CumRiskSEC.csv",col.names=NA,sep=",")
write.table(Potency.order,file="PotencySEC.csv",col.names=NA,sep=",")
###################
# FIGURES
###################
if (do.figures == "Y"){
# La mise en forme des étiquettes des CIE (doivent
correspondre à l'ordre des objets CIE qui se trouvent ligne 162):
SEC.format = c("Bocaccio
Rockfish","Corralline Algae","Corals","Pacific
Halibut",
expression(paste(italic("Isidella"),sep="")),"Macroalgae",
expression(paste(italic("Primnoa"),sep="")),"Prowfish",
"Rougheye Rockfish",
"Sablefish","Sponges","Squat
Lobster","Widow Rockfish", "Yelloweye Rockfish")
####################################
# FIGURE : cote de risque médian et diagrammes de
double projection de l'exposition et des conséquences pour chaque CIE.
####################################
# Réaménager l'ordre des CIE afin qu'il concorde avec
les figures du document de recherche (en fonction des
# objets CIE créés à la ligne 162):
fig.SEC.order = c(1,6,9,10,7,3,5,8,4,2,12,14,11,13)
secs.pp = c(4,3,3,4) # Nombre de CIE par page
biplot.dat = cbind(Risk.sc.Stats.table,Exp.Con.Stat)
for (page in
1:length(secs.pp)){
jpeg(filename=paste("Fig",page,"_ExpConRisk.jpg",sep=""),width=8.5,height=11,units="in",res=300)
par(mar=c(2.5,4,1,1),oma=c(0.5,1,1,1))
no.figs = secs.pp[page]*2
layout(matrix(c(1:no.figs,no.figs+1,no.figs+1),secs.pp[page]+1,2,byrow=TRUE),heights=c(rep.int(1,secs.pp[page]),0.9))
pt.colour =
c("black","red","dodgerblue","grey")
for (j in 1:secs.pp[page]){
t.sec =
SEC[fig.SEC.order[j+sum(secs.pp[0:(page-1)])]]
mrisk = biplot.dat[biplot.dat$SEC==t.sec,]
# Diagramme des cotes de risque médian
assorties de barres d'erreur:
plot(mrisk$A.S.no, mrisk$"Median
Risk", ylab="", xlab='',
xlim=c(1,max(biplot.dat$A.S.no)),
xaxt='n', pch=21, col="black",
ylim=c(0,(max(biplot.dat[,"MR.90"]))), cex=1.5,
bg =
pt.colour[as.integer(mrisk$Activity)])
arrows(mrisk$A.S.no,
mrisk$"MR.10", mrisk$A.S.no,mrisk$"MR.90",
length=0,col="grey70")
points(mrisk$A.S.no, mrisk$"Median
Risk", pch=21, col="black",
bg =
pt.colour[as.integer(mrisk$Activity)], cex=2)
mtext(SEC.format[fig.SEC.order[j+sum(secs.pp[0:(page-1)])]], side=3, adj=1)
mtext("Median Risk", side=2,
line=2.5, cex=0.9)
if(j==secs.pp[page]){mtext("Stressors",side=1,line=2,cex=0.9)}
axis(side=1,
at=seq(1,max(biplot.dat$A.S.no), by=2),
labels=seq(1,max(biplot.dat$A.S.no),by=2))
#abline(v=c(6.5,10.5),lty=2)
if(j==1){legend("topright",c("Vessel
Traffic","Seismic
Surveys","Research","Fishing"),col="black",fill=pt.colour,cex=0.8,inset=0.05)}
# Diagramme de double projection de
l'exposition et des conséquences:
jitter =
rnorm(1,0.001,0.001) #modifiera
l'emplacement des étiquettes de manière aléatoire pour éviter les
chevauchements
mrisk.order
= mrisk[order(-mrisk$"Median Risk"),] #avec la matrice nécessaire pour
schématiser les quatre CIE principales dont le risque médian est le plus élevé.
plot(mrisk$"Mean
Exposure",mrisk$"Mean
Consequence",xlab="",ylab="",pch=1,xlim=c(0,18),ylim=c(0,16),axes=FALSE)
arrows(mrisk$"ME.10",mrisk$"Mean
Consequence",mrisk$"ME.90",
mrisk$"Mean Consequence",length=0,col="lightgrey")
arrows(mrisk$"Mean
Exposure",mrisk$"MC.10",mrisk$"Mean Exposure",
mrisk$"MC.90",length=0,col="lightgrey")
points(mrisk$"Mean
Exposure",mrisk$"Mean Consequence",pch=21,col="black",
bg = pt.colour[as.integer(mrisk$Activity)])
text((mrisk.order$"Mean
Exposure"[1:4]+0.4+jitter),(mrisk.order$"Mean
Consequence"[1:4]+0.4+jitter),
mrisk.order$A.S.no[1:4], cex=0.5)
mtext(SEC.format[fig.SEC.order[j+sum(secs.pp[0:(page-1)])]], side=3,
adj=1)
mtext("Consequence", side=2,
line=2, cex=0.9)
axis(side=1, at=seq(0,18,2),
labels=seq(0,18,2), cex=0.7)
axis(side=2, at=seq(0,16,2),
labels=seq(0,16,2), las=1, cex=0.7)
box()
if(j==secs.pp[page]){mtext("Exposure", side=1, line=2,
cex=0.9)}
}
n.col = 3
S.cols =
ceiling(length(SA.S)/n.col)
plot(1:20,1:20,xlim=c(0,15),ylim=c(0,5),type="n",axes=FALSE,xlab="",ylab="")
if(length(SA.S)>9){leg.pch
= c(paste(1:9,"
",SA.S[1:9],sep=""),paste(10:length(SA.S)," ",SA.S[10:length(SA.S)],sep=""))}
if(length(SA.S)<=9){leg.pch
= c(paste(1:length(SA.S),"
",SA.S,sep=""))}
legend("center",SA.S,pch="",legend=leg.pch,fill="white",col="black",border="white",cex=0.7,ncol=3)
dev.off()
}
cat("Median
Risk/Exposure-Consequence figures complete.\n")
###################################
# FIGURE : diagramme des cotes de risque cumulatif des
CIE assorties de barres d'erreur.
###################################
numsp =
c(1:length(SEC))
cumdat =
CRisk.order
cumsp =
SEC.format[order(CRisk.Stats[,"Mean CRisk"],decreasing=TRUE)]
jpeg(filename=paste("Fig",length(secs.pp)+1,"_CumulativeRisk.jpg",sep=""),width=8.5,height=7,units="in",res=300)
par(mar=c(12,6,2,2))
barcum<-barplot(cumdat[,"Mean
CRisk"],names.arg=cumsp,las=2,ylim=c(0,max(cumdat[,3])+50),cex.names=0.9)
arrows(barcum,(cumdat[,2]),barcum,(cumdat[,3]),length=0)
text(barcum+0.25,cumdat[,1]+25,cumdat[,4])
title(ylab="Cumulative
Risk",line=2.5)
box()
dev.off()
cat("Cumulative
risk figure complete.\n")
###################################
# FIGURE 3.5: diagramme des cotes de puissance des
agents de stress assortis de barres d'erreur.
###################################
numstr =
c(1:length(SA.S))
potdat = Potency.order
potstr = rownames(Potency.order)
# Dresser la liste des activités par ordre décroissant
en fonction des résultats de puissance :
pot.Act.l =
unlist(gregexpr(pattern = " - ",potstr))
pot.Act =
paste(" (",substring(potstr,1,pot.Act.l-1),")
",sep="")
pot.Str =
paste(substring(potstr,pot.Act.l+3,nchar(potstr)),"
",sep="")
#Créer une
figure:
jpeg(filename=paste("Fig",length(secs.pp)+2,"_Potency.jpg",sep=""),width=8.5,height=7,units="in",res=300)
par(mar=c(15,6,2,2))
barpot<-barplot(potdat[,"Mean
Potency"],ylim=c(0,max(potdat[,3])+50),axisnames=FALSE,cex.names=0.8)
arrows(barpot,(potdat[,2]),barpot,(potdat[,3]),length=0)
text(barpot+0.35,potdat[,1]+25,potdat[,4],cex=0.8)
mtext(pot.Str,side=1,
at=barpot-0.25,las=2,cex=0.8)
mtext(pot.Act,side=1,
at=barpot+0.25,las=2, col= "grey30",cex=0.6)
title(ylab="Stressor
Potency",line=2.5)
box()
dev.off()
cat("Potency
figure complete.\n")
}
setwd("../")