-
Notifications
You must be signed in to change notification settings - Fork 35
/
mcmc.out.R
420 lines (392 loc) · 14.7 KB
/
mcmc.out.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
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
#' Summarize, analyze and plot key MCMC output.
#'
#' Makes four panel plot showing trace plots, moving average, autocorrelations,
#' and densities for chosen parameters from MCMC output.
#'
#'
#' @param directory Directory where all results are located, one level above
#' directory for particular run.
#' @param run Directory with files from a particular run.
#' @template file
#' @templateVar file_t posteriors
#' @param namefile The (optional) file name of the dimension and names of
#' posteriors.
#' @param names Read in names file (T) or use generic naming (F).
#' @param headernames Use the names in the header of `file`?
#' @param numparams The number of parameters to analyze.
#' @param closeall By default close all open devices.
#' @param burn Optional burn-in value to apply on top of the option in the
#' starter file and [SSgetMCMC()].
#' @param thin Optional thinning value to apply on top of the option in the
#' starter file, in the `-mcsave` runtime command, and in
#' [SSgetMCMC()].
#' @param scatter Can add a scatter-plot of all params at end, default is none.
#' @param surface Add a surface plot of 2-way correlations.
#' @param surf1 The first parameter for the surface plot.
#' @param surf2 The second parameter for the surface plot.
#' @param stats Print stats if desired.
#' @param plots Show plots or not.
#' @param header Data file with header?
#' @param sep Separator for data file passed to the `read.table` function.
#' @param print Send to screen unless asked to print.
#' @param new Logical whether or not to open a new plot window before plotting
#' @param colNames Specific names of the `file` to extract and work with. `NULL`
#' keeps all columns
#' @author Ian Stewart, Allan Hicks (modifications)
#' @return `directory`, because this function is used for its plotting side effects
#' @export
#' @seealso [mcmc.nuisance()], [SSgetMCMC()]
#' @examples
#' \dontrun{
#' mcmc.df <- SSgetMCMC(
#' dir = "mcmcRun", writecsv = T,
#' keystrings = c("NatM", "R0", "steep", "Q_extraSD"),
#' nuisancestrings = c("Objective_function", "SSB_", "InitAge", "RecrDev")
#' )
#' mcmc.out("mcmcRun", run = "", numparams = 4, closeall = F)
#'
#' # Or for more control
#' par(mar = c(5, 3.5, 0, 0.5), oma = c(0, 2.5, 0.2, 0))
#' mcmc.out("mcmcRun",
#' run = "",
#' numparams = 1,
#' closeall = F,
#' new = F,
#' colNames = c("NatM_p_1_Fem_GP_1")
#' )
#' mtext("M (natural mortality)", side = 2, outer = T, line = 1.5, cex = 1.1)
#' }
mcmc.out <- function(directory = "c:/mydirectory/",
run = "mymodel/",
file = "keyposteriors.csv",
namefile = "postplotnames.sso",
names = FALSE,
headernames = TRUE,
numparams = 1,
closeall = TRUE,
burn = 0,
thin = 1,
scatter = FALSE,
surface = FALSE,
surf1 = 1,
surf2 = 2,
stats = FALSE,
plots = TRUE,
header = TRUE,
sep = ",",
print = FALSE,
new = T,
colNames = NULL)
##############################################################################################################
{
# add section to set up for printing or display to screen (default)
if (print == TRUE) {} # not implemented
if (closeall == TRUE) { # see if the user asked to retain open graphics devices
# useful to compare multiple runs
### Note: the following line has been commented out because it was identified
### by Brian Ripley as "against CRAN policies".
# rm(.SavedPlots,pos=1) # remove any plotting history
}
filename <- file.path(directory, run, file) # put directory,run and file names together for use
# warning if file does not exist
if (!file.exists(filename)) {
stop("file doesn't exist:\n", filename)
}
mcmcdata <- read.table(filename, # make data table of whole file
header = header, # choice of header or not
sep = sep, # space delimited
fill = TRUE
) # fill empty cells to make a symmetrical array
#### Naming section ####
if (names == TRUE) {
nameout <- file.path(directory, run, namefile) # put directory,run and file names together for use
namedata <- read.table(nameout, # make data table of whole file
header = FALSE, # no headers
sep = "", # space delimited
colClasses = "character", # don't convert to factors
fill = TRUE
) # fill empty cells to make a symmetrical array
numparams <- as.numeric(namedata[1, 1]) # get the dimension of the output
# add names to the data table, only used in the scatterplot of all parameters
for (j in 1:numparams) { # loop over the moveparam columns
names(mcmcdata)[j] <- namedata[(j + 1), 1] # name each column
}
}
if (!is.null(colNames)) {
if (length(colNames) != numparams) {
warning("numparams argument overidden by length of colNames argument")
}
numparams <- length(colNames)
mcmcdata <- mcmcdata[, colNames]
if (length(colNames) == 1) {
mcmcdata <- data.frame(mcmcdata)
names(mcmcdata) <- colNames
}
}
##### change to mcmc object for coda #####
mcmcfirst <- mcmc(mcmcdata) # make the mcmc object from the data table
mcmctemp <- window(mcmcfirst, thin = thin, start = (1 + burn)) # thin the chain and remove burn in
mcthinned <- as.matrix(mcmctemp) # get rid of iteration labels
mcmcobject <- mcmc(mcthinned) # send back to mcmc object
draws <- length(mcmcobject[, 1]) # define the post thinning and burn in length of the chain
##### plotting section #####
if (plots == TRUE) { # have plots been activated by user
if (new) dev.new(record = TRUE) # keep the window open for each parameter
if (numparams == 5 || numparams == 9 || numparams == 13 || numparams == 17) { # plots a blank plot if 5,9,13, or 17 plots created
# this avoids the loss of plot n-1 in history (is this an R bug??)
plot(0, 0,
xlab = "",
ylab = "",
frame.plot = FALSE,
yaxt = "n", # turn off this axis
xaxt = "n",
type = "n"
) # plot nothing
}
for (i in 1:numparams) { # loop over the number of parameters
par(
new = FALSE, # make sure to use a new graphical window
mfrow = c(2, 2), # set up "cells" to graph into
ann = TRUE
) # annotate plots
##### Trace plot section #####
traceplot(mcmcobject[, i], # trace plot of parameters
smooth = TRUE
) # add a smoothing line
mtext("Value", # label for y-axis
side = 2, # place it on left of the graph
line = 3, # set the distance above the graph
font = 1, # make the font regular
cex = 0.8
) # scale the text size
if (names | headernames) {
mtext(names(mcmcdata)[i], # label for whole plotting page
side = 3, # place it on left of the graph
adj = 0, # left adjust the text
line = 2, # set the distance above the graph
font = 2, # make the font regular
cex = 1
) # scale the text size
} else {
mtext(paste("param", i), # label for whole plotting page
side = 3, # place it on left of the graph
adj = 0, # left adjust the text
line = 2, # set the distance above the graph
font = 2, # make the font regular
cex = 1
) # scale the text size
}
#### Cumulative plot section ####
lowest <- min(mcmcobject[, i]) # minimum value in chain
highest <- max(mcmcobject[, i]) # maximum value in chain, for graphing
plot(c(seq(1, draws, by = 1)), # the x values
c(lowest, rep(c(highest), (draws - 1))), # the y values
xlab = "Iterations",
ylab = "",
yaxt = "n", # turn off this axis
type = "n"
) # plot nothing
if (!requireNamespace("gtools", quietly = TRUE)) {
warning("Package \"gtools\" needed for the running average plot. Please install it.",
call. = FALSE
)
} else {
lines(gtools::running(mcmcobject[, i],
fun = median, # plot the mean
allow.fewer = TRUE, # begin calculating from the first point
width = draws
)) # Averages the last - observations (all in this case)
fun <- function(x, prob) quantile(x, probs = prob, names = FALSE) # the quantile to use in next 2 lines
lines(
gtools::running(mcmcobject[, i],
fun = fun, # function to use
prob = 0.05,
allow.fewer = TRUE, # begin calculating from the first point
width = draws
),
col = "GREY"
)
lines(
gtools::running(mcmcobject[, i],
fun = fun, # function to use
prob = 0.95,
allow.fewer = TRUE, # begin calculating from the first point
width = draws
),
col = "GREY"
)
} # end temporary turning off
#### Autocorrelation plot section ####
par(ann = FALSE) # turn off default annotation
autocorr.plot(mcmcobject[, i],
auto.layout = FALSE, # do not make a new graph sheet
lag.max = 20, # set the maximum lag
ask = FALSE
) # do not require prompt to move through plots
mtext("Autocorrelation", # label for y-axis
side = 2, # place it on left of the graph
line = 3, # set the distance above the graph
font = 1, # make the font regular
cex = 0.8
) # scale the text size
mtext("Lag", # label for y-axis
side = 1, # place it on left of the graph
line = 3, # set the distance above the graph
font = 1, # make the font regular
cex = 0.8
) # scale the text size
lines(seq(1, 20, by = 1), rep(0.1, 20), col = "GREY") # plot the 0.1 line
lines(seq(1, 20, by = 1), rep(-0.1, 20), col = "GREY") # plot the -0.1 line
##### Density plot section #####
densplot(mcmcobject[, i],
show.obs = TRUE
) # show the x axis
mtext("Density", # label for y-axis
side = 2, # place it on left of the graph
line = 3, # set the distance above the graph
font = 1, # make the font regular
cex = 0.8
) # scale the text size
mtext("Value", # label for y-axis
side = 1, # place it on left of the graph
line = 3, # set the distance above the graph
font = 1, # make the font regular
cex = 0.8
) # scale the text size
} # end parameter loop
} # end plotting loop
#### Statistics section #####
if (stats == TRUE) {
dev.new() # opens a new graphics device
par(mar = c(0, 0, 3, 0)) # makes the margins large
plot(0, # plot a graph of a single point
ylab = "", # label the y-axis for whole screen
xlab = "", # label the x-axis for whole screen
type = "n", # plot nothing in the graph
xlim = c(0, 25),
ylim = c(0, 25),
main = "Summary statistics for key parameters", # label the output
axes = FALSE
) # do not plot the axis
# set up the titles
text(0.001, # the x-axis location
25, # the y-axis location
"Parameter", # what to write
font = 2, # bold font
cex = 0.9, # font size
adj = 0
)
text(4, # the x-axis location
25, # the y-axis location
"Median (0.05-0.95)", # what to write
font = 2, # bold font
cex = 0.9, # font size
adj = 0
)
text(13, # the x-axis location
25, # the y-axis location
"AC Lag 1", # what to write
font = 2, # bold font
cex = 0.9, # font size
adj = 0
)
text(16.5, # the x-axis location
25, # the y-axis location
"Eff. N", # what to write
font = 2, # bold font
cex = 0.9, # font size
adj = 0
)
text(19, # the x-axis location
25, # the y-axis location
"Geweke-Z", # what to write
font = 2, # bold font
cex = 0.9, # font size
adj = 0
)
text(22.5, # the x-axis location
25, # the y-axis location
"Heidel-W", # what to write
font = 2, # bold font
cex = 0.9, # font size
adj = 0
)
# loop over parameters and fill in the values
for (i in 1:numparams)
{
text(0, # the x-axis location
(25 - i), # the y-axis location
paste("param", i), # what to write
font = 1, # normal font
cex = 0.9, # font size
adj = 0
) # make the text left-aligned
med <- quantile(mcmcobject[, i], probs = 0.5, names = FALSE)
range <- quantile(mcmcobject[, i], probs = c(0.05, 0.95), names = FALSE)
text(3.2, 25 - i,
paste(
signif(round(med, 6), 6),
"(",
paste(
signif(round(range[1], 6), 6),
"-",
signif(round(range[2], 6), 6)
),
")"
),
font = 1, cex = 0.9, adj = 0
)
l1.ac <- acf(mcmcobject[, i], lag.max = 1, type = "correlation", plot = F)
acoruse <- round(l1.ac[["acf"]][2], 6)
text(13, 25 - i, acoruse, font = 1, cex = 0.9, adj = 0)
effsize <- effectiveSize(mcmcobject[, i])
text(16.5, 25 - i, round(min(effsize, draws), 0), font = 1, cex = 0.9, adj = 0)
if (acoruse > 0.4) {
gewuse <- "None"
}
if (acoruse <= 0.4) {
geweke <- geweke.diag(mcmcobject[, i], frac1 = 0.1, frac2 = 0.5)
gewuse <- round(geweke[["z"]], 3)
}
text(19, 25 - i, gewuse, font = 1, cex = 0.9, adj = 0)
if (acoruse > 0.4) {
send <- "None"
}
if (acoruse <= 0.4) {
hw <- as.list(heidel.diag(mcmcobject[, i], pvalue = 0.05))
if (hw[1] == 0) {
send <- "Failed"
}
if (hw[1] == 1) {
send <- "Passed"
}
}
text(22.5, 25 - i, send, font = 1, cex = 0.9, adj = 0)
} # close the parameter loop
} # end stats statement
##### Scatter plot section #####
if (scatter == TRUE) {
dev.new()
par(xaxt = "n", yaxt = "n") # suppress the axis labels
pairs(mcmcdata[1:numparams], # make the scatterplot
cex = 0.1,
gap = 0
)
} # end the scatter loop
##### Surface plot section #####
# note: depends on gplots which is no longer installed by default
if (surface == TRUE) {
dev.new()
par(new = FALSE) # use a new window
gplots::hist2d(mcmcobject[, surf1], # x data as a vector
mcmcobject[, surf2], # y data as a vector
nbins = 100,
na.rm = TRUE,
xlab = paste("parameter", surf1),
ylab = paste("parameter", surf2),
show = TRUE, # make figure or not
col = c("GREY", topo.colors(20))
)
} # close surface plot loop'
return(directory)
} # end function