-
Notifications
You must be signed in to change notification settings - Fork 0
/
addSSsummarize.R
300 lines (280 loc) · 17.1 KB
/
addSSsummarize.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
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
#' Add a model to the list of models to compare
#'
#' Adds specified quantities from any model to the list of models returned from
#' \code{\link{SSsummarize}} for further comparison.
#'
#' @param origModels A list of models created by \code{\link{SSsummarize}}.
#' @param newModels A list of models to add to the originals models list. Each
#' new model is an element of the list, and is a list itself with possible
#' components described in the details below.
#' @return Returns list as is returned from \code{\link{SSsummarize}}, but
#' contains additions for the new models.
#' @note This function was made to compare TINSS results and SS results, and
#' assumed that you would always start with a list of SS models output from
#' \code{\link{SSsummarize}}. It has not been tested to see how it works when
#' starting with an empty list.
#' @author Allan Hicks
#' @export
#' @seealso \code{\link{SSsummarize}} \code{\link{SSplotComparisons}}
#' @examples
#'
#' \dontrun{
#' ######################################
#' #DO NOT RUN
#' tinss1 <- list(npars=A$fit$npar,maxgrad=A$fit$maxgrad,nsexes=1,
#' #note, there is an estimated parameter called sd_sbt,
#' # but it is a single value
#' SpawnBio=data.frame(c(1964,1965,A$yrs),
#' c(A$sbo,A$sbo,A$sbt)*1e6,0,
#' qnorm(0.025,c(A$so,A$so,A$sbt)*1e6,0),
#' qnorm(0.975,c(A$so,A$so,A$sbt)*1e6,0)),
#' Bratio=data.frame(A$yrs,A$sbt/A$sbo,0,
#' qnorm(0.025,A$sbt/A$sbo,0),
#' qnorm(0.975,A$sbt/A$sbo,0)),
#' SPRratio=data.frame(A$yr,A$spr,0,qnorm(0.025,A$spr,0),
#' qnorm(0.975,A$spr,0)),
#' recruits=data.frame(A$yr,A$nt[,1]*1e6,0,qnorm(0.025,A$nt[,1]*1e5,0),
#' qnorm(0.975,A$nt[,1]*1e6,0)),
#' #I'm not sure exactly what wt are,
#' # but it is important to line them up correctly
#' recdevs=data.frame(A$recYrs,A$wt),
#' indices = data.frame(A$iyr,1e6*A$yt,1e6*A$qbt,
#' rep(A$q,length(A$iyr)),rep(0.4,length(A$iyr)),
#' rep(0,length(A$iyr)),rep(1,length(A$iyr)))
#' )
#' tinss <- list(tinss1,tinss1) #can add more models here
#'
#'
#' #add TINSS model to SS models already summarized
#' SSnTINSS <- addSSsummarize(models,tinss)
#' mcmcInd <- seq(burnin+1,nrow(A$mc.sbt),thin)
#' SSnTINSS$mcmc[[2]] <- data.frame(A$mc.sb0[mcmcInd],
#' A$mc.sbt[mcmcInd,],
#' A$mc.depl[mcmcInd,],
#' A$mc.spr[mcmcInd,],
#' A$mc.rt[mcmcInd,],
#' log(A$mcmc[mcmcInd,"Ro"]*1e6),
#' A$mcmc[mcmcInd,"msy"]*1e6)
#' names(SSnTINSS$mcmc[[2]]) <-
#' c("SPB_Virgin",paste("SPB",A$yrs,sep="_"),
#' paste("Bratio",A$yrs,sep="_"),
#' paste("SPRratio",A$yr,sep="_"),
#' paste("Recr",A$yr,sep="_"),"SR_R0","TotYield_MSY")
#' modelnames <- c("SS", "TINSS","TINSS.MLE")
#'
#' SSplotComparisons(SSnTINSS, legendlabels=modelnames,
#' subplot=2,endyr=2011,mcmcVec=c(T,T,F))
#' title(main="MCMC")
#' SSplotComparisons(SSnTINSS, legendlabels=modelnames,
#' subplot=4,endyr=2011,mcmcVec=c(T,T,F))
#' title(main="MCMC")
#' ###############################################
#' }
#'
addSSsummarize <- function(origModels,newModels) {
########################################################################
# This function adds the specified components to a model summary object
# If there are no original models, it will create an object that can be used with r4ss model comparisons
#
# Written by Allan Hicks, 1/12/2011, [email protected]
#
# origModels is the list of original models that are being added to
# it is the result of the SSsumarize function
# if origModels is NULL, then a new list will be returned that can be used in r4ss model comparisons
#
# newModels is a list of models to be added to the origModels
# each new model is an element of the list, and is a list itself
# each new model is a list with possible components named:
# npars: the number of parameters in the model
# maxgrad: the maximum gradient component (if used)
# nsexes: the number of sexes
# likelihoods: likelihoods from the model. A data.frame with the 2nd column as names, which matches on names from origModels
# SS uses the following names
# TOTAL
# Equil_catch
# Survey
# Length_comp
# Age_comp
# Recruitment
# recast_Recruitment
# Parm_priors
# Parm_softbounds
# Parm_devs
# Crash_Pen
# Size_at_age
# likelambdas NOT IMPLEMENTED
# #pars#: NOT IMPLEMENTED YET FOR DIFFICULTY IN MATCHING PARAMETERS
# #parsSD#: NOT IMPLEMENTED YET FOR DIFFICULTY IN MATCHING PARAMETERS
# #parsphases#: NOT IMPLEMENTED YET FOR DIFFICULTY IN MATCHING PARAMETERS
# SpawnBio: Spawning biomass matrix
# 1st column is year
# 2nd column is spawning biomass in same units as original models (SS reports female spawning biomass)
# 3rd column is the standard deviation of estimated spawning biomass
# 4th column is a lower bound of the confidence interval to be plotted (say from an MCMC)
# 5th column is an upper bound of the confidence interval to be plotted (say from an MCMC)
# Bratio: Depletion matrix
# 1st column is year
# 2nd column is depletion
# 3rd column is the standard deviation of depletion (optional)
# 4th column is a lower bound of the confidence interval to be plotted (say from an MCMC)
# 5th column is an upper bound of the confidence interval to be plotted (say from an MCMC)
# SPRratio: SPR ratio matrix
# 1st column is year
# 2nd column is depletion
# 3rd column is the standard deviation (optional)
# 4th column is a lower bound of the confidence interval to be plotted (say from an MCMC)
# 5th column is an upper bound of the confidence interval to be plotted (say from an MCMC)
# recruits: Recruitment matrix
# 1st column is year
# 2nd column is recruitment as in original models (SS reports age-0 recruits)
# 3rd column is the standard deviation (optional)
# 4th column is a lower bound of the confidence interval to be plotted (say from an MCMC)
# 5th column is an upper bound of the confidence interval to be plotted (say from an MCMC)
# recdevs: Recruitment deviate matrix
# 1st column is year
# 2nd column is deviate (matched with original models)
# growth: #NOT IMPLEMENTED
# indices: Matrix of fits to indices
# 1st column is year
# 2nd column is observed index (data)
# 3rd column is expected index (prediction)
# 4th column is catchability coefficient (q)
# 5th column is standard error of index (total used in fitting)
# 6th column is a likelihood for this point, or enter any value to make sure it plots, or enter NA not to plot the estimate
#
# InitAgeYrs #NOT IMPLEMENTED
if(!is.list(newModels[[1]])) stop("The new models needs to be a list of lists. A list of new models, each with a list of quantities to enter.")
cat("Adding",length(newModels),"new models to",origModels$n,"original models\n")
if(is.null(origModels)) {
stop("Not implemented yet for no original models\n")
#origModels <- list()
}
addColumn <- function(origDF,newDF,matchOrig,matchNew,nOrig,valNew=1) {
#adds a column for the new model to the dataframe and matches using the original column labeled with matchOrig and the new
#nOrig is the number of original models
if(is.null(newDF)) {
newDF <- data.frame(rep(NA,nrow(origDF)),origDF[,matchOrig])
valNew <- 1
matchNew <- 2
}
x.add <- newDF[!(newDF[,matchNew] %in% origDF[,matchOrig]),]
if(nrow(x.add)>0) { #add additional rows for missing names
origDFnrow <- nrow(origDF)
for(i in 1:nrow(x.add)) {origDF <- rbind(origDF,NA)} #add neessary number of rows of NA
origDF[,matchOrig] <- c(origDF[1:origDFnrow,matchOrig],as.character(x.add[,matchNew]))
}
x.ind <- match(origDF[,matchOrig],newDF[,matchNew]) #should have all rows of newDF and origDF
out <- cbind(origDF[,1:nOrig],newDF[x.ind,valNew],origDF[,(nOrig+1):ncol(origDF)])
colnames(out) <- c(1:(nOrig+1),colnames(origDF)[(nOrig+1):ncol(origDF)])
out
}
addNewModel.fn <- function(x,models) {
#adds a single model to the origModels
#models is the set of models formatted as in SSsummarize
#written as a function to use lapply, avoid loops and make it easier when origModels=NULL
n <- models$n
models$n <- n+1
if(!is.null(x$npars)) { models$npars[n+1] <- x$npars }else{ models$npars[n+1] <- NA }
if(!is.null(x$listnames)) { models$npars[n+1] <- x$listnames }else{ models$npars[n+1] <- NA }
if(!is.null(x$keyvec)) { models$npars[n+1] <- x$keyvec }else{ models$npars[n+1] <- NA }
if(!is.null(x$maxgrad)) { models$maxgrad[n+1] <- x$maxgrad }else{ models$maxgrad[n+1] <- NA }
if(!is.null(x$nsexes)) { models$nsexes[n+1] <- x$nsexes }else{ cat("nsexes must be defined! Inserting a value of 3.\n"); models$nsexes[n+1] <- 3}
models$pars <- addColumn(models$pars,NULL,"Label",2,n) ###NEED TO IMPLEMENT
models$parsSD <- addColumn(models$parsSD,NULL,"Label",2,n) ###NEED TO IMPLEMENT
models$parphases <- addColumn(models$parphases,NULL,"Label",2,n) ###NEED TO IMPLEMENT
models$quants <- addColumn(models$quants,NULL,"Label",2,n) ###NEED TO IMPLEMENT
models$quantsSD <- addColumn(models$quantsSD,NULL,"Label",2,n) ###NEED TO IMPLEMENT
models$likelihoods <- addColumn(models$likelihoods,x$likelihoods,"Label",2,n)
models$likelambdas <- addColumn(models$likelambdas,NULL,"Label",2,n) ###NEED TO IMPLEMENT
#add Spawning Biomasses, note that it is only the column that changes
if(!is.null(x$SpawnBio)) print("adding SpawnBio")
if(!is.null(x$SpawnBio)){if(ncol(x$SpawnBio)<5) {stop("There must be at least 5 columns in SpawnBio: Yr, SB, SD, Lower, Upper\n")}}
models$SpawnBio <- addColumn(models$SpawnBio,x$SpawnBio,"Yr",1,n,2)
models$SpawnBio$Yr <- as.numeric(models$SpawnBio$Yr)
models$SpawnBio <- models$SpawnBio[order(models$SpawnBio[,"Yr"]),]
models$SpawnBioSD <- addColumn(models$SpawnBioSD,x$SpawnBio,"Yr",1,n,3)
models$SpawnBioSD$Yr <- as.numeric(models$SpawnBioSD$Yr)
models$SpawnBioSD <- models$SpawnBioSD[order(models$SpawnBioSD[,"Yr"]),]
models$SpawnBioLower <- addColumn(models$SpawnBioLower,x$SpawnBio,"Yr",1,n,4)
models$SpawnBioLower$Yr <- as.numeric(models$SpawnBioLower$Yr)
models$SpawnBioLower <- models$SpawnBioLower[order(models$SpawnBioLower[,"Yr"]),]
models$SpawnBioUpper <- addColumn(models$SpawnBioUpper,x$SpawnBio,"Yr",1,n,5)
models$SpawnBioUpper$Yr <- as.numeric(models$SpawnBioUpper$Yr)
models$SpawnBioUpper <- models$SpawnBioUpper[order(models$SpawnBioUpper[,"Yr"]),]
#add Bratio (depletion), note that it is only the column that changes
if(!is.null(x$Bratio)) {
if(ncol(x$Bratio)<5) {stop("There must be at least 5 columns in Bratio: Yr, SB, SD, Lower, Upper\n")}
print("adding Bratio")
}
models$Bratio <- addColumn(models$Bratio,x$Bratio,"Yr",1,n,2)
models$Bratio$Yr <- as.numeric(models$Bratio$Yr)
models$Bratio <- models$Bratio[order(models$Bratio[,"Yr"]),]
models$BratioSD <- addColumn(models$BratioSD,x$Bratio,"Yr",1,n,3)
models$BratioSD$Yr <- as.numeric(models$BratioSD$Yr)
models$BratioSD <- models$BratioSD[order(models$BratioSD[,"Yr"]),]
models$BratioLower <- addColumn(models$BratioLower,x$Bratio,"Yr",1,n,4)
models$BratioLower$Yr <- as.numeric(models$BratioLower$Yr)
models$BratioLower <- models$BratioLower[order(models$BratioLower[,"Yr"]),]
models$BratioUpper <- addColumn(models$BratioUpper,x$Bratio,"Yr",1,n,5)
models$BratioUpper$Yr <- as.numeric(models$BratioUpper$Yr)
models$BratioUpper <- models$BratioUpper[order(models$BratioUpper[,"Yr"]),]
#add recruits
models$recruits <- addColumn(models$recruits,x$recruits,"Yr",1,n,2)
models$recruitsSD <- addColumn(models$recruitsSD,x$recruits,"Yr",1,n,3)
models$recruitsLower <- addColumn(models$recruitsLower,x$recruits,"Yr",1,n,4)
models$recruitsUpper <- addColumn(models$recruitsUpper,x$recruits,"Yr",1,n,5)
models$recdevs <- addColumn(models$recdevs,x$recdevs,"Yr",1,n,2)
#add SPR
models$SPRratio <- addColumn(models$SPRratio,x$SPRratio,"Yr",1,n,2)
models$SPRratioSD <- addColumn(models$SPRratioSD,x$SPRratio,"Yr",1,n,3)
models$SPRratioLower <- addColumn(models$SPRratioLower,x$SPRratio,"Yr",1,n,4)
models$SPRratioUpper <- addColumn(models$SPRratioUpper,x$SPRratio,"Yr",1,n,5)
models$growth <- cbind(models$growth,NA) #NEED TO IMPLEMENT
#Indices
if(!is.null(x$indices)) {
newIndices <- as.data.frame(cbind(NA,x$indices[,1],NA,NA,x$indices[,2],x$indices[,3],x$indices[,4],NA,x$indices[,5],NA,x$indices[,6],NA,NA,NA,NA,n+1,n+1))
names(newIndices) <- names(models$indices)
models$indices <- rbind(models$indices,newIndices)
}
#models$InitAgeYrs <- cbind(models$InitAgeYrs,NA) #NEED TO IMPLEMENT
return(models)
}
for(ii in 1:length(newModels)) { #lapply doesn't work because it returns a separate list for each newModel
origModels <- addNewModel.fn(newModels[[ii]],models=origModels)
}
return(origModels)
}
if(F) {
likes <- c(4242,5555,9933,3377)
likes <- data.frame(likes=likes,label=c("Survey","Other","TOTAL","Other2"))
SB <- data.frame(c(1971,1972,1973),c(3333,5555,7777),c(100,200,300),c(1111,3333,5555),c(5555,7777,9999))
BR <- data.frame(c(1971,1972,1973),c(0.3333,0.5555,0.7777),c(0.100,.200,.300),c(.1111,.3333,.5555),c(.5555,.7777,.9999))
newModel1 <- list(npars=42,likelihoods=likes,SpawnBio=SB,Bratio=BR) #data.frame(1,"NatM_p_1_Fem_GP_1"))
likes <- c(4233,5544,9988,3344)
likes <- data.frame(likes=likes,label=c("Recruitment","Other","TOTAL","Other2"))
newModel2 <- list(npars=42,likelihoods=likes)
newModels <- list(newModel1,newModel2)
tmp <- addSSsummarize(mysummary,newModels)
x.ind <- match(models$likelihoods[,"Label"],names(x$likelihoods))
models$likelihoods <- cbind(models$likelihoods[,1:n],x$likelihoods[x.ind],models$likelihoods[,n+1])
colnames(models$likelihoods) <- c(1:(n+1),"Label")
x.add <- x$likelihoods[!(names(x$likelihoods) %in% models$likelihoods[,"Label"])]
if(length(x.add)>0) { #add additional rows for missing names
tmp <- as.data.frame(matrix(NA,nrow=length(x.add),ncol=ncol(models$likelihoods),dimnames=list(NULL,colnames(models$likelihoods))))
tmp[,n+1] <- x.add
tmp[,"Label"] <- names(x.add)
models$likelihoods <- rbind(models$likelihoods,tmp)
}
addNewModel.fn <- function(x,models) {
#adds a single model to the origModels
#models is the set of models formatted as in SSsummarize
#written as a function to use lapply, avoid loops and make it easier when origModels=NULL
n <- models$n
models$n <- n+1
if(!is.null(x$npars)) { models$npars[n+1] <- x$npars }else{ models$npars[n+1] <- NA }
if(!is.null(x$maxgrad)) { models$maxgrad[n+1] <- x$maxgrad }else{ models$maxgrad[n+1] <- NA }
if(!is.null(x$nsexes)) { models$nsexes[n+1] <- x$nsexes }else{ models$nsexes[n+1] <- NA }
models$likelihoods <- addColumn(models$likelihoods,x$likelihoods,"Label",2,n)
return(models)
}
}