-
Notifications
You must be signed in to change notification settings - Fork 0
/
pgrcalc.R
175 lines (153 loc) · 6.27 KB
/
pgrcalc.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
170
171
172
173
174
175
# Calculate population parameters for a range of start
# times with varying forward years. Saves output for
# Total, 1SW and MSW. Separate file for max forward
# year.
#
# Filename: pgrcalc.R
# Author: Darren Kavanagh
# Created: 30.07.2009
# Last Modified: 06.08.2008
# ======================================================
#
# load libraries and functions
source('calcp.R')
L <- 5 # set L value
mints = 15 # minimum number of time series counts
stryr <- 1971 # start year
endyr <- 2008 # end year
maxstyr <- (endyr - mints) # max start year
strtrng <- c(stryr:maxstyr) # start year range
notss <- length(strtrng) # no of years in range
# ----------------------------- #
# Load & Run Analysis #
# ----------------------------- #
# spawner data
cdata <- read.csv("./data/nea_spawners.csv", header = TRUE, sep = ",", quote="\"", dec=".", fill = TRUE, comment.char="")
cdata$dataset<-factor(cdata$dataset, levels=unique(cdata$dataset))
# cdata <- subset(cdata, dataset == "NE Atlantic") # limit for testing
# PFA data
cpfa <- read.csv("./data/nea_pfa.csv", header = TRUE, sep = ",", quote="\"", dec=".", fill = TRUE, comment.char="")
# Returns data
crtn <- read.csv("./data/nea_returns.csv", header = TRUE, sep = ",", quote="\"", dec=".", fill = TRUE, comment.char="")
ordds <- unique(cdata$dataset)
orgrp <- unique(cdata$dsgroup)
nods <- length(ordds) # number of datasets
# find number of records to allocate
nrec <- 0 # rec counter
for (j in 1:nods) {
for (k in 1:notss) {
nofy <- (maxstyr+1) - strtrng[k]
for (i in 1:nofy) {
nrec <- nrec + 1
}
}
}
nomx <- notss*nods
# output data frames
# total
totgr <- data.frame(dataset=factor(character(nrec), levels=(ordds)),dsgroup=factor(character(nrec),levels=orgrp), syear=numeric(nrec), fyear=numeric(nrec), ctoratio=numeric(nrec), ctmratio=numeric(nrec), apsize=numeric(nrec), tmrate=numeric(nrec), dmrate=numeric(nrec), lmrate=numeric(nrec), L=numeric(nrec), mu=numeric(nrec), sigmasq=numeric(nrec), lambda=numeric(nrec), lambdal95=numeric(nrec), lambdah95=numeric(nrec))
totmxgr <- totgr[1:nomx,]
# 1SW
oswgr <- totgr
oswmxgr <- totmxgr
# MSW
mswgr <- totgr
mswmxgr <- totmxgr
cnt <- 0; mxcnt <- 0 # counters
print(paste("No of datasets:", nods))
# run analysis for each dataset
for (j in 1:nods) {
curds <- as.character(ordds[j])
cdssub <- subset(cdata, dataset == curds) # spawners
cpfasub <- subset(cpfa, dataset == curds) # PFA
crtnsub <- subset(crtn, dataset == curds) # returns
print(paste("Dataset",j,curds,"processing"))
# vary start year
for (k in 1:notss) {
mxcnt <- mxcnt + 1
csyear <- strtrng[k]
stsds <- subset(cdssub, year >= csyear) # spawners
stpfa <- subset(cpfasub, year >= csyear) # PFA
strtn <- subset(crtnsub, year >= csyear) # returns
# vary no years into the future
nofy <- (maxstyr+1) - csyear
for (i in 1:nofy) {
cnt <- cnt + 1
cfyear <- mints+i-1
ftsds <- subset(stsds, year <= (csyear+cfyear)) # spawners
ftspfa <- subset(stpfa, year <= (csyear+cfyear)) # PFA
ftsrtn <- subset(strtn, year <= (csyear+cfyear)) # returns
ftslnt <- length(ftsds$dataset)
cds <- as.character(ftsds$dataset[1])
cdsgrp <- as.character(ftsds$dsgroup[1])
# 1SW/MSW ratios
ctoratio <- mean(ftsds$osw/ftsds$total)
ctmratio <- mean(ftsds$msw/ftsds$total)
curinp <- data.frame(cds, cdsgrp, csyear, cfyear, ctoratio, ctmratio)
# total
ctotap <- mean(ftsds$total) # total average pop size
# total mortality rate
tctmrate <- mean((ftspfa$totpfa-ftsds$total)/ftspfa$totpfa)
# distant mortality rate
tcdmrate <- mean((ftspfa$totpfa-ftsrtn$totret)/ftspfa$totpfa)
# local mortality rate
tclmrate <- mean((ftsrtn$totret-ftsds$total)/ftspfa$totpfa)
# calc growth rate parameters
totp <- calcp(ftsds$total, L)
# populate total data frame
totgr[cnt,] <- data.frame(curinp, ctotap, tctmrate, tcdmrate, tclmrate, L, totp)
# 1SW
coswap <- mean(ftsds$osw) # osw average pop size
# total mortality rate
octmrate <- mean((ftspfa$oswpfa-ftsds$osw)/ftspfa$oswpfa)
# distant mortality rate
ocdmrate <- mean((ftspfa$oswpfa-ftsrtn$oswret)/ftspfa$oswpfa)
# local mortality rate
oclmrate <- mean((ftsrtn$oswret-ftsds$osw)/ftspfa$oswpfa)
# calc growth rate parameters
oswp <- calcp(ftsds$osw, L)
# populate total data frame
oswgr[cnt,] <- data.frame(curinp, coswap, octmrate, ocdmrate, oclmrate, L, oswp)
# MSW
cmswap <- mean(ftsds$msw) # msw average pop size
# total mortality rate
mctmrate <- mean((ftspfa$mswpfa-ftsds$msw)/ftspfa$mswpfa)
# distant mortality rate
mcdmrate <- mean((ftspfa$mswpfa-ftsrtn$mswret)/ftspfa$mswpfa)
# local mortality rate
mclmrate <- mean((ftsrtn$mswret-ftsds$msw)/ftspfa$mswpfa)
# calc growth rate parameters
mswp <- calcp(ftsds$msw, L)
# populate total data frame
mswgr[cnt,] <- data.frame(curinp, cmswap, mctmrate, mcdmrate, mclmrate, L, mswp)
} # end of future years loop
# allocate max forward years
totmxgr[mxcnt,] <- totgr[cnt,] # total
oswmxgr[mxcnt,] <- oswgr[cnt,] # 1SW
mswmxgr[mxcnt,] <- mswgr[cnt,] # MSW
} # end of start year loop
} # end of dataset loop
# ----------------------------- #
# Generate Output #
# ----------------------------- #
# total
fileoutnm <- paste("./data/", "totgr.csv", sep = "")
write.table(totgr, file = fileoutnm, sep = ",", col.names = NA,
qmethod = "double")
fileoutnm <- paste("./data/", "totmxgr.csv", sep = "")
write.table(totmxgr, file = fileoutnm, sep = ",", col.names = NA,
qmethod = "double")
# 1SW
fileoutnm <- paste("./data/", "oswgr.csv", sep = "")
write.table(oswgr, file = fileoutnm, sep = ",", col.names = NA,
qmethod = "double")
fileoutnm <- paste("./data/", "oswmxgr.csv", sep = "")
write.table(oswmxgr, file = fileoutnm, sep = ",", col.names = NA,
qmethod = "double")
# MSW
fileoutnm <- paste("./data/", "mswgr.csv", sep = "")
write.table(mswgr, file = fileoutnm, sep = ",", col.names = NA,
qmethod = "double")
fileoutnm <- paste("./data/", "mswmxgr.csv", sep = "")
write.table(mswmxgr, file = fileoutnm, sep = ",", col.names = NA,
qmethod = "double")