Dear all,
I have a hierarchical time series model (see below) in which the
temporal evolution of an environmental variable (denoted by GPS) is
related to a climatic one (CLIM); this equation lacks the Markov
property, but it is assumed to be measured with sampling error
(GPS.tau2). A second set of equations describes the evolution of the
size of a bird population (POPSIZE) according to a state-space model,
plus an effect of the environmental variable GPS. The problem is that
there are missing observations in both the POPSIZE and GPS variables.
My goal is to use the "cut()" function to prevent the model from fitting
the parameters of the process equations to the missing data, as
suggested by the winbugs manual, but I'm not sure of how to implement
this in the correct way. Below is the code using the cut function, but I
want to check if this is the correct way of using it to this end. Does
anyone have experience with this function in a similar setting?
Many thanks in advance
Pablo
model; {
# Observation models
for(i in 1:27){
POPSIZE[i+1] ~ dnorm(POPSIZE.state[i+1] , POPSIZE.itau2)
GPS[i+1] ~ dnorm(GPS.exp[i+1], GPS.itau2)
}
# Initial states
POPSIZE.state[1] ~ dnorm(POPSIZE[1], POPSIZE.state.isigma2)
GPS.exp[1] ~ dnorm(GPS[1], GPS.isigma2)
GPS.cut[k, 1] <- cut(GPS.exp[k, 1]) # (ˇˇ)Is this OK?
# Process model
for(i in 1:27){
GPS.exp.mu[i+1] <- rho*CLIM[i+1]
GPS.exp[i+1] ~ dnorm(GPS.exp.mu[i+1], GPS.isigma2)
GPS.cut[i+1] <- cut(GPS.exp[i+1]) # (ˇˇ) Is this OK?
POPSIZE.state.mu[i+1] <- beta0 + b1*POPSIZE.state[i] + alfa*GPS.cut[i+1]
POPSIZE.state[i+1] ~ dnorm(POPSIZE.state.mu[i+1], POPSIZE.state.isigma2)
}
# Parameters and priors
beta0 ~ dnorm(0,1.E-10)
b1 ~ dnorm(0,1.E-10)
alfa ~ dnorm(0,1.E-10)
rho ~ dnorm(0,1.E-10)
POPSIZE.state.isigma2 ~ dalfa(1.E-3,1.E-3)
POPSIZE.state.sigma2 <-1/POPSIZE.state.isigma2
GPS.isigma2 ~ dalfa(1.E-3,1.E-3)
GPS.sigma2 <- 1/GPS.isigma2
POPSIZE.itau2 ~ dalfa(1.E-3,1.E-3)
POPSIZE.tau2 <- 1/POPSIZE.itau2
GPS.itau2 ~ dalfa(1.E-3,1.E-3)
GPS.tau2 <- 1/GPS.itau2
}
# data
list(
CLIM=c(6.108,5.92,5.169,5.237,5.782,5.248,6.029,6.063,5.435,5.75,6.629,5.651,6.648,5.121,5.689,5.22,5.549,4.927,6.62,6.613,6.305,5.191,5.745,6.261,5.977,5.988,6.224,4.626),
GPS=c(8.574,NA,NA,NA,10.015,NA,10.369,NA,9.563,8.884,10.512,10.233,10.595,NA,9.857,8.689,9.2,8.467,NA,NA,NA,9.074,10.246,10.911,10.412,10.733,10.392,NA),
POPSIZE=c(10.662,11.13,11.145,10.909,10.561,11.018,NA,11.116,10.974,10.775,10.575,11.173,10.47,10.596,10.95,10.752,10.739,NA,10.087,10.439,9.814,10.78,10.596,10.296,11.102,10.519,10.317,10.985))
# inits
list(
beta0=14.84,
b1=-0.215,
rho=1.702,
GPS.isigma2=8.6730,
GPS.itau2=7.855,
alfa=-0.1873,
POPSIZE.itau2=18.5359,
POPSIZE.state.isigma2=36.865,
POPSIZE.state=c(10.66,10.809,10.78,11.01,10.69,10.65,10.48,10.76,10.93,10.66,10.47,10.93,10.36,10.91,10.94,11,10.685,11.16,10.105,10.77,10.51,10.92,10.58,10.29,10.84,10.66,10.451,11.13),
GPS.exp=c(9.134,9.638,8.625,8.72,9.464,8.731,9.792,9.827,8.996,9.395,10.578,9.29,10.604,8.568,9.336,8.698,9.139,8.3,10.57,10.561,10.157,8.665,9.415,10.103,9.722,9.742,10.05,7.893),
POPSIZE=c(NA,NA,NA,NA,NA,NA,11.161,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,11.279,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA),
GPS=c(NA,9.638,8.625,8.72,NA,8.731,NA,9.827,NA,NA,NA,NA,NA,8.568,NA,NA,NA,NA,10.57,10.561,10.157,NA,NA,NA,NA,NA,NA,7.893))
__________ Información de ESET NOD32 Antivirus, versión de la base de firmas de virus 4647 (20091129) __________
ESET NOD32 Antivirus ha comprobado este mensaje.
http://www.eset.com
-------------------------------------------------------------------
This list is for discussion of modelling issues and the BUGS software.
For help with crashes and error messages, first mail [log in to unmask]
To mail the BUGS list, mail to [log in to unmask]
Before mailing, please check the archive at www.jiscmail.ac.uk/lists/bugs.html
Please do not mail attachments to the list.
To leave the BUGS list, send LEAVE BUGS to [log in to unmask]
If this fails, mail [log in to unmask], NOT the whole list
|