# Nest survival model using JAGS
# Only fixed effects at this time
#
# 2019-06-30 CHJS First Edition
#
library("R2jags") # used for call to JAGS
## Loading required package: rjags
## Loading required package: coda
## Linked to JAGS 4.3.0
## Loaded modules: basemod,bugs
##
## Attaching package: 'R2jags'
## The following object is masked from 'package:coda':
##
## traceplot
library(coda)
library(ggplot2)
## Registered S3 methods overwritten by 'ggplot2':
## method from
## [.quosures rlang
## c.quosures rlang
## print.quosures rlang
library(reshape2)
options(width=200)
source(file.path("..","..","jags-nest-survival-fixed-effects.r"))
# The input dataframe must contain the following fields with the following names
#
# NestID: id code of the nest (alpha numeric)
# FirstFound: day the nest was first found
# LastPresent: last day that a chick was present in the nest
# LastChecked: last day the nest was checked
# Fate: fate of the nest; 0 = success; 1=fail
# AgheDay1 = age of the nest on day 1 (if you are fitting age of nest models)
#
# You could also have a nest level covariates, survey level covariates, and
# next x survey time covariates as well
nestdata <- readxl::read_excel(file.path("..","Killdeer.xlsx"),
sheet="killdeer-age")
nestdata <- plyr::rename(nestdata, c("id"="NestId"))
head(nestdata)
## # A tibble: 6 x 7
## NestId FirstFound LastPresent LastChecked Fate Freq AgeDay1
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 /*A*/ 1 9 9 0 1 0
## 2 /*B*/ 5 5 9 1 1 -2
## 3 /*C*/ 5 40 40 0 1 -3
## 4 /*D*/ 9 32 32 0 1 -4
## 5 /*E*/ 7 8 8 0 1 -4
## 6 /*F*/ 3 15 15 0 1 1
# Unfortunately, JAGS cannot deal with alpha numeric code and
# so we need to convert the alphanumberic NestID to numeric codes
# by declaring NestId as a factor and extracting the level values
nestdata$NestId.num <- as.numeric(factor(nestdata$NestId))
# We must create a file with every combination of next x day nest was "active"
# being every day from FirstCound to LastChecked-1
nesttime <- plyr::adply(nestdata, 1, function(x){
nesttime <- expand.grid(NestId.num=x$NestId.num,
Day=x$FirstFound:(x$LastChecked-1),
stringsAsFactors=FALSE)
nesttime
})
# Extract the nest level covariates (including AgeNest1)
# The next level covariates should be indexed using NestId
# If AgeNest1 variable is present then the age of the nest is computed
#
nest.covariates <- NULL
if( !is.null(nest.covariates)){
nesttime <- merge(nesttime, nest.covariates, by="NestId")
}
# Extract any survey time covariates such as time, time^2, early/late
# weather covariates ect.
# All of these covariates will affect all nests simultaneouls
nesttime $Day2 <- (nesttime$Day-20)^2 # day^2 for quadratic trends
nesttime $Period <- car::recode(nesttime$Day,
paste("lo:", (max(nesttime$Day)+min(nesttime$Day))/2, "='Early';",
"else='Late'"))
xtabs(~Period+Day, data=nesttime, exclude=NULL, na.action=na.pass)
## Day
## Period 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
## Early 1 1 2 2 4 4 5 6 5 5 5 5 6 6 6 8 7 7 7 7 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## Late 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 8 8 9 9 9 9 10 10 9 9 9 6 4 3 3 3 2 2 2
# if there is a AgeDay1 variable, we compute the nest age for each time for each nest
if( !is.null(nesttime$AgeDay1)){
nesttime$NestAge <- nesttime$AgeDay1 + nesttime$Day -1
}
head(nesttime)
## NestId FirstFound LastPresent LastChecked Fate Freq AgeDay1 NestId.num Day Day2 Period NestAge
## 1 /*A*/ 1 9 9 0 1 0 1 1 361 Early 0
## 2 /*A*/ 1 9 9 0 1 0 1 2 324 Early 1
## 3 /*A*/ 1 9 9 0 1 0 1 3 289 Early 2
## 4 /*A*/ 1 9 9 0 1 0 1 4 256 Early 3
## 5 /*A*/ 1 9 9 0 1 0 1 5 225 Early 4
## 6 /*A*/ 1 9 9 0 1 0 1 6 196 Early 5
# Add any next x day survey covariates to the nesttime data
#
# there is nothing here for this example
# Set up the design matrix for the fixed effects
fe.design <- model.matrix( ~ 1, data=nesttime)
head(fe.design)
## (Intercept)
## 1 1
## 2 1
## 3 1
## 4 1
## 5 1
## 6 1
# Finally, load and run the JAGS model
fitted.model <- jags.nest.survival.fixed.effects(
nestdata=nestdata, # nest data
nesttime=nesttime, # daily nest values with nest, time, nest x time covariates
fe.design=fe.design, # fixed effects design matrix
init.seed=12321312) # initial seed)
## module glm loaded
## Compiling data graph
## Declaring variables
## Resolving undeclared variables
## Allocating nodes
## Initializing
## Reading data back into data table
## Compiling model graph
## Declaring variables
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 22
## Unobserved stochastic nodes: 1
## Total graph size: 803
##
## Initializing model
# the nesttime dataframe has the estimated DSR for every combination of NestId.num and Day
# in this case, we fit a S~1 model, so all rows of nesttime have the same estimated DSR
fitted.model$nesttime[1,]
## NestId.num Day NestId FirstFound LastPresent LastChecked Fate Freq AgeDay1 Day2 Period NestAge mean sd X2.5. X25. X50. X75. X97.5. Rhat n.eff
## 1 1 1 /*A*/ 1 9 9 0 1 0 361 Early 0 0.9722 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
# the results list has lots of other stuff
results <- fitted.model$results
names(results)
## [1] "model" "BUGSoutput" "parameters.to.save" "model.file" "n.iter" "DIC"
names(results$BUGSoutput)
## [1] "n.chains" "n.iter" "n.burnin" "n.thin" "n.keep" "n.sims" "sims.array" "sims.list" "sims.matrix" "summary"
## [11] "mean" "sd" "median" "root.short" "long.short" "dimension.short" "indexes.short" "last.values" "program" "model.file"
## [21] "isDIC" "DICbyR" "pD" "DIC"
# we can also look at the beta estimates
# in this case this is the logit DSR which is the same for all nest x days
results$BUGSoutput$summary[ grepl("beta", row.names(results$BUGSoutput$summary)),,drop=FALSE]
## mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff
## beta 3.636509 0.4282816 2.884788 3.334836 3.601548 3.90071 4.571069 1.00093 4500
#######################################
# get the full summary table
results$BUGSoutput$summary
## mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff
## S[1,1] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[1,2] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[1,3] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[6,3] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[1,4] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[6,4] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[1,5] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[2,5] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[3,5] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[6,5] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[1,6] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[2,6] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[3,6] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[6,6] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[1,7] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[2,7] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[3,7] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[5,7] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[6,7] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[1,8] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[2,8] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[3,8] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[6,8] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[7,8] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[9,8] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[3,9] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[4,9] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[6,9] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[7,9] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[9,9] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[3,10] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[4,10] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[6,10] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[7,10] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[9,10] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[3,11] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[4,11] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[6,11] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[7,11] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[9,11] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[3,12] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[4,12] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[6,12] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[7,12] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[9,12] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[3,13] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[4,13] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[6,13] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[7,13] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[9,13] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[10,13] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[3,14] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[4,14] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[6,14] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[7,14] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[8,14] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[11,14] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[3,15] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[4,15] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[7,15] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[8,15] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[11,15] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[12,15] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[3,16] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[4,16] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[7,16] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[11,16] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[12,16] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[13,16] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[14,16] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[15,16] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[3,17] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[4,17] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[7,17] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[11,17] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[12,17] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[13,17] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[14,17] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[3,18] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[4,18] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[7,18] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[11,18] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[12,18] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[13,18] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[14,18] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[3,19] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[4,19] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[7,19] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[11,19] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[12,19] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[13,19] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[14,19] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[3,20] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[4,20] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[7,20] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[11,20] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[12,20] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[13,20] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[14,20] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[3,21] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[4,21] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[7,21] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[11,21] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[12,21] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[13,21] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[14,21] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[16,21] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[3,22] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[4,22] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[7,22] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[11,22] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[12,22] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[13,22] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[14,22] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[16,22] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[3,23] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[4,23] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[7,23] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[11,23] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[12,23] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[13,23] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[14,23] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[16,23] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[17,23] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[3,24] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[4,24] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[7,24] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[11,24] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[12,24] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[13,24] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[14,24] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[16,24] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[17,24] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[3,25] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[4,25] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[7,25] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[11,25] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[12,25] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[13,25] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[14,25] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[16,25] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[17,25] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[3,26] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[4,26] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[7,26] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[11,26] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[12,26] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[13,26] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[14,26] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[16,26] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[17,26] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[3,27] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[4,27] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[7,27] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[11,27] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[12,27] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[13,27] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[14,27] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[16,27] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[17,27] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[18,27] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[3,28] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[4,28] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[7,28] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[11,28] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[12,28] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[13,28] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[14,28] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[16,28] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[17,28] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[18,28] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[3,29] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[4,29] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[7,29] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[11,29] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[12,29] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[13,29] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[14,29] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[16,29] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[17,29] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[3,30] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[4,30] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[7,30] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[11,30] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[12,30] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[13,30] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[14,30] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[16,30] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[17,30] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[3,31] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[4,31] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[7,31] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[11,31] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[12,31] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[13,31] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[14,31] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[16,31] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[17,31] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[3,32] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[11,32] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[12,32] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[13,32] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[16,32] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[17,32] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[3,33] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[12,33] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[13,33] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[17,33] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[3,34] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[12,34] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[13,34] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[3,35] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[12,35] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[13,35] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[3,36] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[12,36] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[13,36] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[3,37] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[13,37] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[3,38] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[13,38] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[3,39] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## S[13,39] 0.972200 0.01111064 0.9470893 0.9656048 0.973443 0.9801735 0.9897591 1.000883 4500
## beta 3.636509 0.42828160 2.8847878 3.3348363 3.601548 3.9007101 4.5710693 1.000930 4500
## deviance 43.521831 1.48866734 42.5113204 42.6147795 42.955421 43.7840171 47.8513413 1.001939 1600
results$BUGSoutput$summary[grepl("beta",rownames(results$BUGSoutput$summary)),
c("mean", "sd", "2.5%","97.5%","Rhat", "n.eff")]
## mean sd 2.5% 97.5% Rhat n.eff
## 3.6365091 0.4282816 2.8847878 4.5710693 1.0009298 4500.0000000
# get just the means
results$BUGSoutput$mean
## $S
## [,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]
## [1,] 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722
## [2,] 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722
## [3,] 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722
## [4,] 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722
## [5,] 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722
## [6,] 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722
## [7,] 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722
## [8,] 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722
## [9,] 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722
## [10,] 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722
## [11,] 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722
## [12,] 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722
## [13,] 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722
## [14,] 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722
## [15,] 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722
## [16,] 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722
## [17,] 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722
## [18,] 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722
## [,28] [,29] [,30] [,31] [,32] [,33] [,34] [,35] [,36] [,37] [,38] [,39]
## [1,] 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722
## [2,] 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722
## [3,] 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722
## [4,] 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722
## [5,] 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722
## [6,] 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722
## [7,] 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722
## [8,] 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722
## [9,] 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722
## [10,] 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722
## [11,] 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722
## [12,] 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722
## [13,] 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722
## [14,] 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722
## [15,] 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722
## [16,] 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722
## [17,] 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722
## [18,] 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722 0.9722
##
## $beta
## [1] 3.636509
##
## $deviance
## [1] 43.52183
results$BUGSoutput$mean$parm
## NULL
# the results$BUGSoutput$sims.array is a 3-d object [iterations, chains, variables]
# that can be used for posterior plots, diagnostics etc
parm <- "beta"
dim(results$BUGSoutput$sims.array)
## [1] 1500 3 225
#results$BUGSoutput$sims.array[1:5,,]
results$BUGSoutput$sims.array[1:5,1,parm,drop=FALSE]
## , , beta
##
## [,1]
## [1,] 3.295938
## [2,] 3.607257
## [3,] 3.432225
## [4,] 3.643801
## [5,] 3.494060
# the results$BUGSoutput$sims.matrix is a 2-d object [iterations, variables] with chains stacked
# on top of each other
dim(results$BUGSoutput$sims.matrix)
## [1] 4500 225
#results$BUGSoutput$sims.matrix[1:5,]
results$BUGSoutput$sims.matrix[1:5,parm,drop=FALSE]
## beta
## [1,] 2.678632
## [2,] 3.280311
## [3,] 3.909404
## [4,] 3.530549
## [5,] 3.057137
# make a posterior density plot of the parameter
plotdata <- data.frame(parm=results$BUGSoutput$sims.matrix[,parm], stringsAsFactors=FALSE)
head(plotdata)
## parm
## 1 2.678632
## 2 3.280311
## 3 3.909404
## 4 3.530549
## 5 3.057137
## 6 3.586671
postplot.parm <- ggplot2::ggplot( data=plotdata, aes(x=parm, y=..density..))+
geom_histogram(alpha=0.3)+
geom_density()+
ggtitle(paste("Posterior density plot for ", parm))
postplot.parm
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#ggsave(plot=postplot.parm, file='...posterior.png', h=4, w=6, units="in", dpi=300)
# make a trace plot (notice we use the sims.array here)
plotdata <- data.frame(results$BUGSoutput$sims.array[,,parm,drop=FALSE], stringsAsFactors=FALSE)
plotdata$iteration <- 1:nrow(plotdata)
head(plotdata)
## X1.beta X2.beta X3.beta iteration
## 1 3.295938 3.782221 3.263847 1
## 2 3.607257 3.127559 2.876422 2
## 3 3.432225 3.343433 3.089447 3
## 4 3.643801 3.565407 3.000498 4
## 5 3.494060 4.021540 3.810255 5
## 6 4.014514 3.170581 3.674709 6
# convert from wide to long format
plotdata2 <- reshape2::melt(data=plotdata,
id.vars="iteration",
measure.vars=paste("X",1:results$BUGSoutput$n.chains,".",parm,sep=""),
variable.name="chain",
value.name=parm)
head(plotdata2)
## iteration chain beta
## 1 1 X1.beta 3.295938
## 2 2 X1.beta 3.607257
## 3 3 X1.beta 3.432225
## 4 4 X1.beta 3.643801
## 5 5 X1.beta 3.494060
## 6 6 X1.beta 4.014514
traceplot.parm <- ggplot2::ggplot(data=plotdata2, aes_string(x="iteration", y=parm, color="chain"))+
ggtitle("Trace plot")+
geom_line(alpha=.2)
traceplot.parm
#ggsave(plot=traceplot.parm, file='....-trace-parm.png', h=4, w=6, units="in", dpi=300)
# autocorrelation plot
# First compute the autocorrelation plot
acf.parm <-acf( results$BUGSoutput$sims.matrix[,parm,drop=FALSE], plot=FALSE)
acf.parm
##
## Autocorrelations of series 'results$BUGSoutput$sims.matrix[, parm, drop = FALSE]', by lag
##
## 0 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
## 1.000 0.024 0.007 -0.024 -0.019 0.005 -0.012 0.008 0.012 0.008 0.046 0.007 0.027 -0.026 0.008 -0.021 -0.010 0.010 0.012 0.022 -0.006 -0.009 0.021 0.010 -0.017 -0.036 0.034 -0.008
## 28 29 30 31 32 33 34 35 36
## 0.017 -0.011 -0.005 0.001 -0.007 0.029 0.019 -0.001 0.004
acfplot.parm <- ggplot(data=with(acf.parm, data.frame(lag, acf)), aes(x = lag, y = acf)) +
ggtitle(paste("Autocorrelation plot for ",parm))+
geom_hline(aes(yintercept = 0)) +
geom_segment(aes(xend = lag, yend = 0))
acfplot.parm
#ggsave(plot=acfplot.parm, file="...-acf-parm.png",h=4, w=6, units="in", dpi=300)
# You can save the results information to an r-object for retreival later using a load# command
#save(list=c("results"), file="results.Rdata")
#load(file="results.Rdata")