# Single Species Single Season Occupancy models
# Salamander Data
# using JAGS
# 2018-08-15 Code contributed by Carl James Schwarz (cschwarz.stat.sfu.cs@gmail.com)
library("R2jags") # used for call to JAGS
## Loading required package: rjags
## Loading required package: coda
## Linked to JAGS 4.3.2
## 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="CompleteData",
col_names=FALSE) # notice no column names in row 1 of data file.
## New names:
## • `` -> `...1`
## • `` -> `...2`
## • `` -> `...3`
## • `` -> `...4`
## • `` -> `...5`
head(input.data)
## # A tibble: 6 × 5
## ...1 ...2 ...3 ...4 ...5
## <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 0 0 0 1 1
## 2 0 1 0 0 0
## 3 0 1 0 0 0
## 4 1 1 1 1 0
## 5 0 0 1 0 0
## 6 0 0 1 0 0
# Extract the history records
input.history <- input.data[, 1:5] # the history extracted
head(input.history)
## # A tibble: 6 × 5
## ...1 ...2 ...3 ...4 ...5
## <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 0 0 0 1 1
## 2 0 1 0 0 0
## 3 0 1 0 0 0
## 4 1 1 1 1 0
## 5 0 0 1 0 0
## 6 0 0 1 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] 0
# 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] 0
Survey <- Survey[!is.na(Survey$History),]
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 39 0 0 0
## 2 0 0 39 39 39
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] 195
##
## $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 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 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 0 1 1 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 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 0 1 1 0 0 1 0 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 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 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 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 3 4 5 6 7 8 9 10 11 12 13 14
## [132] 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 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
##
## $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 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 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 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 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
## 42 1 0
## 43 1 0
## 44 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
## 80 1 1
## 81 1 1
## 82 1 1
## 83 1 1
## 84 1 1
## 85 1 1
## 86 1 1
## 87 1 1
## 88 1 1
## 89 1 1
## 90 1 1
## 91 1 1
## 92 1 1
## 93 1 1
## 94 1 1
## 95 1 1
## 96 1 1
## 97 1 1
## 98 1 1
## 99 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
## 143 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
## 163 1 1
## 164 1 1
## 165 1 1
## 166 1 1
## 167 1 1
## 168 1 1
## 169 1 1
## 170 1 1
## 171 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
## 188 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: 195
## Unobserved stochastic nodes: 42
## Total graph size: 1154
##
## 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.7350115 0.45471403 -2.67336414 -2.0322590 -1.7175884 -1.4260594 -0.8918359 1.003311 760
## beta.p[2] 1.0044813 0.49602644 0.04536792 0.6696130 0.9948865 1.3182221 1.9928630 1.004493 520
## beta.psi 0.4964369 0.60231027 -0.43316009 0.0881291 0.4169006 0.7896980 1.9429472 1.004361 990
## deviance 130.7800617 11.55545430 112.25539357 121.7866896 129.8162840 138.0986086 156.4732395 1.003223 1100
## occ.sites 23.9204444 4.05358417 18.00000000 21.0000000 23.0000000 26.0000000 34.0000000 1.003570 940
## p[1] 0.1588179 0.05854104 0.06456349 0.1158573 0.1521821 0.1937134 0.2907309 1.003294 760
## p[2] 0.1588179 0.05854104 0.06456349 0.1158573 0.1521821 0.1937134 0.2907309 1.003294 760
## p[3] 0.1588179 0.05854104 0.06456349 0.1158573 0.1521821 0.1937134 0.2907309 1.003294 760
## p[4] 0.1588179 0.05854104 0.06456349 0.1158573 0.1521821 0.1937134 0.2907309 1.003294 760
## p[5] 0.1588179 0.05854104 0.06456349 0.1158573 0.1521821 0.1937134 0.2907309 1.003294 760
## p[6] 0.1588179 0.05854104 0.06456349 0.1158573 0.1521821 0.1937134 0.2907309 1.003294 760
## p[7] 0.1588179 0.05854104 0.06456349 0.1158573 0.1521821 0.1937134 0.2907309 1.003294 760
## p[8] 0.1588179 0.05854104 0.06456349 0.1158573 0.1521821 0.1937134 0.2907309 1.003294 760
## p[9] 0.1588179 0.05854104 0.06456349 0.1158573 0.1521821 0.1937134 0.2907309 1.003294 760
## p[10] 0.1588179 0.05854104 0.06456349 0.1158573 0.1521821 0.1937134 0.2907309 1.003294 760
## p[11] 0.1588179 0.05854104 0.06456349 0.1158573 0.1521821 0.1937134 0.2907309 1.003294 760
## p[12] 0.1588179 0.05854104 0.06456349 0.1158573 0.1521821 0.1937134 0.2907309 1.003294 760
## p[13] 0.1588179 0.05854104 0.06456349 0.1158573 0.1521821 0.1937134 0.2907309 1.003294 760
## p[14] 0.1588179 0.05854104 0.06456349 0.1158573 0.1521821 0.1937134 0.2907309 1.003294 760
## p[15] 0.1588179 0.05854104 0.06456349 0.1158573 0.1521821 0.1937134 0.2907309 1.003294 760
## p[16] 0.1588179 0.05854104 0.06456349 0.1158573 0.1521821 0.1937134 0.2907309 1.003294 760
## p[17] 0.1588179 0.05854104 0.06456349 0.1158573 0.1521821 0.1937134 0.2907309 1.003294 760
## p[18] 0.1588179 0.05854104 0.06456349 0.1158573 0.1521821 0.1937134 0.2907309 1.003294 760
## p[19] 0.1588179 0.05854104 0.06456349 0.1158573 0.1521821 0.1937134 0.2907309 1.003294 760
## p[20] 0.1588179 0.05854104 0.06456349 0.1158573 0.1521821 0.1937134 0.2907309 1.003294 760
## p[21] 0.1588179 0.05854104 0.06456349 0.1158573 0.1521821 0.1937134 0.2907309 1.003294 760
## p[22] 0.1588179 0.05854104 0.06456349 0.1158573 0.1521821 0.1937134 0.2907309 1.003294 760
## p[23] 0.1588179 0.05854104 0.06456349 0.1158573 0.1521821 0.1937134 0.2907309 1.003294 760
## p[24] 0.1588179 0.05854104 0.06456349 0.1158573 0.1521821 0.1937134 0.2907309 1.003294 760
## p[25] 0.1588179 0.05854104 0.06456349 0.1158573 0.1521821 0.1937134 0.2907309 1.003294 760
## p[26] 0.1588179 0.05854104 0.06456349 0.1158573 0.1521821 0.1937134 0.2907309 1.003294 760
## p[27] 0.1588179 0.05854104 0.06456349 0.1158573 0.1521821 0.1937134 0.2907309 1.003294 760
## p[28] 0.1588179 0.05854104 0.06456349 0.1158573 0.1521821 0.1937134 0.2907309 1.003294 760
## p[29] 0.1588179 0.05854104 0.06456349 0.1158573 0.1521821 0.1937134 0.2907309 1.003294 760
## p[30] 0.1588179 0.05854104 0.06456349 0.1158573 0.1521821 0.1937134 0.2907309 1.003294 760
## p[31] 0.1588179 0.05854104 0.06456349 0.1158573 0.1521821 0.1937134 0.2907309 1.003294 760
## p[32] 0.1588179 0.05854104 0.06456349 0.1158573 0.1521821 0.1937134 0.2907309 1.003294 760
## p[33] 0.1588179 0.05854104 0.06456349 0.1158573 0.1521821 0.1937134 0.2907309 1.003294 760
## p[34] 0.1588179 0.05854104 0.06456349 0.1158573 0.1521821 0.1937134 0.2907309 1.003294 760
## p[35] 0.1588179 0.05854104 0.06456349 0.1158573 0.1521821 0.1937134 0.2907309 1.003294 760
## p[36] 0.1588179 0.05854104 0.06456349 0.1158573 0.1521821 0.1937134 0.2907309 1.003294 760
## p[37] 0.1588179 0.05854104 0.06456349 0.1158573 0.1521821 0.1937134 0.2907309 1.003294 760
## p[38] 0.1588179 0.05854104 0.06456349 0.1158573 0.1521821 0.1937134 0.2907309 1.003294 760
## p[39] 0.1588179 0.05854104 0.06456349 0.1158573 0.1521821 0.1937134 0.2907309 1.003294 760
## p[40] 0.1588179 0.05854104 0.06456349 0.1158573 0.1521821 0.1937134 0.2907309 1.003294 760
## p[41] 0.1588179 0.05854104 0.06456349 0.1158573 0.1521821 0.1937134 0.2907309 1.003294 760
## p[42] 0.1588179 0.05854104 0.06456349 0.1158573 0.1521821 0.1937134 0.2907309 1.003294 760
## p[43] 0.1588179 0.05854104 0.06456349 0.1158573 0.1521821 0.1937134 0.2907309 1.003294 760
## p[44] 0.1588179 0.05854104 0.06456349 0.1158573 0.1521821 0.1937134 0.2907309 1.003294 760
## p[45] 0.1588179 0.05854104 0.06456349 0.1158573 0.1521821 0.1937134 0.2907309 1.003294 760
## p[46] 0.1588179 0.05854104 0.06456349 0.1158573 0.1521821 0.1937134 0.2907309 1.003294 760
## p[47] 0.1588179 0.05854104 0.06456349 0.1158573 0.1521821 0.1937134 0.2907309 1.003294 760
## p[48] 0.1588179 0.05854104 0.06456349 0.1158573 0.1521821 0.1937134 0.2907309 1.003294 760
## p[49] 0.1588179 0.05854104 0.06456349 0.1158573 0.1521821 0.1937134 0.2907309 1.003294 760
## p[50] 0.1588179 0.05854104 0.06456349 0.1158573 0.1521821 0.1937134 0.2907309 1.003294 760
## p[51] 0.1588179 0.05854104 0.06456349 0.1158573 0.1521821 0.1937134 0.2907309 1.003294 760
## p[52] 0.1588179 0.05854104 0.06456349 0.1158573 0.1521821 0.1937134 0.2907309 1.003294 760
## p[53] 0.1588179 0.05854104 0.06456349 0.1158573 0.1521821 0.1937134 0.2907309 1.003294 760
## p[54] 0.1588179 0.05854104 0.06456349 0.1158573 0.1521821 0.1937134 0.2907309 1.003294 760
## p[55] 0.1588179 0.05854104 0.06456349 0.1158573 0.1521821 0.1937134 0.2907309 1.003294 760
## p[56] 0.1588179 0.05854104 0.06456349 0.1158573 0.1521821 0.1937134 0.2907309 1.003294 760
## p[57] 0.1588179 0.05854104 0.06456349 0.1158573 0.1521821 0.1937134 0.2907309 1.003294 760
## p[58] 0.1588179 0.05854104 0.06456349 0.1158573 0.1521821 0.1937134 0.2907309 1.003294 760
## p[59] 0.1588179 0.05854104 0.06456349 0.1158573 0.1521821 0.1937134 0.2907309 1.003294 760
## p[60] 0.1588179 0.05854104 0.06456349 0.1158573 0.1521821 0.1937134 0.2907309 1.003294 760
## p[61] 0.1588179 0.05854104 0.06456349 0.1158573 0.1521821 0.1937134 0.2907309 1.003294 760
## p[62] 0.1588179 0.05854104 0.06456349 0.1158573 0.1521821 0.1937134 0.2907309 1.003294 760
## p[63] 0.1588179 0.05854104 0.06456349 0.1158573 0.1521821 0.1937134 0.2907309 1.003294 760
## p[64] 0.1588179 0.05854104 0.06456349 0.1158573 0.1521821 0.1937134 0.2907309 1.003294 760
## p[65] 0.1588179 0.05854104 0.06456349 0.1158573 0.1521821 0.1937134 0.2907309 1.003294 760
## p[66] 0.1588179 0.05854104 0.06456349 0.1158573 0.1521821 0.1937134 0.2907309 1.003294 760
## p[67] 0.1588179 0.05854104 0.06456349 0.1158573 0.1521821 0.1937134 0.2907309 1.003294 760
## p[68] 0.1588179 0.05854104 0.06456349 0.1158573 0.1521821 0.1937134 0.2907309 1.003294 760
## p[69] 0.1588179 0.05854104 0.06456349 0.1158573 0.1521821 0.1937134 0.2907309 1.003294 760
## p[70] 0.1588179 0.05854104 0.06456349 0.1158573 0.1521821 0.1937134 0.2907309 1.003294 760
## p[71] 0.1588179 0.05854104 0.06456349 0.1158573 0.1521821 0.1937134 0.2907309 1.003294 760
## p[72] 0.1588179 0.05854104 0.06456349 0.1158573 0.1521821 0.1937134 0.2907309 1.003294 760
## p[73] 0.1588179 0.05854104 0.06456349 0.1158573 0.1521821 0.1937134 0.2907309 1.003294 760
## p[74] 0.1588179 0.05854104 0.06456349 0.1158573 0.1521821 0.1937134 0.2907309 1.003294 760
## p[75] 0.1588179 0.05854104 0.06456349 0.1158573 0.1521821 0.1937134 0.2907309 1.003294 760
## p[76] 0.1588179 0.05854104 0.06456349 0.1158573 0.1521821 0.1937134 0.2907309 1.003294 760
## p[77] 0.1588179 0.05854104 0.06456349 0.1158573 0.1521821 0.1937134 0.2907309 1.003294 760
## p[78] 0.1588179 0.05854104 0.06456349 0.1158573 0.1521821 0.1937134 0.2907309 1.003294 760
## p[79] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[80] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[81] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[82] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[83] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[84] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[85] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[86] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[87] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[88] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[89] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[90] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[91] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[92] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[93] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[94] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[95] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[96] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[97] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[98] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[99] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[100] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[101] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[102] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[103] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[104] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[105] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[106] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[107] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[108] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[109] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[110] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[111] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[112] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[113] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[114] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[115] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[116] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[117] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[118] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[119] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[120] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[121] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[122] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[123] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[124] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[125] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[126] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[127] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[128] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[129] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[130] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[131] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[132] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[133] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[134] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[135] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[136] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[137] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[138] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[139] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[140] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[141] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[142] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[143] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[144] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[145] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[146] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[147] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[148] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[149] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[150] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[151] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[152] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[153] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[154] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[155] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[156] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[157] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[158] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[159] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[160] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[161] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[162] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[163] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[164] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[165] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[166] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[167] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[168] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[169] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[170] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[171] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[172] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[173] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[174] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[175] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[176] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[177] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[178] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[179] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[180] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[181] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[182] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[183] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[184] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[185] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[186] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[187] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[188] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[189] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[190] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[191] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[192] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[193] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[194] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## p[195] 0.3297387 0.07605704 0.19078940 0.2747733 0.3273642 0.3824401 0.4849671 1.003663 710
## prob.psi.greater.50 0.8124444 0.39040023 0.00000000 1.0000000 1.0000000 1.0000000 1.0000000 1.002017 1500
## psi[1] 0.6105124 0.12288945 0.39337198 0.5220180 0.6027413 0.6877665 0.8746755 1.002649 1100
## psi[2] 0.6105124 0.12288945 0.39337198 0.5220180 0.6027413 0.6877665 0.8746755 1.002649 1100
## psi[3] 0.6105124 0.12288945 0.39337198 0.5220180 0.6027413 0.6877665 0.8746755 1.002649 1100
## psi[4] 0.6105124 0.12288945 0.39337198 0.5220180 0.6027413 0.6877665 0.8746755 1.002649 1100
## psi[5] 0.6105124 0.12288945 0.39337198 0.5220180 0.6027413 0.6877665 0.8746755 1.002649 1100
## psi[6] 0.6105124 0.12288945 0.39337198 0.5220180 0.6027413 0.6877665 0.8746755 1.002649 1100
## psi[7] 0.6105124 0.12288945 0.39337198 0.5220180 0.6027413 0.6877665 0.8746755 1.002649 1100
## psi[8] 0.6105124 0.12288945 0.39337198 0.5220180 0.6027413 0.6877665 0.8746755 1.002649 1100
## psi[9] 0.6105124 0.12288945 0.39337198 0.5220180 0.6027413 0.6877665 0.8746755 1.002649 1100
## psi[10] 0.6105124 0.12288945 0.39337198 0.5220180 0.6027413 0.6877665 0.8746755 1.002649 1100
## psi[11] 0.6105124 0.12288945 0.39337198 0.5220180 0.6027413 0.6877665 0.8746755 1.002649 1100
## psi[12] 0.6105124 0.12288945 0.39337198 0.5220180 0.6027413 0.6877665 0.8746755 1.002649 1100
## psi[13] 0.6105124 0.12288945 0.39337198 0.5220180 0.6027413 0.6877665 0.8746755 1.002649 1100
## psi[14] 0.6105124 0.12288945 0.39337198 0.5220180 0.6027413 0.6877665 0.8746755 1.002649 1100
## psi[15] 0.6105124 0.12288945 0.39337198 0.5220180 0.6027413 0.6877665 0.8746755 1.002649 1100
## psi[16] 0.6105124 0.12288945 0.39337198 0.5220180 0.6027413 0.6877665 0.8746755 1.002649 1100
## psi[17] 0.6105124 0.12288945 0.39337198 0.5220180 0.6027413 0.6877665 0.8746755 1.002649 1100
## psi[18] 0.6105124 0.12288945 0.39337198 0.5220180 0.6027413 0.6877665 0.8746755 1.002649 1100
## psi[19] 0.6105124 0.12288945 0.39337198 0.5220180 0.6027413 0.6877665 0.8746755 1.002649 1100
## psi[20] 0.6105124 0.12288945 0.39337198 0.5220180 0.6027413 0.6877665 0.8746755 1.002649 1100
## psi[21] 0.6105124 0.12288945 0.39337198 0.5220180 0.6027413 0.6877665 0.8746755 1.002649 1100
## psi[22] 0.6105124 0.12288945 0.39337198 0.5220180 0.6027413 0.6877665 0.8746755 1.002649 1100
## psi[23] 0.6105124 0.12288945 0.39337198 0.5220180 0.6027413 0.6877665 0.8746755 1.002649 1100
## psi[24] 0.6105124 0.12288945 0.39337198 0.5220180 0.6027413 0.6877665 0.8746755 1.002649 1100
## psi[25] 0.6105124 0.12288945 0.39337198 0.5220180 0.6027413 0.6877665 0.8746755 1.002649 1100
## psi[26] 0.6105124 0.12288945 0.39337198 0.5220180 0.6027413 0.6877665 0.8746755 1.002649 1100
## psi[27] 0.6105124 0.12288945 0.39337198 0.5220180 0.6027413 0.6877665 0.8746755 1.002649 1100
## psi[28] 0.6105124 0.12288945 0.39337198 0.5220180 0.6027413 0.6877665 0.8746755 1.002649 1100
## psi[29] 0.6105124 0.12288945 0.39337198 0.5220180 0.6027413 0.6877665 0.8746755 1.002649 1100
## psi[30] 0.6105124 0.12288945 0.39337198 0.5220180 0.6027413 0.6877665 0.8746755 1.002649 1100
## psi[31] 0.6105124 0.12288945 0.39337198 0.5220180 0.6027413 0.6877665 0.8746755 1.002649 1100
## psi[32] 0.6105124 0.12288945 0.39337198 0.5220180 0.6027413 0.6877665 0.8746755 1.002649 1100
## psi[33] 0.6105124 0.12288945 0.39337198 0.5220180 0.6027413 0.6877665 0.8746755 1.002649 1100
## psi[34] 0.6105124 0.12288945 0.39337198 0.5220180 0.6027413 0.6877665 0.8746755 1.002649 1100
## psi[35] 0.6105124 0.12288945 0.39337198 0.5220180 0.6027413 0.6877665 0.8746755 1.002649 1100
## psi[36] 0.6105124 0.12288945 0.39337198 0.5220180 0.6027413 0.6877665 0.8746755 1.002649 1100
## psi[37] 0.6105124 0.12288945 0.39337198 0.5220180 0.6027413 0.6877665 0.8746755 1.002649 1100
## psi[38] 0.6105124 0.12288945 0.39337198 0.5220180 0.6027413 0.6877665 0.8746755 1.002649 1100
## psi[39] 0.6105124 0.12288945 0.39337198 0.5220180 0.6027413 0.6877665 0.8746755 1.002649 1100
## z[1] 1.0000000 0.00000000 1.00000000 1.0000000 1.0000000 1.0000000 1.0000000 1.000000 1
## z[2] 1.0000000 0.00000000 1.00000000 1.0000000 1.0000000 1.0000000 1.0000000 1.000000 1
## z[3] 1.0000000 0.00000000 1.00000000 1.0000000 1.0000000 1.0000000 1.0000000 1.000000 1
## z[4] 1.0000000 0.00000000 1.00000000 1.0000000 1.0000000 1.0000000 1.0000000 1.000000 1
## z[5] 1.0000000 0.00000000 1.00000000 1.0000000 1.0000000 1.0000000 1.0000000 1.000000 1
## z[6] 1.0000000 0.00000000 1.00000000 1.0000000 1.0000000 1.0000000 1.0000000 1.000000 1
## z[7] 1.0000000 0.00000000 1.00000000 1.0000000 1.0000000 1.0000000 1.0000000 1.000000 1
## z[8] 1.0000000 0.00000000 1.00000000 1.0000000 1.0000000 1.0000000 1.0000000 1.000000 1
## z[9] 1.0000000 0.00000000 1.00000000 1.0000000 1.0000000 1.0000000 1.0000000 1.000000 1
## z[10] 1.0000000 0.00000000 1.00000000 1.0000000 1.0000000 1.0000000 1.0000000 1.000000 1
## z[11] 1.0000000 0.00000000 1.00000000 1.0000000 1.0000000 1.0000000 1.0000000 1.000000 1
## z[12] 1.0000000 0.00000000 1.00000000 1.0000000 1.0000000 1.0000000 1.0000000 1.000000 1
## z[13] 1.0000000 0.00000000 1.00000000 1.0000000 1.0000000 1.0000000 1.0000000 1.000000 1
## z[14] 1.0000000 0.00000000 1.00000000 1.0000000 1.0000000 1.0000000 1.0000000 1.000000 1
## z[15] 0.2920000 0.45473284 0.00000000 0.0000000 0.0000000 1.0000000 1.0000000 1.001036 4500
## z[16] 0.2831111 0.45056001 0.00000000 0.0000000 0.0000000 1.0000000 1.0000000 1.001507 2400
## z[17] 0.2826667 0.45034575 0.00000000 0.0000000 0.0000000 1.0000000 1.0000000 1.000829 4500
## z[18] 0.2928889 0.45513846 0.00000000 0.0000000 0.0000000 1.0000000 1.0000000 1.000871 4500
## z[19] 0.2866667 0.45225473 0.00000000 0.0000000 0.0000000 1.0000000 1.0000000 1.001136 4200
## z[20] 0.2875556 0.45267305 0.00000000 0.0000000 0.0000000 1.0000000 1.0000000 1.002740 970
## z[21] 0.2757778 0.44695502 0.00000000 0.0000000 0.0000000 1.0000000 1.0000000 1.001396 2700
## z[22] 0.2795556 0.44883072 0.00000000 0.0000000 0.0000000 1.0000000 1.0000000 1.000867 4500
## z[23] 0.2744444 0.44628349 0.00000000 0.0000000 0.0000000 1.0000000 1.0000000 1.001368 2800
## z[24] 0.2800000 0.44904878 0.00000000 0.0000000 0.0000000 1.0000000 1.0000000 1.002707 980
## z[25] 0.2782222 0.44817324 0.00000000 0.0000000 0.0000000 1.0000000 1.0000000 1.001240 3500
## z[26] 0.2891111 0.45340000 0.00000000 0.0000000 0.0000000 1.0000000 1.0000000 1.000672 4500
## z[27] 0.2735556 0.44583302 0.00000000 0.0000000 0.0000000 1.0000000 1.0000000 1.001391 2800
## z[28] 0.2771111 0.44762157 0.00000000 0.0000000 0.0000000 1.0000000 1.0000000 1.001779 1800
## z[29] 0.2871111 0.45246416 0.00000000 0.0000000 0.0000000 1.0000000 1.0000000 1.001223 3600
## z[30] 0.2733333 0.44572006 0.00000000 0.0000000 0.0000000 1.0000000 1.0000000 1.000992 4500
## z[31] 0.2937778 0.45554200 0.00000000 0.0000000 0.0000000 1.0000000 1.0000000 1.000907 4500
## z[32] 0.2795556 0.44883072 0.00000000 0.0000000 0.0000000 1.0000000 1.0000000 1.001427 2600
## z[33] 0.2777778 0.44795298 0.00000000 0.0000000 0.0000000 1.0000000 1.0000000 1.001449 2600
## z[34] 0.2806667 0.44937486 0.00000000 0.0000000 0.0000000 1.0000000 1.0000000 1.002059 1400
## z[35] 0.2755556 0.44684344 0.00000000 0.0000000 0.0000000 1.0000000 1.0000000 1.002701 990
## z[36] 1.0000000 0.00000000 1.00000000 1.0000000 1.0000000 1.0000000 1.0000000 1.000000 1
## z[37] 1.0000000 0.00000000 1.00000000 1.0000000 1.0000000 1.0000000 1.0000000 1.000000 1
## z[38] 1.0000000 0.00000000 1.00000000 1.0000000 1.0000000 1.0000000 1.0000000 1.000000 1
## z[39] 1.0000000 0.00000000 1.00000000 1.0000000 1.0000000 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.735011 1.004481
##
## $beta.psi
## [1] 0.4964369
##
## $deviance
## [1] 130.7801
##
## $occ.sites
## [1] 23.92044
##
## $p
## [1] 0.1588179 0.1588179 0.1588179 0.1588179 0.1588179 0.1588179 0.1588179 0.1588179 0.1588179 0.1588179 0.1588179 0.1588179 0.1588179 0.1588179 0.1588179 0.1588179 0.1588179 0.1588179 0.1588179 0.1588179 0.1588179 0.1588179 0.1588179 0.1588179 0.1588179 0.1588179 0.1588179 0.1588179 0.1588179 0.1588179 0.1588179 0.1588179 0.1588179 0.1588179 0.1588179 0.1588179 0.1588179 0.1588179 0.1588179
## [40] 0.1588179 0.1588179 0.1588179 0.1588179 0.1588179 0.1588179 0.1588179 0.1588179 0.1588179 0.1588179 0.1588179 0.1588179 0.1588179 0.1588179 0.1588179 0.1588179 0.1588179 0.1588179 0.1588179 0.1588179 0.1588179 0.1588179 0.1588179 0.1588179 0.1588179 0.1588179 0.1588179 0.1588179 0.1588179 0.1588179 0.1588179 0.1588179 0.1588179 0.1588179 0.1588179 0.1588179 0.1588179 0.1588179 0.1588179
## [79] 0.3297387 0.3297387 0.3297387 0.3297387 0.3297387 0.3297387 0.3297387 0.3297387 0.3297387 0.3297387 0.3297387 0.3297387 0.3297387 0.3297387 0.3297387 0.3297387 0.3297387 0.3297387 0.3297387 0.3297387 0.3297387 0.3297387 0.3297387 0.3297387 0.3297387 0.3297387 0.3297387 0.3297387 0.3297387 0.3297387 0.3297387 0.3297387 0.3297387 0.3297387 0.3297387 0.3297387 0.3297387 0.3297387 0.3297387
## [118] 0.3297387 0.3297387 0.3297387 0.3297387 0.3297387 0.3297387 0.3297387 0.3297387 0.3297387 0.3297387 0.3297387 0.3297387 0.3297387 0.3297387 0.3297387 0.3297387 0.3297387 0.3297387 0.3297387 0.3297387 0.3297387 0.3297387 0.3297387 0.3297387 0.3297387 0.3297387 0.3297387 0.3297387 0.3297387 0.3297387 0.3297387 0.3297387 0.3297387 0.3297387 0.3297387 0.3297387 0.3297387 0.3297387 0.3297387
## [157] 0.3297387 0.3297387 0.3297387 0.3297387 0.3297387 0.3297387 0.3297387 0.3297387 0.3297387 0.3297387 0.3297387 0.3297387 0.3297387 0.3297387 0.3297387 0.3297387 0.3297387 0.3297387 0.3297387 0.3297387 0.3297387 0.3297387 0.3297387 0.3297387 0.3297387 0.3297387 0.3297387 0.3297387 0.3297387 0.3297387 0.3297387 0.3297387 0.3297387 0.3297387 0.3297387 0.3297387 0.3297387 0.3297387 0.3297387
##
## $prob.psi.greater.50
## [1] 0.8124444
##
## $psi
## [1] 0.6105124 0.6105124 0.6105124 0.6105124 0.6105124 0.6105124 0.6105124 0.6105124 0.6105124 0.6105124 0.6105124 0.6105124 0.6105124 0.6105124 0.6105124 0.6105124 0.6105124 0.6105124 0.6105124 0.6105124 0.6105124 0.6105124 0.6105124 0.6105124 0.6105124 0.6105124 0.6105124 0.6105124 0.6105124 0.6105124 0.6105124 0.6105124 0.6105124 0.6105124 0.6105124 0.6105124 0.6105124 0.6105124 0.6105124
##
## $z
## [1] 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 0.2920000 0.2831111 0.2826667 0.2928889 0.2866667 0.2875556 0.2757778 0.2795556 0.2744444 0.2800000 0.2782222 0.2891111 0.2735556 0.2771111 0.2871111 0.2733333 0.2937778 0.2795556 0.2777778 0.2806667 0.2755556 1.0000000 1.0000000 1.0000000 1.0000000
results$BUGSoutput$mean$psi
## [1] 0.6105124 0.6105124 0.6105124 0.6105124 0.6105124 0.6105124 0.6105124 0.6105124 0.6105124 0.6105124 0.6105124 0.6105124 0.6105124 0.6105124 0.6105124 0.6105124 0.6105124 0.6105124 0.6105124 0.6105124 0.6105124 0.6105124 0.6105124 0.6105124 0.6105124 0.6105124 0.6105124 0.6105124 0.6105124 0.6105124 0.6105124 0.6105124 0.6105124 0.6105124 0.6105124 0.6105124 0.6105124 0.6105124 0.6105124
# the results$BUGSoutput$sims.array is a 3-d object [iterations, chains, variables]
dim(results$BUGSoutput$sims.array)
## [1] 1500 3 279
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.228669 0.3551734 0.61608755 122.9064 21 0.2264145 0.2264145 0.2264145 0.2264145 0.2264145
## [2,] -1.500943 0.6382737 0.76369368 130.7159 24 0.1822850 0.1822850 0.1822850 0.1822850 0.1822850
## [3,] -1.448151 0.2129412 1.20810563 153.1177 32 0.1902862 0.1902862 0.1902862 0.1902862 0.1902862
## [4,] -1.694131 0.6413470 1.29780410 141.3284 28 0.1552333 0.1552333 0.1552333 0.1552333 0.1552333
## [5,] -1.813665 0.5675678 0.09957935 144.2208 29 0.1401958 0.1401958 0.1401958 0.1401958 0.1401958
results$BUGSoutput$sims.array[1:5,1,"psi[1]", drop=FALSE]
## , , psi[1]
##
## [,1]
## [1,] 0.6493282
## [2,] 0.6821551
## [3,] 0.7699636
## [4,] 0.7854652
## [5,] 0.5248743
# 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 279
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.077035 0.6024707 0.5978249 122.4838 21 0.25406761 0.25406761 0.25406761 0.25406761 0.25406761
## [2,] -2.117245 1.6087376 0.5685131 131.7492 24 0.10743195 0.10743195 0.10743195 0.10743195 0.10743195
## [3,] -2.538809 1.6877522 0.4624011 130.7495 23 0.07318191 0.07318191 0.07318191 0.07318191 0.07318191
## [4,] -1.681517 0.5642576 0.9174025 146.3643 30 0.15689465 0.15689465 0.15689465 0.15689465 0.15689465
## [5,] -2.212409 1.2695266 1.4197743 146.0190 30 0.09864172 0.09864172 0.09864172 0.09864172 0.09864172
results$BUGSoutput$sims.matrix[1:5,"psi[1]", drop=FALSE]
## psi[1]
## [1,] 0.6451585
## [2,] 0.6384200
## [3,] 0.6135836
## [4,] 0.7145126
## [5,] 0.8053030
# make a posterior density plot
plotdata <- data.frame(parm=results$BUGSoutput$sims.matrix[,"psi[1]"], stringsAsFactors=FALSE)
head(plotdata)
## parm
## 1 0.6451585
## 2 0.6384200
## 3 0.6135836
## 4 0.7145126
## 5 0.8053030
## 6 0.6510562
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
## Warning: The dot-dot notation (`..density..`) was deprecated in ggplot2 3.4.0.
## ℹ Please use `after_stat(density)` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.
## `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.6493282 0.5986172 0.5229432 1
## 2 0.6821551 0.5987142 0.5624629 2
## 3 0.7699636 0.7240780 0.5749332 3
## 4 0.7854652 0.7195014 0.4860518 4
## 5 0.5248743 0.5421892 0.4504520 5
## 6 0.6721033 0.6248019 0.5417386 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.6493282
## 2 2 psi.1 0.6821551
## 3 3 psi.1 0.7699636
## 4 4 psi.1 0.7854652
## 5 5 psi.1 0.5248743
## 6 6 psi.1 0.6721033
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.014 -0.028 0.008 -0.007 0.018 0.008 -0.005 -0.001 0.007 -0.014 -0.010 0.028 -0.014 0.012 -0.007 -0.006 0.002 -0.035 -0.016 -0.015 0.017 -0.024 0.009 -0.017 -0.014 -0.002 -0.016 0.004 0.002 -0.014 0.016 0.007 0.024 -0.017 -0.006 -0.011
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)