This repository was archived by the owner on Jun 13, 2023. It is now read-only.
-
Notifications
You must be signed in to change notification settings - Fork 5
/
Copy pathComputeIndices.R
145 lines (137 loc) · 9.66 KB
/
ComputeIndices.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
#' Compute index of abundance
#'
#' @param Data The data to be fit
#' @param Model The model object (as a list)
#' @param FileName The name of the file/directory
#' @param maxDims aesthetics argument for plot as dimensions of rows/columns, defaults to 6
#' @param Folder Optional argument for named folder to store results (in current working directory)
#' @param Weights Defaults to "StrataAreas", other option is "Equal"
#' @param StrataTable Table of strata, required
#' @param PlotStrataYearMcmc Whether to plot MCMC chains for each Strata:Year, defaults to \code{TRUE}
#'
#' @return Returns a list of output, by year and by Strata:Year
#' @import R2jags
#' @export
#'
ComputeIndices <- function(Data, Model, FileName, maxDims = 6, Folder = NA,
Weights = "StrataAreas", StrataTable, PlotStrataYearMcmc = TRUE) {
if (is.na(Folder)) Folder <- paste(getwd(), "/", sep = "")
# Attach stuff -- listed by search()
attach(Model$BUGSoutput$sims.list, warn.conflicts = FALSE)
on.exit(detach(Model$BUGSoutput$sims.list))
modelStructure <- Model$modelStructure
Dist <- Model$likelihood
# Compute average AreaSwept - This AreaSwept used to calculate densities, and matters if Catch is a nonlinear function of AreaSwept
MeanLogAreaSwept <- mean(Data$logeffort)
# New objects
Chains <- array(NA, dim = c(nrow(Sdev), length(unique(StrataTable[, "year"])), length(unique(StrataTable[, "strata"])), 2))
Year <- Strata <- Area <- PosMedian <- PresMedian <- IndexMedian <- IndexMedianWeighted <- PosMean <- PresMean <- IndexMean <- IndexMeanWeighted <- CvMedian <- SdLog <- RawPos <- RawPres <- Raw <- RawWeighted <- RawVar <- RawVarWeighted <- RawCV <- matrix(NA, nrow = length(unique(StrataTable[, "year"])), ncol = length(unique(StrataTable[, "strata"])))
CvMedianYear <- SdLogYear <- rep(NA, length(unique(StrataTable[, "year"])))
# Calculate values
for (YearI in 1:nrow(PosMean)) {
for (StratI in 1:ncol(PosMean)) {
# Derived indicators
StrataYearI <- which(levels(strataYear) == paste(levels(strata)[StratI], ":", levels(year)[YearI], sep = ""))
AreaI <- which(StrataTable[, "strataYear"] == paste(levels(strata)[StratI], ":", levels(year)[YearI], sep = ""))
Which <- which(strata == levels(strata)[StratI] & year == levels(year)[YearI])
# Year, strata, and area
Year[YearI, StratI] <- levels(year)[YearI]
Strata[YearI, StratI] <- levels(strata)[StratI]
if (Weights == "StrataAreas") Area[YearI, StratI] <- StrataTable[AreaI, "Area_Hectares"]
if (Weights == "Equal") Area[YearI, StratI] <- 1
# Save Positive catch chain in normal-space and correct for transformation biases
Chains[, YearI, StratI, 1] <- exp(cMx(Sdev)[, StratI] + cMx(Ydev)[, YearI] + cMx(SYdev)[, StrataYearI] + log(1) * B.pos[, 1] + log(1)^2 * B.pos[, 2]) # wardJAGS uses logeffort offset
# Lognormal -- Bias correction
if (Dist == "lognormal") {
Sigma <- sqrt(log(CV[, 1]^2 + 1))
Chains[, YearI, StratI, 1] <- Chains[, YearI, StratI, 1] * exp(Sigma^2 / 2)
}
# LognormalECE -- Bias correction + incorporate ECE
if (Dist == "lognormalECE" | Dist == "lognormalECE2") {
Sigma <- sqrt(log(CV[, 1]^2 + 1))
Sigma2 <- sqrt(log(CV[, 2]^2 + 1))
Chains[, YearI, StratI, 1] <- Chains[, YearI, StratI, 1] * exp(Sigma^2 / 2) * p.ece[, 1] + ratio * Chains[, YearI, StratI, 1] * exp(Sigma2^2 / 2) * p.ece[, 2]
}
# GammaECE -- incorporate ECE
if (Dist == "gammaECE" | Dist == "gammaECE2") {
Chains[, YearI, StratI, 1] <- Chains[, YearI, StratI, 1] * p.ece[, 1] + ratio * Chains[, YearI, StratI, 1] * p.ece[, 2]
}
# Don't make mean-unbiased for unobserved vessel
# if(modelStructure$VesselYear.positiveTows=="random") Chains[,YearI,StratI,1] = Chains[,YearI,StratI,1] * exp(sigmaVY[,1]^2/2)
# Save Presence/Absence chain in normal-space and correct for transformation biases
Chains[, YearI, StratI, 2] <- plogis(cMx(pSdev)[, StratI] + cMx(pYdev)[, YearI] + cMx(pSYdev)[, StrataYearI] + (1) * B.zero[, 1] + (1)^2 * B.zero[, 2]) # wardJAGS predicts the probability of 0 catch, and uses offset as effort, not logeffort
# Don't make mean-unbiased for unobserved vessel (this was done incorrectly anyway)
# if(modelStructure$VesselYear.zeroTows=="random") Chains[,YearI,StratI,2] = mean(plogis(rnorm(n=dim(Chains)[1]*1e3, mean=qlogis(Chains[,YearI,StratI,2]), sd=sigmaVY[,2])))
# Index (median)
PosMedian[YearI, StratI] <- median(Chains[, YearI, StratI, 1]) # / 2e4 # Convert kilograms to metric tons
PresMedian[YearI, StratI] <- median(Chains[, YearI, StratI, 2])
IndexMedian[YearI, StratI] <- median(Chains[, YearI, StratI, 1] * Chains[, YearI, StratI, 2])
IndexMedianWeighted[YearI, StratI] <- IndexMedian[YearI, StratI] * Area[YearI, StratI]
# Index (mean)
PosMean[YearI, StratI] <- mean(Chains[, YearI, StratI, 1]) # / 2e4 # Convert kilograms to metric tons
PresMean[YearI, StratI] <- mean(Chains[, YearI, StratI, 2])
IndexMean[YearI, StratI] <- mean(Chains[, YearI, StratI, 1] * Chains[, YearI, StratI, 2])
IndexMeanWeighted[YearI, StratI] <- IndexMean[YearI, StratI] * Area[YearI, StratI]
# CV of median (from J. Wallace "Survey.Biomass.GlmmBUGS.ver.3.00.r)
Temp <- Area[YearI, StratI] * Chains[, YearI, StratI, 1] * Chains[, YearI, StratI, 2]
CvMedian[YearI, StratI] <- sqrt(var(Temp)) / median(Temp)
# SD of log of index (from J. Wallace "Survey.Biomass.GlmmBUGS.ver.3.00.r")
SdLog[YearI, StratI] <- sd(log(Temp))
# Raw
RawPos[YearI, StratI] <- mean(ifelse(Data[Which, "HAUL_WT_KG"] > 0, Data[Which, "HAUL_WT_KG"] / Data[Which, "effort"], NA), na.rm = TRUE)
RawPres[YearI, StratI] <- mean(Data[Which, "HAUL_WT_KG"] > 0)
Raw[YearI, StratI] <- RawPos[YearI, StratI] * RawPres[YearI, StratI]
RawWeighted[YearI, StratI] <- Raw[YearI, StratI] * Area[YearI, StratI]
RawVar[YearI, StratI] <- var(Data[Which, "HAUL_WT_KG"] / Data[Which, "effort"], na.rm = TRUE) / length(Which)
RawVarWeighted[YearI, StratI] <- RawVar[YearI, StratI] * Area[YearI, StratI]^2
RawCV[YearI, StratI] <- sqrt(RawVarWeighted[YearI, StratI]) / RawWeighted[YearI, StratI]
}
# CV of median (from J. Wallace "Survey.Biomass.GlmmBUGS.ver.3.00.r")
# Temp = Area[YearI,StratI] * rowSums( cMx(Chains[, YearI, ,1]) * cMx(Chains[,YearI,,2]) ) # By using StratI outside of the Strata loop only the final area estimate is obained.
Temp <- (cMx(Chains[, YearI, , 1]) * cMx(Chains[, YearI, , 2])) %*% array(Area[YearI, ]) # array(Area[YearI,], c(3,1))
CvMedianYear[YearI] <- sqrt(var(Temp)) / median(Temp)
# SD of log of index (from J. Wallace "Survey.Biomass.GlmmBUGS.ver.3.00.r")
SdLogYear[YearI] <- sd(log(Temp))
} # 1085-115
# Plot MCMC chains for each Strata:Year
maxDims2 <- maxDims^2
if (PlotStrataYearMcmc == TRUE) {
for (Type in c("Trace", "ACF")) {
for (ChainI in 1:2) {
Nstrat <- length(unique(strataYear[nonZeros]))
Nplots <- ceiling(Nstrat / maxDims2)
for (PlotI in 1:Nplots) {
ParSet <- (PlotI - 1) * maxDims^2 + 1:min(Nstrat - (PlotI - 1) * maxDims2, maxDims2)
Ncol <- ceiling(sqrt(length(ParSet)))
Nrow <- ceiling(length(ParSet) / Ncol)
jpeg(paste(Folder, "/", FileName, "Chain_", c("Positive", "Presence")[ChainI], "_", Type, "_by_StrataYear_", PlotI, ".jpg", sep = ""), width = Ncol * 1.5, height = Nrow * 1.5, res = 150, units = "in")
par(mfrow = c(Nrow, Ncol), mar = c(2, 2, 2, 0), mgp = c(1.25, 0.25, 0), tck = -0.02)
for (StrataYearI in ParSet) {
StrataI <- strata[which(strataYear == unique(strataYear)[StrataYearI])[1]]
YearI <- year[which(strataYear == unique(strataYear)[StrataYearI])[1]]
Mat <- matrix(Chains[, YearI, StrataI, ChainI], ncol = Model$mcmc.control$chains, byrow = FALSE)
if (Type == "Trace") {
matplot(Mat, type = "l", lty = "solid", col = rainbow(Model$mcmc.control$chains, alpha = 0.4), main = paste(YearI, StrataI), xaxt = "n", xlab = "", ylab = "")
}
if (Type == "ACF") {
Acf <- apply(Mat, MARGIN = 2, FUN = function(Vec) {
acf(Vec, plot = FALSE)$acf
})
matplot(Acf, type = "h", lty = "solid", col = rainbow(Model$mcmc.control$chains, alpha = 0.4), main = paste(YearI, StrataI), xaxt = "n", xlab = "", ylab = "", ylim = c(0, 1), lwd = 2)
}
}
dev.off()
}
}
}
}
# Compile into matrices
Results1 <- data.frame(Year = as.vector(Year), Strata = as.vector(Strata), Raw = as.vector(RawWeighted), RawCV = as.vector(RawCV), IndexMedian = as.vector(IndexMedianWeighted), IndexMean = as.vector(IndexMeanWeighted), CvMedian = as.vector(CvMedian), SdLog = as.vector(SdLog), Area = as.vector(Area), PosMedian = as.vector(PosMedian), PresMedian = as.vector(PresMedian), PosMean = as.vector(PosMean), PresMean = as.vector(PresMean), RawPos = as.vector(RawPos), RawPres = as.vector(RawPres), RawSD = ifelse(as.vector(RawVarWeighted) == 0, NA, as.vector(sqrt(RawVarWeighted))))
Results2 <- data.frame(Year = Year[, 1], Raw = rowSums(RawWeighted, na.rm = TRUE), RawCV = sqrt(rowSums(RawVarWeighted, na.rm = TRUE)) / rowSums(RawWeighted, na.rm = TRUE), IndexMedian = rowSums(IndexMedianWeighted, na.rm = TRUE), IndexMean = rowSums(IndexMeanWeighted, na.rm = TRUE), CvMedian = CvMedianYear, SdLog = SdLogYear)
# Write and print output
write.csv(Results1, file = paste(Folder, "/", FileName, "ResultsByYearAndStrata.csv", sep = ""))
write.csv(Results2, file = paste(Folder, "/", FileName, "ResultsByYear.csv", sep = ""))
# Return output
Return <- list(byYearAndStrata = Results1, byYear = Results2)
return(Return)
}