-
Notifications
You must be signed in to change notification settings - Fork 1
/
do.pssd.troph.ag.r
169 lines (133 loc) · 6.83 KB
/
do.pssd.troph.ag.r
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
# Function for generating NOEC distributions for each species taking into account trophic levels
# --- EXAMPLE OF A SPECIAL CASE FOR SILVER
#
# Yields a matrix with dimensions number of species * number of iterations in the simulation
#
# Arguments: - DP : matrix of data points
# - UFt : matrix of uncertainty factors for the exposure time
# - UFdd : matrix of uncertainty factors for the dose-descriptor
# - SIM : number of iterations in the simulation
# - CV.DP : coefficient of variation for the interlaboratory variation
# - CV.UF : coefficient of variation for the use of non-substance-specific
# uncertainty factors
# - TrophLevel : Trophic levels of each species. 1 for primary producers, 2 for herbivores and 3 for carnivores
# - Troph.Perc.Typ : Typical fractions in which each trophic level is represented in freshwater communities.
#
# Date of last modification: 03.07.2019
#
# Associated publication: Systematic consideration of parameter uncertainty and variability in
# probabilistic species sensitivity distributions
#
# Authors: Henning Wigger, Delphine Kawecki, Bernd Nowack and Veronique Adam
#
# Institute: Empa, Swiss Federal Laboratories for Materials Science and Technology,
# Technology and Society Laboratory, Lerchenfeldstrasse 5, 9014 St. Gallen, Switzerland
#
# submitted to Integrated Environmental Assessment and Management in July 2019
# -------------------------------------------------------------------------------------------------
do.pSSD.troph.Ag <- function(DP,
UFt,
UFdd,
SIM,
CV.DP,
CV.UF,
TrophLevel,
Troph.Perc.Typ = c(0.65, 0.25, 0.1)){
# test if there is no data available for one species
if(any(apply(DP,2,function(x) length(which(!is.na(x)))) == 0)){
warning("No data is available for one or more species, it/they won't contribute to the PSSD calculation.")
# find which species has no data
ind.sp.rem <- which(apply(DP,2,function(x) length(which(!is.na(x)))) == 0)
# remove those columns
DP <- DP[,-ind.sp.rem]
UFt <- UFt[,-ind.sp.rem]
UFdd <- UFdd[,-ind.sp.rem]
}
# Create the step distributions (or triangular or trapezoidal) for each species
# Create an empty matrix in which step distributions will be compiled
NOEC_comb <- matrix(NA, ncol(DP), SIM,
dimnames = list(colnames(DP), NULL))
# Fill in the matrix. If there is only one data point, NOEC stays the same. If there are
# 2 endpoints, a uniform distribution is produced. If there are more than 2 endpoints, a step
# distribution is produced. One line is for one species.
require(trapezoid)
require(mc2d)
# store the corrected endpoints
corr.endpoints <- DP/(UFdd*UFt)
sort.endpoints <- apply(corr.endpoints,2, sort)
for (sp in colnames(DP)){
# store the indices of the minimal and maximal data point
ind.min <- which.min(corr.endpoints[,sp])
ind.max <- which.max(corr.endpoints[,sp])
# calculate the theoretical minimum and maximum of the distribution we are looking for
if (UFt[ind.min, sp] == 10){
sp.min <- corr.endpoints[ind.min,sp]*(1-(sqrt((CV.DP/2.45)^2)))/100
} else {
sp.min <- corr.endpoints[ind.min,sp]*(1-(sqrt((CV.DP/2.45)^2 + (CV.UF/2.45)^2 + (CV.UF/2.45)^2)*2.45))
}
if (UFt[ind.max, sp] == 10){
sp.max <- corr.endpoints[ind.max,sp]*(1+(sqrt((CV.DP/2.45)^2)))/1
} else {
sp.max <- corr.endpoints[ind.max,sp]*(1+(sqrt((CV.DP/2.45)^2 + (CV.UF/2.45)^2 + (CV.UF/2.45)^2)*2.45))
}
# For species with one unique data point, NOEC stays the same:
if(length(unique(sort.endpoints[[sp]])) == 1){
NOEC_comb[sp,] <- rtrunc("rtriang", min = sp.min,
mode = sort.endpoints[[sp]][1],
max = sp.max,
n = SIM, linf = 0)
# For species with two endpoints:
} else if(length(sort.endpoints[[sp]]) == 2){
# Create a trapezoidal distribution including both endpoints
NOEC_comb[sp,] <- rtrunc("rtrapezoid", SIM,
mode1 = sort.endpoints[[sp]][1],
mode2 = sort.endpoints[[sp]][2],
min = sp.min, max = sp.max,
linf = 0)
# For species with three endpoints or more:
} else {
# Sample from this step distribution for each species
NOEC_comb[sp,] <- rmore(values = sort.endpoints[[sp]], max = sp.max, min = sp.min, N = SIM, linf = 0)
}
}
# Sample species with typical proportions representative of a typical community (if terrestrial, Troph.Perc need to be changed)
Nb_Troph <- c(length(which(TrophLevel == 1)),
length(which(TrophLevel == 2)),
length(which(TrophLevel == 3)))
# calculate the shares of the different trophic levels in the data
Troph.Perc.Data <- Nb_Troph/ncol(DP)
# identify the trophic level that will be unchanged
if(length(which(Troph.Perc.Data < Troph.Perc.Typ)) == 1){
lev.fix <- which(Troph.Perc.Data < Troph.Perc.Typ)
} else {
lev.fix <- which.max((Troph.Perc.Typ - Troph.Perc.Data)/Troph.Perc.Typ)
}
# create an empty vector to contain the corrected number of species for the species being corrected
Nb_points <- rep(NA,3)
# create an empty list with an element per trophic level
NOEC_lev <- vector("list", 3)
for(lev in 1:3){
# for fix level, take as such
if(lev == lev.fix){
Nb_points[lev] <- Nb_Troph[lev.fix]
# identify species corresponding to level lev
ind.sp <- which(TrophLevel == lev) # vector of all indices for which the species level is the fixed one
# store in list
NOEC_lev[[lev]] <- NOEC_comb[ind.sp,]
} else {
# change the remaining levels
Nb_points[lev] <- round(Nb_Troph[lev.fix] * Troph.Perc.Typ[lev] / Troph.Perc.Data[lev.fix]) # Nouveau nombre de lignes
ind.sp <- which(TrophLevel == lev) # vector of all indices of the level
# store in temporary matrix the NOEC distributions of all species of the level
temp <- NOEC_comb[ind.sp,]
NOEC_lev[[lev]] <- matrix(NA, Nb_points[lev],SIM)
for(i in 1:SIM){
NOEC_lev[[lev]][,i] <- sample(temp[,i], Nb_points[lev])
}
}
}
# combine that stuff
NOEC_TrophCorr <- do.call(rbind,NOEC_lev)
# return the whole matrix
return(NOEC_TrophCorr)
}