# Single Species Single Season Occupancy models
# Salamander Data with some missing data
# Using JAGS
# The only real change is the you need to delete any site-survey data with missing
# data before sending to JAGS. This is currently done automatically.
#
# Be sure to specify the missing data code when reading in the data
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(car)
## Loading required package: carData
library(coda)
library(ggplot2)
library(reshape2)
options(width=400) # make html output wider
# The BUGS model is specified as a text file.
# The model file.
# The cat() command is used to save the model to the working directory.
# Notice that you CANNOT have any " (double quotes) in the bugs code
# between the start and end of the cat("...",) command.
# Inputs to the model are
# Nsites - number of sites
# Nvisits - (max) number of visits over all sites.
# Nsites.visits - number of sites x number of visits
# if there is missing data (no visits), simply drop the corresponding row
# History - vector of 1 or 0 corresponding to Site-Visit pair
# Site - vector indicating which site the row corresponds to
# Visit - vector indicating which visit the row corresponds to
#
# dmatrix.psi - design matrix for psi
# Nbeta.psi - number of columns in design matrix for psi
# dmatrix.p - design matrix for p
# Nbeta.p - number of columns of design matrix for p
#
cat(file="model.txt", "
############################################################
model {
# estimate psi for each site from the design matrix
for(i in 1:Nsites){
logit(psi[i]) = inprod( dmatrix.psi[i, 1:Nbeta.psi], beta.psi[1:Nbeta.psi])
}
# estimate p for each observation
for(i in 1:Nsites.visits){
logit(p[i]) = inprod( dmatrix.p[i, 1:Nbeta.p], beta.p[1:Nbeta.p])
}
# set up the state model, i.e. is the site actually occupied or not
for(i in 1:Nsites){
z[i] ~ dbern(psi[i])
}
# the observation model.
for(j in 1:Nsites.visits){
p.z[j] <- z[Site[j]]*p[j]
History[j] ~ dbern(p.z[j])
}
# priors on the betas
beta.psi[1] ~ dnorm(0, .25) # intercept we want basically flat on regular scale
for(i in 2:Nbeta.psi){
beta.psi[i] ~ dnorm(0, .0001)
}
beta.p[1] ~ dnorm(0, .25)
for(i in 2:Nbeta.p){
beta.p[i] ~ dnorm(0, .0001)
}
# derived variables
# number of occupied sites
occ.sites <- sum(z[1:Nsites])
# belief that psi is above some value
prob.psi.greater.50 <- ifelse( psi[1] > 0.5, 1, 0)
}
") # End of the model
# get the data read in
# Data for detections should be a data frame with rows corresponding to sites
# and columns to visits.
# The usual 1=detected; 0=not detected; NA=not visited is used.
input.data <- readxl::read_excel("../salamander.xls",
sheet="MissingData",
na="-",
col_names=FALSE) # notice no column names in row 1 of data file.
head(input.data)
## # A tibble: 6 x 5
## X__1 X__2 X__3 X__4 X__5
## <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 0 0 0 1 1
## 2 0 1 NA 0 0
## 3 0 NA NA 0 0
## 4 1 NA NA 1 0
## 5 0 NA NA 0 0
## 6 0 0 NA 0 0
# Extract the history records
input.history <- input.data[, 1:5] # the history extracted
head(input.history)
## # A tibble: 6 x 5
## X__1 X__2 X__3 X__4 X__5
## <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 0 0 0 1 1
## 2 0 1 NA 0 0
## 3 0 NA NA 0 0
## 4 1 NA NA 1 0
## 5 0 NA NA 0 0
## 6 0 0 NA 0 0
# do some basic checks on your data
# e.g. check number of sites; number of visits etc
nrow(input.history)
## [1] 39
ncol(input.history)
## [1] 5
range(input.history, na.rm=TRUE) # check that all values are either 0 or 1
## [1] 0 1
sum(is.na(input.history)) # are there any missing values?
## [1] 14
# get the data in the right format
Survey <- data.frame(
History = as.vector(unlist(input.history)), # stacks the columns
Site = rep(1:nrow(input.history), ncol(input.history)),
Visit = rep(1:ncol(input.history), each=nrow(input.history)),
stringsAsFactors=FALSE)
Survey[1:10,]
## History Site Visit
## 1 0 1 1
## 2 0 2 1
## 3 0 3 1
## 4 1 4 1
## 5 0 5 1
## 6 0 6 1
## 7 0 7 1
## 8 0 8 1
## 9 0 9 1
## 10 1 10 1
# Remove any rows with missing history value (i.e missing)
sum(is.na(Survey$History))
## [1] 14
dim(Survey)
## [1] 195 3
Survey <- Survey[!is.na(Survey$History),]
dim(Survey)
## [1] 181 3
Nsites <- nrow(input.history)
Nvisits <- ncol(input.history)
Nsites.visits <- nrow(Survey)
# Get the design matrix for model p~1 psi~1
dmatrix.psi <- model.matrix(~1, data=input.data) # constant psi
Nbeta.psi <- ncol(dmatrix.psi)
dmatrix.p <- model.matrix(~1, data=Survey)
Nbeta.p <- ncol(dmatrix.p)
model.name <- "pdot-psidot"
# Get the design matrix for model p~time psi~1
dmatrix.psi <- model.matrix(~1, data=input.data) # constant psi
Nbeta.psi <- ncol(dmatrix.psi)
dmatrix.p <- model.matrix(~as.factor(Visit), data=Survey)
Nbeta.p <- ncol(dmatrix.p)
model.name <- "pt-psidot"
# Get the design matrix for model p~custom psi~1
dmatrix.psi <- model.matrix(~1, data=input.data) # constant psi
Nbeta.psi <- ncol(dmatrix.psi)
Survey$D <- car::recode(Survey$Visit,
"1:2=1; 3:5=2")
xtabs(~D+Visit, data=Survey)
## Visit
## D 1 2 3 4 5
## 1 39 36 0 0 0
## 2 0 0 32 38 36
dmatrix.p <- model.matrix(~as.factor(D), data=Survey)
Nbeta.p <- ncol(dmatrix.p)
model.name <- "pcustom-psidot"
# The datalist will be passed to JAGS with the names of the data
# values.
data.list <- list(Nsites=Nsites,
Nvisit=Nvisits,
Nsites.visits=Nsites.visits,
History = Survey$History,
Site = Survey$Site,
Visit = Survey$Visit,
dmatrix.psi=dmatrix.psi, Nbeta.psi=Nbeta.psi,
dmatrix.p =dmatrix.p, Nbeta.p =Nbeta.p)
# check the list
data.list
## $Nsites
## [1] 39
##
## $Nvisit
## [1] 5
##
## $Nsites.visits
## [1] 181
##
## $History
## [1] 0 0 0 1 0 0 0 0 0 1 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 0 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 1 0 0 0 0 0 0 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 1 0 0 0 0 0 0 0 0 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1
##
## $Site
## [1] 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 1 2 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 1 7 8 9 10 11 12 14 15 16 17 18 19 20 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24
## [132] 25 27 28 29 30 31 32 33 34 35 36 37 38 39 1 2 3 4 5 6 8 9 10 11 12 13 14 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 33 34 35 36 37 38 39
##
## $Visit
## [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
##
## $dmatrix.psi
## (Intercept)
## 1 1
## 2 1
## 3 1
## 4 1
## 5 1
## 6 1
## 7 1
## 8 1
## 9 1
## 10 1
## 11 1
## 12 1
## 13 1
## 14 1
## 15 1
## 16 1
## 17 1
## 18 1
## 19 1
## 20 1
## 21 1
## 22 1
## 23 1
## 24 1
## 25 1
## 26 1
## 27 1
## 28 1
## 29 1
## 30 1
## 31 1
## 32 1
## 33 1
## 34 1
## 35 1
## 36 1
## 37 1
## 38 1
## 39 1
## attr(,"assign")
## [1] 0
##
## $Nbeta.psi
## [1] 1
##
## $dmatrix.p
## (Intercept) as.factor(D)2
## 1 1 0
## 2 1 0
## 3 1 0
## 4 1 0
## 5 1 0
## 6 1 0
## 7 1 0
## 8 1 0
## 9 1 0
## 10 1 0
## 11 1 0
## 12 1 0
## 13 1 0
## 14 1 0
## 15 1 0
## 16 1 0
## 17 1 0
## 18 1 0
## 19 1 0
## 20 1 0
## 21 1 0
## 22 1 0
## 23 1 0
## 24 1 0
## 25 1 0
## 26 1 0
## 27 1 0
## 28 1 0
## 29 1 0
## 30 1 0
## 31 1 0
## 32 1 0
## 33 1 0
## 34 1 0
## 35 1 0
## 36 1 0
## 37 1 0
## 38 1 0
## 39 1 0
## 40 1 0
## 41 1 0
## 45 1 0
## 46 1 0
## 47 1 0
## 48 1 0
## 49 1 0
## 50 1 0
## 51 1 0
## 52 1 0
## 53 1 0
## 54 1 0
## 55 1 0
## 56 1 0
## 57 1 0
## 58 1 0
## 59 1 0
## 60 1 0
## 61 1 0
## 62 1 0
## 63 1 0
## 64 1 0
## 65 1 0
## 66 1 0
## 67 1 0
## 68 1 0
## 69 1 0
## 70 1 0
## 71 1 0
## 72 1 0
## 73 1 0
## 74 1 0
## 75 1 0
## 76 1 0
## 77 1 0
## 78 1 0
## 79 1 1
## 85 1 1
## 86 1 1
## 87 1 1
## 88 1 1
## 89 1 1
## 90 1 1
## 92 1 1
## 93 1 1
## 94 1 1
## 95 1 1
## 96 1 1
## 97 1 1
## 98 1 1
## 100 1 1
## 101 1 1
## 102 1 1
## 103 1 1
## 104 1 1
## 105 1 1
## 106 1 1
## 107 1 1
## 108 1 1
## 109 1 1
## 110 1 1
## 111 1 1
## 112 1 1
## 113 1 1
## 114 1 1
## 115 1 1
## 116 1 1
## 117 1 1
## 118 1 1
## 119 1 1
## 120 1 1
## 121 1 1
## 122 1 1
## 123 1 1
## 124 1 1
## 125 1 1
## 126 1 1
## 127 1 1
## 128 1 1
## 129 1 1
## 130 1 1
## 131 1 1
## 132 1 1
## 133 1 1
## 134 1 1
## 135 1 1
## 136 1 1
## 137 1 1
## 138 1 1
## 139 1 1
## 140 1 1
## 141 1 1
## 142 1 1
## 144 1 1
## 145 1 1
## 146 1 1
## 147 1 1
## 148 1 1
## 149 1 1
## 150 1 1
## 151 1 1
## 152 1 1
## 153 1 1
## 154 1 1
## 155 1 1
## 156 1 1
## 157 1 1
## 158 1 1
## 159 1 1
## 160 1 1
## 161 1 1
## 162 1 1
## 164 1 1
## 165 1 1
## 166 1 1
## 167 1 1
## 168 1 1
## 169 1 1
## 170 1 1
## 172 1 1
## 173 1 1
## 174 1 1
## 175 1 1
## 176 1 1
## 177 1 1
## 178 1 1
## 179 1 1
## 180 1 1
## 181 1 1
## 182 1 1
## 183 1 1
## 184 1 1
## 185 1 1
## 186 1 1
## 187 1 1
## 189 1 1
## 190 1 1
## 191 1 1
## 192 1 1
## 193 1 1
## 194 1 1
## 195 1 1
## attr(,"assign")
## [1] 0 1
## attr(,"contrasts")
## attr(,"contrasts")$`as.factor(D)`
## [1] "contr.treatment"
##
##
## $Nbeta.p
## [1] 2
# Next create the initial values.
# If you are using more than one chain, you need to create a function
# that returns initial values for each chain.
# We define the initial value of z as 1 if any visit resulted in a detection, other wise 0
init.z <- apply(input.history, 1, max, na.rm=TRUE)
# we will start at the same initial starting point for each chain even though this
# is not recommended.
init.list <- list(
list(z=init.z, beta.psi=rep(0,Nbeta.psi), beta.p=rep(0,Nbeta.p) ),
list(z=init.z, beta.psi=rep(0,Nbeta.psi), beta.p=rep(0,Nbeta.p) ),
list(z=init.z, beta.psi=rep(0,Nbeta.psi), beta.p=rep(0,Nbeta.p) )
) # end of list of lists of initial values
# Next create the list of parameters to monitor.
# The deviance is automatically monitored.
#
monitor.list <- c("z","occ.sites", "prob.psi.greater.50",
"psi", "p",
"beta.psi", "beta.p") # parameters to monitor
# Finally, the actual call to JAGS
set.seed(234234) # intitalize seed for MCMC
results <- R2jags::jags(
data =data.list, # list of data variables
inits =init.list, # list/function for initial values
parameters=monitor.list,# list of parameters to monitor
model.file="model.txt", # file with bugs model
n.chains=3,
n.iter =5000, # total iterations INCLUDING burn in
n.burnin=2000, # number of burning iterations
n.thin=2, # how much to thin
DIC=TRUE # is DIC to be computed?
)
## module glm loaded
## Warning in jags.model(model.file, data = data, inits = init.values, n.chains = n.chains, : Unused variable "Nvisit" in data
## Warning in jags.model(model.file, data = data, inits = init.values, n.chains = n.chains, : Unused variable "Visit" in data
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 181
## Unobserved stochastic nodes: 42
## Total graph size: 1084
##
## Initializing model
#######################################
# extract some of the usual stuff and use R code directly
# use the standard print method
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" "mean" "sd" "median" "root.short" "long.short" "dimension.short" "indexes.short" "last.values" "program" "model.file" "isDIC" "DICbyR"
## [23] "pD" "DIC"
# get the summary table
results$BUGSoutput$summary
## mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff
## beta.p[1] -1.83572077 0.50330204 -2.87007814 -2.1750834 -1.821147703 -1.4844257 -0.8866843 1.002257 1300
## beta.p[2] 1.35879456 0.54646693 0.28677305 1.0011160 1.340884921 1.7077881 2.4626376 1.001523 2300
## beta.psi 0.04604139 0.55944841 -0.84122208 -0.3040633 -0.008909423 0.3156635 1.2592332 1.002657 4500
## deviance 100.90999273 10.69374388 84.31738986 92.9207679 99.229969602 107.2373105 126.2717924 1.001046 4500
## occ.sites 19.70466667 3.78878042 15.00000000 17.0000000 19.000000000 21.0000000 29.5207718 1.001079 4500
## p[1] 0.14803092 0.06130976 0.05365268 0.1020104 0.139296215 0.1847599 0.2917944 1.002344 1200
## p[2] 0.14803092 0.06130976 0.05365268 0.1020104 0.139296215 0.1847599 0.2917944 1.002344 1200
## p[3] 0.14803092 0.06130976 0.05365268 0.1020104 0.139296215 0.1847599 0.2917944 1.002344 1200
## p[4] 0.14803092 0.06130976 0.05365268 0.1020104 0.139296215 0.1847599 0.2917944 1.002344 1200
## p[5] 0.14803092 0.06130976 0.05365268 0.1020104 0.139296215 0.1847599 0.2917944 1.002344 1200
## p[6] 0.14803092 0.06130976 0.05365268 0.1020104 0.139296215 0.1847599 0.2917944 1.002344 1200
## p[7] 0.14803092 0.06130976 0.05365268 0.1020104 0.139296215 0.1847599 0.2917944 1.002344 1200
## p[8] 0.14803092 0.06130976 0.05365268 0.1020104 0.139296215 0.1847599 0.2917944 1.002344 1200
## p[9] 0.14803092 0.06130976 0.05365268 0.1020104 0.139296215 0.1847599 0.2917944 1.002344 1200
## p[10] 0.14803092 0.06130976 0.05365268 0.1020104 0.139296215 0.1847599 0.2917944 1.002344 1200
## p[11] 0.14803092 0.06130976 0.05365268 0.1020104 0.139296215 0.1847599 0.2917944 1.002344 1200
## p[12] 0.14803092 0.06130976 0.05365268 0.1020104 0.139296215 0.1847599 0.2917944 1.002344 1200
## p[13] 0.14803092 0.06130976 0.05365268 0.1020104 0.139296215 0.1847599 0.2917944 1.002344 1200
## p[14] 0.14803092 0.06130976 0.05365268 0.1020104 0.139296215 0.1847599 0.2917944 1.002344 1200
## p[15] 0.14803092 0.06130976 0.05365268 0.1020104 0.139296215 0.1847599 0.2917944 1.002344 1200
## p[16] 0.14803092 0.06130976 0.05365268 0.1020104 0.139296215 0.1847599 0.2917944 1.002344 1200
## p[17] 0.14803092 0.06130976 0.05365268 0.1020104 0.139296215 0.1847599 0.2917944 1.002344 1200
## p[18] 0.14803092 0.06130976 0.05365268 0.1020104 0.139296215 0.1847599 0.2917944 1.002344 1200
## p[19] 0.14803092 0.06130976 0.05365268 0.1020104 0.139296215 0.1847599 0.2917944 1.002344 1200
## p[20] 0.14803092 0.06130976 0.05365268 0.1020104 0.139296215 0.1847599 0.2917944 1.002344 1200
## p[21] 0.14803092 0.06130976 0.05365268 0.1020104 0.139296215 0.1847599 0.2917944 1.002344 1200
## p[22] 0.14803092 0.06130976 0.05365268 0.1020104 0.139296215 0.1847599 0.2917944 1.002344 1200
## p[23] 0.14803092 0.06130976 0.05365268 0.1020104 0.139296215 0.1847599 0.2917944 1.002344 1200
## p[24] 0.14803092 0.06130976 0.05365268 0.1020104 0.139296215 0.1847599 0.2917944 1.002344 1200
## p[25] 0.14803092 0.06130976 0.05365268 0.1020104 0.139296215 0.1847599 0.2917944 1.002344 1200
## p[26] 0.14803092 0.06130976 0.05365268 0.1020104 0.139296215 0.1847599 0.2917944 1.002344 1200
## p[27] 0.14803092 0.06130976 0.05365268 0.1020104 0.139296215 0.1847599 0.2917944 1.002344 1200
## p[28] 0.14803092 0.06130976 0.05365268 0.1020104 0.139296215 0.1847599 0.2917944 1.002344 1200
## p[29] 0.14803092 0.06130976 0.05365268 0.1020104 0.139296215 0.1847599 0.2917944 1.002344 1200
## p[30] 0.14803092 0.06130976 0.05365268 0.1020104 0.139296215 0.1847599 0.2917944 1.002344 1200
## p[31] 0.14803092 0.06130976 0.05365268 0.1020104 0.139296215 0.1847599 0.2917944 1.002344 1200
## p[32] 0.14803092 0.06130976 0.05365268 0.1020104 0.139296215 0.1847599 0.2917944 1.002344 1200
## p[33] 0.14803092 0.06130976 0.05365268 0.1020104 0.139296215 0.1847599 0.2917944 1.002344 1200
## p[34] 0.14803092 0.06130976 0.05365268 0.1020104 0.139296215 0.1847599 0.2917944 1.002344 1200
## p[35] 0.14803092 0.06130976 0.05365268 0.1020104 0.139296215 0.1847599 0.2917944 1.002344 1200
## p[36] 0.14803092 0.06130976 0.05365268 0.1020104 0.139296215 0.1847599 0.2917944 1.002344 1200
## p[37] 0.14803092 0.06130976 0.05365268 0.1020104 0.139296215 0.1847599 0.2917944 1.002344 1200
## p[38] 0.14803092 0.06130976 0.05365268 0.1020104 0.139296215 0.1847599 0.2917944 1.002344 1200
## p[39] 0.14803092 0.06130976 0.05365268 0.1020104 0.139296215 0.1847599 0.2917944 1.002344 1200
## p[40] 0.14803092 0.06130976 0.05365268 0.1020104 0.139296215 0.1847599 0.2917944 1.002344 1200
## p[41] 0.14803092 0.06130976 0.05365268 0.1020104 0.139296215 0.1847599 0.2917944 1.002344 1200
## p[42] 0.14803092 0.06130976 0.05365268 0.1020104 0.139296215 0.1847599 0.2917944 1.002344 1200
## p[43] 0.14803092 0.06130976 0.05365268 0.1020104 0.139296215 0.1847599 0.2917944 1.002344 1200
## p[44] 0.14803092 0.06130976 0.05365268 0.1020104 0.139296215 0.1847599 0.2917944 1.002344 1200
## p[45] 0.14803092 0.06130976 0.05365268 0.1020104 0.139296215 0.1847599 0.2917944 1.002344 1200
## p[46] 0.14803092 0.06130976 0.05365268 0.1020104 0.139296215 0.1847599 0.2917944 1.002344 1200
## p[47] 0.14803092 0.06130976 0.05365268 0.1020104 0.139296215 0.1847599 0.2917944 1.002344 1200
## p[48] 0.14803092 0.06130976 0.05365268 0.1020104 0.139296215 0.1847599 0.2917944 1.002344 1200
## p[49] 0.14803092 0.06130976 0.05365268 0.1020104 0.139296215 0.1847599 0.2917944 1.002344 1200
## p[50] 0.14803092 0.06130976 0.05365268 0.1020104 0.139296215 0.1847599 0.2917944 1.002344 1200
## p[51] 0.14803092 0.06130976 0.05365268 0.1020104 0.139296215 0.1847599 0.2917944 1.002344 1200
## p[52] 0.14803092 0.06130976 0.05365268 0.1020104 0.139296215 0.1847599 0.2917944 1.002344 1200
## p[53] 0.14803092 0.06130976 0.05365268 0.1020104 0.139296215 0.1847599 0.2917944 1.002344 1200
## p[54] 0.14803092 0.06130976 0.05365268 0.1020104 0.139296215 0.1847599 0.2917944 1.002344 1200
## p[55] 0.14803092 0.06130976 0.05365268 0.1020104 0.139296215 0.1847599 0.2917944 1.002344 1200
## p[56] 0.14803092 0.06130976 0.05365268 0.1020104 0.139296215 0.1847599 0.2917944 1.002344 1200
## p[57] 0.14803092 0.06130976 0.05365268 0.1020104 0.139296215 0.1847599 0.2917944 1.002344 1200
## p[58] 0.14803092 0.06130976 0.05365268 0.1020104 0.139296215 0.1847599 0.2917944 1.002344 1200
## p[59] 0.14803092 0.06130976 0.05365268 0.1020104 0.139296215 0.1847599 0.2917944 1.002344 1200
## p[60] 0.14803092 0.06130976 0.05365268 0.1020104 0.139296215 0.1847599 0.2917944 1.002344 1200
## p[61] 0.14803092 0.06130976 0.05365268 0.1020104 0.139296215 0.1847599 0.2917944 1.002344 1200
## p[62] 0.14803092 0.06130976 0.05365268 0.1020104 0.139296215 0.1847599 0.2917944 1.002344 1200
## p[63] 0.14803092 0.06130976 0.05365268 0.1020104 0.139296215 0.1847599 0.2917944 1.002344 1200
## p[64] 0.14803092 0.06130976 0.05365268 0.1020104 0.139296215 0.1847599 0.2917944 1.002344 1200
## p[65] 0.14803092 0.06130976 0.05365268 0.1020104 0.139296215 0.1847599 0.2917944 1.002344 1200
## p[66] 0.14803092 0.06130976 0.05365268 0.1020104 0.139296215 0.1847599 0.2917944 1.002344 1200
## p[67] 0.14803092 0.06130976 0.05365268 0.1020104 0.139296215 0.1847599 0.2917944 1.002344 1200
## p[68] 0.14803092 0.06130976 0.05365268 0.1020104 0.139296215 0.1847599 0.2917944 1.002344 1200
## p[69] 0.14803092 0.06130976 0.05365268 0.1020104 0.139296215 0.1847599 0.2917944 1.002344 1200
## p[70] 0.14803092 0.06130976 0.05365268 0.1020104 0.139296215 0.1847599 0.2917944 1.002344 1200
## p[71] 0.14803092 0.06130976 0.05365268 0.1020104 0.139296215 0.1847599 0.2917944 1.002344 1200
## p[72] 0.14803092 0.06130976 0.05365268 0.1020104 0.139296215 0.1847599 0.2917944 1.002344 1200
## p[73] 0.14803092 0.06130976 0.05365268 0.1020104 0.139296215 0.1847599 0.2917944 1.002344 1200
## p[74] 0.14803092 0.06130976 0.05365268 0.1020104 0.139296215 0.1847599 0.2917944 1.002344 1200
## p[75] 0.14803092 0.06130976 0.05365268 0.1020104 0.139296215 0.1847599 0.2917944 1.002344 1200
## p[76] 0.38732380 0.09216694 0.21718393 0.3226051 0.385864752 0.4494486 0.5771274 1.000805 4500
## p[77] 0.38732380 0.09216694 0.21718393 0.3226051 0.385864752 0.4494486 0.5771274 1.000805 4500
## p[78] 0.38732380 0.09216694 0.21718393 0.3226051 0.385864752 0.4494486 0.5771274 1.000805 4500
## p[79] 0.38732380 0.09216694 0.21718393 0.3226051 0.385864752 0.4494486 0.5771274 1.000805 4500
## p[80] 0.38732380 0.09216694 0.21718393 0.3226051 0.385864752 0.4494486 0.5771274 1.000805 4500
## p[81] 0.38732380 0.09216694 0.21718393 0.3226051 0.385864752 0.4494486 0.5771274 1.000805 4500
## p[82] 0.38732380 0.09216694 0.21718393 0.3226051 0.385864752 0.4494486 0.5771274 1.000805 4500
## p[83] 0.38732380 0.09216694 0.21718393 0.3226051 0.385864752 0.4494486 0.5771274 1.000805 4500
## p[84] 0.38732380 0.09216694 0.21718393 0.3226051 0.385864752 0.4494486 0.5771274 1.000805 4500
## p[85] 0.38732380 0.09216694 0.21718393 0.3226051 0.385864752 0.4494486 0.5771274 1.000805 4500
## p[86] 0.38732380 0.09216694 0.21718393 0.3226051 0.385864752 0.4494486 0.5771274 1.000805 4500
## p[87] 0.38732380 0.09216694 0.21718393 0.3226051 0.385864752 0.4494486 0.5771274 1.000805 4500
## p[88] 0.38732380 0.09216694 0.21718393 0.3226051 0.385864752 0.4494486 0.5771274 1.000805 4500
## p[89] 0.38732380 0.09216694 0.21718393 0.3226051 0.385864752 0.4494486 0.5771274 1.000805 4500
## p[90] 0.38732380 0.09216694 0.21718393 0.3226051 0.385864752 0.4494486 0.5771274 1.000805 4500
## p[91] 0.38732380 0.09216694 0.21718393 0.3226051 0.385864752 0.4494486 0.5771274 1.000805 4500
## p[92] 0.38732380 0.09216694 0.21718393 0.3226051 0.385864752 0.4494486 0.5771274 1.000805 4500
## p[93] 0.38732380 0.09216694 0.21718393 0.3226051 0.385864752 0.4494486 0.5771274 1.000805 4500
## p[94] 0.38732380 0.09216694 0.21718393 0.3226051 0.385864752 0.4494486 0.5771274 1.000805 4500
## p[95] 0.38732380 0.09216694 0.21718393 0.3226051 0.385864752 0.4494486 0.5771274 1.000805 4500
## p[96] 0.38732380 0.09216694 0.21718393 0.3226051 0.385864752 0.4494486 0.5771274 1.000805 4500
## p[97] 0.38732380 0.09216694 0.21718393 0.3226051 0.385864752 0.4494486 0.5771274 1.000805 4500
## p[98] 0.38732380 0.09216694 0.21718393 0.3226051 0.385864752 0.4494486 0.5771274 1.000805 4500
## p[99] 0.38732380 0.09216694 0.21718393 0.3226051 0.385864752 0.4494486 0.5771274 1.000805 4500
## p[100] 0.38732380 0.09216694 0.21718393 0.3226051 0.385864752 0.4494486 0.5771274 1.000805 4500
## p[101] 0.38732380 0.09216694 0.21718393 0.3226051 0.385864752 0.4494486 0.5771274 1.000805 4500
## p[102] 0.38732380 0.09216694 0.21718393 0.3226051 0.385864752 0.4494486 0.5771274 1.000805 4500
## p[103] 0.38732380 0.09216694 0.21718393 0.3226051 0.385864752 0.4494486 0.5771274 1.000805 4500
## p[104] 0.38732380 0.09216694 0.21718393 0.3226051 0.385864752 0.4494486 0.5771274 1.000805 4500
## p[105] 0.38732380 0.09216694 0.21718393 0.3226051 0.385864752 0.4494486 0.5771274 1.000805 4500
## p[106] 0.38732380 0.09216694 0.21718393 0.3226051 0.385864752 0.4494486 0.5771274 1.000805 4500
## p[107] 0.38732380 0.09216694 0.21718393 0.3226051 0.385864752 0.4494486 0.5771274 1.000805 4500
## p[108] 0.38732380 0.09216694 0.21718393 0.3226051 0.385864752 0.4494486 0.5771274 1.000805 4500
## p[109] 0.38732380 0.09216694 0.21718393 0.3226051 0.385864752 0.4494486 0.5771274 1.000805 4500
## p[110] 0.38732380 0.09216694 0.21718393 0.3226051 0.385864752 0.4494486 0.5771274 1.000805 4500
## p[111] 0.38732380 0.09216694 0.21718393 0.3226051 0.385864752 0.4494486 0.5771274 1.000805 4500
## p[112] 0.38732380 0.09216694 0.21718393 0.3226051 0.385864752 0.4494486 0.5771274 1.000805 4500
## p[113] 0.38732380 0.09216694 0.21718393 0.3226051 0.385864752 0.4494486 0.5771274 1.000805 4500
## p[114] 0.38732380 0.09216694 0.21718393 0.3226051 0.385864752 0.4494486 0.5771274 1.000805 4500
## p[115] 0.38732380 0.09216694 0.21718393 0.3226051 0.385864752 0.4494486 0.5771274 1.000805 4500
## p[116] 0.38732380 0.09216694 0.21718393 0.3226051 0.385864752 0.4494486 0.5771274 1.000805 4500
## p[117] 0.38732380 0.09216694 0.21718393 0.3226051 0.385864752 0.4494486 0.5771274 1.000805 4500
## p[118] 0.38732380 0.09216694 0.21718393 0.3226051 0.385864752 0.4494486 0.5771274 1.000805 4500
## p[119] 0.38732380 0.09216694 0.21718393 0.3226051 0.385864752 0.4494486 0.5771274 1.000805 4500
## p[120] 0.38732380 0.09216694 0.21718393 0.3226051 0.385864752 0.4494486 0.5771274 1.000805 4500
## p[121] 0.38732380 0.09216694 0.21718393 0.3226051 0.385864752 0.4494486 0.5771274 1.000805 4500
## p[122] 0.38732380 0.09216694 0.21718393 0.3226051 0.385864752 0.4494486 0.5771274 1.000805 4500
## p[123] 0.38732380 0.09216694 0.21718393 0.3226051 0.385864752 0.4494486 0.5771274 1.000805 4500
## p[124] 0.38732380 0.09216694 0.21718393 0.3226051 0.385864752 0.4494486 0.5771274 1.000805 4500
## p[125] 0.38732380 0.09216694 0.21718393 0.3226051 0.385864752 0.4494486 0.5771274 1.000805 4500
## p[126] 0.38732380 0.09216694 0.21718393 0.3226051 0.385864752 0.4494486 0.5771274 1.000805 4500
## p[127] 0.38732380 0.09216694 0.21718393 0.3226051 0.385864752 0.4494486 0.5771274 1.000805 4500
## p[128] 0.38732380 0.09216694 0.21718393 0.3226051 0.385864752 0.4494486 0.5771274 1.000805 4500
## p[129] 0.38732380 0.09216694 0.21718393 0.3226051 0.385864752 0.4494486 0.5771274 1.000805 4500
## p[130] 0.38732380 0.09216694 0.21718393 0.3226051 0.385864752 0.4494486 0.5771274 1.000805 4500
## p[131] 0.38732380 0.09216694 0.21718393 0.3226051 0.385864752 0.4494486 0.5771274 1.000805 4500
## p[132] 0.38732380 0.09216694 0.21718393 0.3226051 0.385864752 0.4494486 0.5771274 1.000805 4500
## p[133] 0.38732380 0.09216694 0.21718393 0.3226051 0.385864752 0.4494486 0.5771274 1.000805 4500
## p[134] 0.38732380 0.09216694 0.21718393 0.3226051 0.385864752 0.4494486 0.5771274 1.000805 4500
## p[135] 0.38732380 0.09216694 0.21718393 0.3226051 0.385864752 0.4494486 0.5771274 1.000805 4500
## p[136] 0.38732380 0.09216694 0.21718393 0.3226051 0.385864752 0.4494486 0.5771274 1.000805 4500
## p[137] 0.38732380 0.09216694 0.21718393 0.3226051 0.385864752 0.4494486 0.5771274 1.000805 4500
## p[138] 0.38732380 0.09216694 0.21718393 0.3226051 0.385864752 0.4494486 0.5771274 1.000805 4500
## p[139] 0.38732380 0.09216694 0.21718393 0.3226051 0.385864752 0.4494486 0.5771274 1.000805 4500
## p[140] 0.38732380 0.09216694 0.21718393 0.3226051 0.385864752 0.4494486 0.5771274 1.000805 4500
## p[141] 0.38732380 0.09216694 0.21718393 0.3226051 0.385864752 0.4494486 0.5771274 1.000805 4500
## p[142] 0.38732380 0.09216694 0.21718393 0.3226051 0.385864752 0.4494486 0.5771274 1.000805 4500
## p[143] 0.38732380 0.09216694 0.21718393 0.3226051 0.385864752 0.4494486 0.5771274 1.000805 4500
## p[144] 0.38732380 0.09216694 0.21718393 0.3226051 0.385864752 0.4494486 0.5771274 1.000805 4500
## p[145] 0.38732380 0.09216694 0.21718393 0.3226051 0.385864752 0.4494486 0.5771274 1.000805 4500
## p[146] 0.38732380 0.09216694 0.21718393 0.3226051 0.385864752 0.4494486 0.5771274 1.000805 4500
## p[147] 0.38732380 0.09216694 0.21718393 0.3226051 0.385864752 0.4494486 0.5771274 1.000805 4500
## p[148] 0.38732380 0.09216694 0.21718393 0.3226051 0.385864752 0.4494486 0.5771274 1.000805 4500
## p[149] 0.38732380 0.09216694 0.21718393 0.3226051 0.385864752 0.4494486 0.5771274 1.000805 4500
## p[150] 0.38732380 0.09216694 0.21718393 0.3226051 0.385864752 0.4494486 0.5771274 1.000805 4500
## p[151] 0.38732380 0.09216694 0.21718393 0.3226051 0.385864752 0.4494486 0.5771274 1.000805 4500
## p[152] 0.38732380 0.09216694 0.21718393 0.3226051 0.385864752 0.4494486 0.5771274 1.000805 4500
## p[153] 0.38732380 0.09216694 0.21718393 0.3226051 0.385864752 0.4494486 0.5771274 1.000805 4500
## p[154] 0.38732380 0.09216694 0.21718393 0.3226051 0.385864752 0.4494486 0.5771274 1.000805 4500
## p[155] 0.38732380 0.09216694 0.21718393 0.3226051 0.385864752 0.4494486 0.5771274 1.000805 4500
## p[156] 0.38732380 0.09216694 0.21718393 0.3226051 0.385864752 0.4494486 0.5771274 1.000805 4500
## p[157] 0.38732380 0.09216694 0.21718393 0.3226051 0.385864752 0.4494486 0.5771274 1.000805 4500
## p[158] 0.38732380 0.09216694 0.21718393 0.3226051 0.385864752 0.4494486 0.5771274 1.000805 4500
## p[159] 0.38732380 0.09216694 0.21718393 0.3226051 0.385864752 0.4494486 0.5771274 1.000805 4500
## p[160] 0.38732380 0.09216694 0.21718393 0.3226051 0.385864752 0.4494486 0.5771274 1.000805 4500
## p[161] 0.38732380 0.09216694 0.21718393 0.3226051 0.385864752 0.4494486 0.5771274 1.000805 4500
## p[162] 0.38732380 0.09216694 0.21718393 0.3226051 0.385864752 0.4494486 0.5771274 1.000805 4500
## p[163] 0.38732380 0.09216694 0.21718393 0.3226051 0.385864752 0.4494486 0.5771274 1.000805 4500
## p[164] 0.38732380 0.09216694 0.21718393 0.3226051 0.385864752 0.4494486 0.5771274 1.000805 4500
## p[165] 0.38732380 0.09216694 0.21718393 0.3226051 0.385864752 0.4494486 0.5771274 1.000805 4500
## p[166] 0.38732380 0.09216694 0.21718393 0.3226051 0.385864752 0.4494486 0.5771274 1.000805 4500
## p[167] 0.38732380 0.09216694 0.21718393 0.3226051 0.385864752 0.4494486 0.5771274 1.000805 4500
## p[168] 0.38732380 0.09216694 0.21718393 0.3226051 0.385864752 0.4494486 0.5771274 1.000805 4500
## p[169] 0.38732380 0.09216694 0.21718393 0.3226051 0.385864752 0.4494486 0.5771274 1.000805 4500
## p[170] 0.38732380 0.09216694 0.21718393 0.3226051 0.385864752 0.4494486 0.5771274 1.000805 4500
## p[171] 0.38732380 0.09216694 0.21718393 0.3226051 0.385864752 0.4494486 0.5771274 1.000805 4500
## p[172] 0.38732380 0.09216694 0.21718393 0.3226051 0.385864752 0.4494486 0.5771274 1.000805 4500
## p[173] 0.38732380 0.09216694 0.21718393 0.3226051 0.385864752 0.4494486 0.5771274 1.000805 4500
## p[174] 0.38732380 0.09216694 0.21718393 0.3226051 0.385864752 0.4494486 0.5771274 1.000805 4500
## p[175] 0.38732380 0.09216694 0.21718393 0.3226051 0.385864752 0.4494486 0.5771274 1.000805 4500
## p[176] 0.38732380 0.09216694 0.21718393 0.3226051 0.385864752 0.4494486 0.5771274 1.000805 4500
## p[177] 0.38732380 0.09216694 0.21718393 0.3226051 0.385864752 0.4494486 0.5771274 1.000805 4500
## p[178] 0.38732380 0.09216694 0.21718393 0.3226051 0.385864752 0.4494486 0.5771274 1.000805 4500
## p[179] 0.38732380 0.09216694 0.21718393 0.3226051 0.385864752 0.4494486 0.5771274 1.000805 4500
## p[180] 0.38732380 0.09216694 0.21718393 0.3226051 0.385864752 0.4494486 0.5771274 1.000805 4500
## p[181] 0.38732380 0.09216694 0.21718393 0.3226051 0.385864752 0.4494486 0.5771274 1.000805 4500
## prob.psi.greater.50 0.49288889 0.50000499 0.00000000 0.0000000 0.000000000 1.0000000 1.0000000 1.001149 4100
## psi[1] 0.50786163 0.12084376 0.30127745 0.4245645 0.497772657 0.5782671 0.7788937 1.001579 4500
## psi[2] 0.50786163 0.12084376 0.30127745 0.4245645 0.497772657 0.5782671 0.7788937 1.001579 4500
## psi[3] 0.50786163 0.12084376 0.30127745 0.4245645 0.497772657 0.5782671 0.7788937 1.001579 4500
## psi[4] 0.50786163 0.12084376 0.30127745 0.4245645 0.497772657 0.5782671 0.7788937 1.001579 4500
## psi[5] 0.50786163 0.12084376 0.30127745 0.4245645 0.497772657 0.5782671 0.7788937 1.001579 4500
## psi[6] 0.50786163 0.12084376 0.30127745 0.4245645 0.497772657 0.5782671 0.7788937 1.001579 4500
## psi[7] 0.50786163 0.12084376 0.30127745 0.4245645 0.497772657 0.5782671 0.7788937 1.001579 4500
## psi[8] 0.50786163 0.12084376 0.30127745 0.4245645 0.497772657 0.5782671 0.7788937 1.001579 4500
## psi[9] 0.50786163 0.12084376 0.30127745 0.4245645 0.497772657 0.5782671 0.7788937 1.001579 4500
## psi[10] 0.50786163 0.12084376 0.30127745 0.4245645 0.497772657 0.5782671 0.7788937 1.001579 4500
## psi[11] 0.50786163 0.12084376 0.30127745 0.4245645 0.497772657 0.5782671 0.7788937 1.001579 4500
## psi[12] 0.50786163 0.12084376 0.30127745 0.4245645 0.497772657 0.5782671 0.7788937 1.001579 4500
## psi[13] 0.50786163 0.12084376 0.30127745 0.4245645 0.497772657 0.5782671 0.7788937 1.001579 4500
## psi[14] 0.50786163 0.12084376 0.30127745 0.4245645 0.497772657 0.5782671 0.7788937 1.001579 4500
## psi[15] 0.50786163 0.12084376 0.30127745 0.4245645 0.497772657 0.5782671 0.7788937 1.001579 4500
## psi[16] 0.50786163 0.12084376 0.30127745 0.4245645 0.497772657 0.5782671 0.7788937 1.001579 4500
## psi[17] 0.50786163 0.12084376 0.30127745 0.4245645 0.497772657 0.5782671 0.7788937 1.001579 4500
## psi[18] 0.50786163 0.12084376 0.30127745 0.4245645 0.497772657 0.5782671 0.7788937 1.001579 4500
## psi[19] 0.50786163 0.12084376 0.30127745 0.4245645 0.497772657 0.5782671 0.7788937 1.001579 4500
## psi[20] 0.50786163 0.12084376 0.30127745 0.4245645 0.497772657 0.5782671 0.7788937 1.001579 4500
## psi[21] 0.50786163 0.12084376 0.30127745 0.4245645 0.497772657 0.5782671 0.7788937 1.001579 4500
## psi[22] 0.50786163 0.12084376 0.30127745 0.4245645 0.497772657 0.5782671 0.7788937 1.001579 4500
## psi[23] 0.50786163 0.12084376 0.30127745 0.4245645 0.497772657 0.5782671 0.7788937 1.001579 4500
## psi[24] 0.50786163 0.12084376 0.30127745 0.4245645 0.497772657 0.5782671 0.7788937 1.001579 4500
## psi[25] 0.50786163 0.12084376 0.30127745 0.4245645 0.497772657 0.5782671 0.7788937 1.001579 4500
## psi[26] 0.50786163 0.12084376 0.30127745 0.4245645 0.497772657 0.5782671 0.7788937 1.001579 4500
## psi[27] 0.50786163 0.12084376 0.30127745 0.4245645 0.497772657 0.5782671 0.7788937 1.001579 4500
## psi[28] 0.50786163 0.12084376 0.30127745 0.4245645 0.497772657 0.5782671 0.7788937 1.001579 4500
## psi[29] 0.50786163 0.12084376 0.30127745 0.4245645 0.497772657 0.5782671 0.7788937 1.001579 4500
## psi[30] 0.50786163 0.12084376 0.30127745 0.4245645 0.497772657 0.5782671 0.7788937 1.001579 4500
## psi[31] 0.50786163 0.12084376 0.30127745 0.4245645 0.497772657 0.5782671 0.7788937 1.001579 4500
## psi[32] 0.50786163 0.12084376 0.30127745 0.4245645 0.497772657 0.5782671 0.7788937 1.001579 4500
## psi[33] 0.50786163 0.12084376 0.30127745 0.4245645 0.497772657 0.5782671 0.7788937 1.001579 4500
## psi[34] 0.50786163 0.12084376 0.30127745 0.4245645 0.497772657 0.5782671 0.7788937 1.001579 4500
## psi[35] 0.50786163 0.12084376 0.30127745 0.4245645 0.497772657 0.5782671 0.7788937 1.001579 4500
## psi[36] 0.50786163 0.12084376 0.30127745 0.4245645 0.497772657 0.5782671 0.7788937 1.001579 4500
## psi[37] 0.50786163 0.12084376 0.30127745 0.4245645 0.497772657 0.5782671 0.7788937 1.001579 4500
## psi[38] 0.50786163 0.12084376 0.30127745 0.4245645 0.497772657 0.5782671 0.7788937 1.001579 4500
## psi[39] 0.50786163 0.12084376 0.30127745 0.4245645 0.497772657 0.5782671 0.7788937 1.001579 4500
## z[1] 1.00000000 0.00000000 1.00000000 1.0000000 1.000000000 1.0000000 1.0000000 1.000000 1
## z[2] 1.00000000 0.00000000 1.00000000 1.0000000 1.000000000 1.0000000 1.0000000 1.000000 1
## z[3] 0.26688889 0.44238298 0.00000000 0.0000000 0.000000000 1.0000000 1.0000000 1.002076 1400
## z[4] 1.00000000 0.00000000 1.00000000 1.0000000 1.000000000 1.0000000 1.0000000 1.000000 1
## z[5] 0.26866667 0.44331542 0.00000000 0.0000000 0.000000000 1.0000000 1.0000000 1.002633 1000
## z[6] 0.24911111 0.43254635 0.00000000 0.0000000 0.000000000 0.0000000 1.0000000 1.000703 4500
## z[7] 1.00000000 0.00000000 1.00000000 1.0000000 1.000000000 1.0000000 1.0000000 1.000000 1
## z[8] 1.00000000 0.00000000 1.00000000 1.0000000 1.000000000 1.0000000 1.0000000 1.000000 1
## z[9] 1.00000000 0.00000000 1.00000000 1.0000000 1.000000000 1.0000000 1.0000000 1.000000 1
## z[10] 1.00000000 0.00000000 1.00000000 1.0000000 1.000000000 1.0000000 1.0000000 1.000000 1
## z[11] 1.00000000 0.00000000 1.00000000 1.0000000 1.000000000 1.0000000 1.0000000 1.000000 1
## z[12] 1.00000000 0.00000000 1.00000000 1.0000000 1.000000000 1.0000000 1.0000000 1.000000 1
## z[13] 1.00000000 0.00000000 1.00000000 1.0000000 1.000000000 1.0000000 1.0000000 1.000000 1
## z[14] 1.00000000 0.00000000 1.00000000 1.0000000 1.000000000 1.0000000 1.0000000 1.000000 1
## z[15] 0.24933333 0.43267519 0.00000000 0.0000000 0.000000000 0.0000000 1.0000000 1.000980 4500
## z[16] 0.17088889 0.37645367 0.00000000 0.0000000 0.000000000 0.0000000 1.0000000 1.001168 4000
## z[17] 0.17866667 0.38311553 0.00000000 0.0000000 0.000000000 0.0000000 1.0000000 1.000964 4500
## z[18] 0.17466667 0.37972393 0.00000000 0.0000000 0.000000000 0.0000000 1.0000000 1.000794 4500
## z[19] 0.17311111 0.37838535 0.00000000 0.0000000 0.000000000 0.0000000 1.0000000 1.000968 4500
## z[20] 0.17266667 0.37800085 0.00000000 0.0000000 0.000000000 0.0000000 1.0000000 1.001016 4500
## z[21] 0.23711111 0.42535825 0.00000000 0.0000000 0.000000000 0.0000000 1.0000000 1.000698 4500
## z[22] 0.17911111 0.38348794 0.00000000 0.0000000 0.000000000 0.0000000 1.0000000 1.001049 4500
## z[23] 0.17555556 0.38048387 0.00000000 0.0000000 0.000000000 0.0000000 1.0000000 1.000756 4500
## z[24] 0.17444444 0.37953339 0.00000000 0.0000000 0.000000000 0.0000000 1.0000000 1.001184 3900
## z[25] 0.17466667 0.37972393 0.00000000 0.0000000 0.000000000 0.0000000 1.0000000 1.000689 4500
## z[26] 0.24022222 0.42726581 0.00000000 0.0000000 0.000000000 0.0000000 1.0000000 1.001143 4200
## z[27] 0.16733333 0.37331469 0.00000000 0.0000000 0.000000000 0.0000000 1.0000000 1.000741 4500
## z[28] 0.17355556 0.37876893 0.00000000 0.0000000 0.000000000 0.0000000 1.0000000 1.000932 4500
## z[29] 0.17977778 0.38404492 0.00000000 0.0000000 0.000000000 0.0000000 1.0000000 1.000673 4500
## z[30] 0.17044444 0.37606457 0.00000000 0.0000000 0.000000000 0.0000000 1.0000000 1.000758 4500
## z[31] 0.17888889 0.38330184 0.00000000 0.0000000 0.000000000 0.0000000 1.0000000 1.001566 2200
## z[32] 0.22955556 0.42059376 0.00000000 0.0000000 0.000000000 0.0000000 1.0000000 1.001242 3500
## z[33] 0.17244444 0.37780826 0.00000000 0.0000000 0.000000000 0.0000000 1.0000000 1.001208 3700
## z[34] 0.17266667 0.37800085 0.00000000 0.0000000 0.000000000 0.0000000 1.0000000 1.001188 3800
## z[35] 0.17488889 0.37991425 0.00000000 0.0000000 0.000000000 0.0000000 1.0000000 1.000739 4500
## z[36] 1.00000000 0.00000000 1.00000000 1.0000000 1.000000000 1.0000000 1.0000000 1.000000 1
## z[37] 1.00000000 0.00000000 1.00000000 1.0000000 1.000000000 1.0000000 1.0000000 1.000000 1
## z[38] 1.00000000 0.00000000 1.00000000 1.0000000 1.000000000 1.0000000 1.0000000 1.000000 1
## z[39] 1.00000000 0.00000000 1.00000000 1.0000000 1.000000000 1.0000000 1.0000000 1.000000 1
#results$BUGSoutput$summary[,c("mean", "sd", "2.5%","97.5%","Rhat", "n.eff")]
#results$BUGSoutput$summary[,c("mean", "sd")]
# get just the means
results$BUGSoutput$mean
## $beta.p
## [1] -1.835721 1.358795
##
## $beta.psi
## [1] 0.04604139
##
## $deviance
## [1] 100.91
##
## $occ.sites
## [1] 19.70467
##
## $p
## [1] 0.1480309 0.1480309 0.1480309 0.1480309 0.1480309 0.1480309 0.1480309 0.1480309 0.1480309 0.1480309 0.1480309 0.1480309 0.1480309 0.1480309 0.1480309 0.1480309 0.1480309 0.1480309 0.1480309 0.1480309 0.1480309 0.1480309 0.1480309 0.1480309 0.1480309 0.1480309 0.1480309 0.1480309 0.1480309 0.1480309 0.1480309 0.1480309 0.1480309 0.1480309 0.1480309 0.1480309 0.1480309 0.1480309 0.1480309
## [40] 0.1480309 0.1480309 0.1480309 0.1480309 0.1480309 0.1480309 0.1480309 0.1480309 0.1480309 0.1480309 0.1480309 0.1480309 0.1480309 0.1480309 0.1480309 0.1480309 0.1480309 0.1480309 0.1480309 0.1480309 0.1480309 0.1480309 0.1480309 0.1480309 0.1480309 0.1480309 0.1480309 0.1480309 0.1480309 0.1480309 0.1480309 0.1480309 0.1480309 0.1480309 0.1480309 0.1480309 0.3873238 0.3873238 0.3873238
## [79] 0.3873238 0.3873238 0.3873238 0.3873238 0.3873238 0.3873238 0.3873238 0.3873238 0.3873238 0.3873238 0.3873238 0.3873238 0.3873238 0.3873238 0.3873238 0.3873238 0.3873238 0.3873238 0.3873238 0.3873238 0.3873238 0.3873238 0.3873238 0.3873238 0.3873238 0.3873238 0.3873238 0.3873238 0.3873238 0.3873238 0.3873238 0.3873238 0.3873238 0.3873238 0.3873238 0.3873238 0.3873238 0.3873238 0.3873238
## [118] 0.3873238 0.3873238 0.3873238 0.3873238 0.3873238 0.3873238 0.3873238 0.3873238 0.3873238 0.3873238 0.3873238 0.3873238 0.3873238 0.3873238 0.3873238 0.3873238 0.3873238 0.3873238 0.3873238 0.3873238 0.3873238 0.3873238 0.3873238 0.3873238 0.3873238 0.3873238 0.3873238 0.3873238 0.3873238 0.3873238 0.3873238 0.3873238 0.3873238 0.3873238 0.3873238 0.3873238 0.3873238 0.3873238 0.3873238
## [157] 0.3873238 0.3873238 0.3873238 0.3873238 0.3873238 0.3873238 0.3873238 0.3873238 0.3873238 0.3873238 0.3873238 0.3873238 0.3873238 0.3873238 0.3873238 0.3873238 0.3873238 0.3873238 0.3873238 0.3873238 0.3873238 0.3873238 0.3873238 0.3873238 0.3873238
##
## $prob.psi.greater.50
## [1] 0.4928889
##
## $psi
## [1] 0.5078616 0.5078616 0.5078616 0.5078616 0.5078616 0.5078616 0.5078616 0.5078616 0.5078616 0.5078616 0.5078616 0.5078616 0.5078616 0.5078616 0.5078616 0.5078616 0.5078616 0.5078616 0.5078616 0.5078616 0.5078616 0.5078616 0.5078616 0.5078616 0.5078616 0.5078616 0.5078616 0.5078616 0.5078616 0.5078616 0.5078616 0.5078616 0.5078616 0.5078616 0.5078616 0.5078616 0.5078616 0.5078616 0.5078616
##
## $z
## [1] 1.0000000 1.0000000 0.2668889 1.0000000 0.2686667 0.2491111 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 0.2493333 0.1708889 0.1786667 0.1746667 0.1731111 0.1726667 0.2371111 0.1791111 0.1755556 0.1744444 0.1746667 0.2402222 0.1673333 0.1735556 0.1797778 0.1704444 0.1788889 0.2295556 0.1724444 0.1726667 0.1748889 1.0000000 1.0000000 1.0000000 1.0000000
results$BUGSoutput$mean$psi
## [1] 0.5078616 0.5078616 0.5078616 0.5078616 0.5078616 0.5078616 0.5078616 0.5078616 0.5078616 0.5078616 0.5078616 0.5078616 0.5078616 0.5078616 0.5078616 0.5078616 0.5078616 0.5078616 0.5078616 0.5078616 0.5078616 0.5078616 0.5078616 0.5078616 0.5078616 0.5078616 0.5078616 0.5078616 0.5078616 0.5078616 0.5078616 0.5078616 0.5078616 0.5078616 0.5078616 0.5078616 0.5078616 0.5078616 0.5078616
# the results$BUGSoutput$sims.array is a 3-d object [iterations, chains, variables]
dim(results$BUGSoutput$sims.array)
## [1] 1500 3 265
results$BUGSoutput$sims.array[1:5,1,1:10]
## beta.p[1] beta.p[2] beta.psi deviance occ.sites p[1] p[2] p[3] p[4] p[5]
## [1,] -1.484610 0.9889579 -0.9157445 85.52712 15 0.1847321 0.1847321 0.1847321 0.1847321 0.1847321
## [2,] -1.405711 1.1672109 -0.3318074 94.28395 18 0.1969114 0.1969114 0.1969114 0.1969114 0.1969114
## [3,] -1.438183 0.3681355 0.1766649 106.23510 20 0.1918269 0.1918269 0.1918269 0.1918269 0.1918269
## [4,] -1.412309 0.6270413 0.5079229 110.55165 23 0.1958702 0.1958702 0.1958702 0.1958702 0.1958702
## [5,] -1.552404 0.9934135 -0.6962203 95.26030 18 0.1747393 0.1747393 0.1747393 0.1747393 0.1747393
results$BUGSoutput$sims.array[1:5,1,"psi[1]", drop=FALSE]
## , , psi[1]
##
## [,1]
## [1,] 0.2858258
## [2,] 0.4178009
## [3,] 0.5440517
## [4,] 0.6243194
## [5,] 0.3326508
# 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 265
results$BUGSoutput$sims.matrix[1:5,1:10]
## beta.p[1] beta.p[2] beta.psi deviance occ.sites p[1] p[2] p[3] p[4] p[5]
## [1,] -1.221165 1.1679592 -0.5376714 101.48679 19 0.2277315 0.2277315 0.2277315 0.2277315 0.2277315
## [2,] -1.492434 0.8998334 0.1699359 92.09811 17 0.1835566 0.1835566 0.1835566 0.1835566 0.1835566
## [3,] -2.000118 1.8317410 0.1031967 95.39125 18 0.1191905 0.1191905 0.1191905 0.1191905 0.1191905
## [4,] -1.678803 1.1621381 0.3177595 100.38832 20 0.1572540 0.1572540 0.1572540 0.1572540 0.1572540
## [5,] -2.075308 2.4228936 -0.1869708 89.79337 16 0.1115200 0.1115200 0.1115200 0.1115200 0.1115200
results$BUGSoutput$sims.matrix[1:5,"psi[1]", drop=FALSE]
## psi[1]
## [1,] 0.3687294
## [2,] 0.5423820
## [3,] 0.5257763
## [4,] 0.5787781
## [5,] 0.4533930
# make a posterior density plot
plotdata <- data.frame(parm=results$BUGSoutput$sims.matrix[,"psi[1]"], stringsAsFactors=FALSE)
head(plotdata)
## parm
## 1 0.3687294
## 2 0.5423820
## 3 0.5257763
## 4 0.5787781
## 5 0.4533930
## 6 0.4034227
postplot.parm <- ggplot2::ggplot( data=plotdata, aes(x=parm, y=..density..))+
geom_histogram(alpha=0.3)+
geom_density()+
ggtitle("Posterior density plot for psi[1]")
postplot.parm
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
ggsave(plot=postplot.parm,
file=paste('psi-posterior-',model.name,'.png',sep=""), h=4, w=6, units="in", dpi=300)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
# make a trace plot (notice we use the sims.array here)
plotdata <- data.frame(psi=results$BUGSoutput$sims.array[,,"psi[1]"], stringsAsFactors=FALSE)
plotdata$iteration <- 1:nrow(plotdata)
head(plotdata)
## psi.1 psi.2 psi.3 iteration
## 1 0.2858258 0.4219379 0.4830863 1
## 2 0.4178009 0.4413266 0.4419501 2
## 3 0.5440517 0.5807583 0.3246133 3
## 4 0.6243194 0.4335565 0.4351200 4
## 5 0.3326508 0.4967923 0.5496172 5
## 6 0.2995774 0.2867484 0.5943521 6
# convert from wide to long format
plotdata2 <- reshape2:::melt.data.frame(data=plotdata,
id.vars="iteration",
measure.vars=paste("psi",1:results$BUGSoutput$n.chains,sep="."),
variable.name="chain",
value.name='psi')
head(plotdata2)
## iteration chain psi
## 1 1 psi.1 0.2858258
## 2 2 psi.1 0.4178009
## 3 3 psi.1 0.5440517
## 4 4 psi.1 0.6243194
## 5 5 psi.1 0.3326508
## 6 6 psi.1 0.2995774
traceplot.parm <- ggplot2::ggplot(data=plotdata2, aes(x=iteration, y=psi, color=chain))+
ggtitle("Trace plot for psi")+
geom_line(alpha=.2)
traceplot.parm
ggsave(plot=traceplot.parm,
file=paste('psi-trace-',model.name,'.png',sep=""), h=4, w=6, units="in", dpi=300)
# autocorrelation plot
# First compute the autocorrelation plot
acf.parm <-acf( results$BUGSoutput$sims.matrix[,"psi[1]"], plot=FALSE)
acf.parm
##
## Autocorrelations of series 'results$BUGSoutput$sims.matrix[, "psi[1]"]', 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 28 29 30 31 32 33 34 35 36
## 1.000 0.001 -0.006 0.000 -0.037 0.011 -0.028 0.008 0.004 0.029 -0.019 0.016 -0.004 -0.021 -0.005 0.005 0.016 0.005 0.007 0.003 0.011 0.024 -0.009 0.026 0.011 0.024 0.008 0.005 0.038 -0.017 -0.001 -0.005 0.002 0.047 0.018 -0.025 -0.001
acfplot.parm <- ggplot(data=with(acf.parm, data.frame(lag, acf)), aes(x = lag, y = acf)) +
ggtitle("Autocorrelation plot for psi")+
geom_hline(aes(yintercept = 0)) +
geom_segment(aes(xend = lag, yend = 0))
acfplot.parm
ggsave(plot=acfplot.parm,
file=paste("psi-acf-", model.name,".png",sep=""),h=4, w=6, units="in", dpi=300)