I am estimating the number of adult salmon salmon which spawned in a stream
from carcass tagging data. I estimate the carcasses at large using a Petersen
estimate and carcass persistence using an exponetial decay function. Weekly
estimate is the carcasses at large minus carcasses at large that persist into
the next week, plus the number of marked carcasses removed from the
population. Assuming independence of the periods, and sampling with
replacement I use the following code to estimatethe weekly population and
sum them together to get the total population estimate (Popsum). The code
works but with high autocorrelation.
Model:{
#Vague Priors
for( i in 1 : N ) {
U[i]~dunif(T[i],100000)
q[i] ~dbeta(1,1)
}
S ~dbeta(1,1)
r ~dbeta(1,1)
tau2 ~ dgamma(0.001,0.001)
#Estimate survival (S) from exponetial decay function
for (i in 1:NN){
log_mu[i] <- x[i]*log((1-r)*S)+log(r/(1-r))
ln_mM[i]~dnorm(log_mu[i],tau2)
}
#Use Binomial to estimate weekly abundance U[i]
for( i in 1 : N ) {
R[i] ~dbin(q[i],T[i])
C[i] ~dbin(q[i],U[i])
#Estimate X[i] carcasses that persisted from previous weeks
X[i] <- U[i] * S
#Estimate Weekly Escapement
Pop[i] <- U[i] - X[i] + m[i]
}
#Sum of weekly estimates
Popsum <- sum(Pop[])
}
# data
list(N=18,NN=54,x=c
(1,2,3,1,4,1,2,1,2,3,1,2,3,4,1,2,3,1,2,3,1,2,3,4,1,2,3,4,1,2,3,5,1,2,3,1,2,3,1,2,
3,4,5,1,2,3,4,1,2,3,4,1,2,3), ln_mM=c(-1.38629,-1.79176,-2.48491,-1.09861,-
2.19722,-0.74444,-2.30259,-1.51583,-2.10788,-4.67283,-1.33685,-2.74084,-
3.58814,-4.68675,-1.34249,-3.2884,-4.89784,-1.3273,-2.25406,-3.58906,-
1.00764,-3.16231,-2.84385,-4.5486,-1.46546,-2.41644,-3.30374,-5.24965,-
1.32176,-1.97981,-2.90476,-5.34711,-1.25063,-1.89997,-3.2581,-1.2209,-
2.44957,-2.98856,-1.4929,-2.13726,-3.79549,-4.08317,-5.18178,-0.87184,-
1.81011,-2.75457,-4.00733,-1.76766,-2.46081,-2.00882,-3.30811,-0.87925,-
3.27714,-3.97029))
T[] C[] R[] m[]
12 14 3 6
9 25 3 4
80 118 38 46
214 449 47 75
217 573 57 79
134 491 35 41
181 456 48 72
189 671 69 90
381 956 88 138
420 1222 112 193
234 1058 67 111
139 1388 41 60
178 1311 40 68
110 622 46 73
82 434 14 35
53 218 22 25
47 167 20 22
38 151 9 10
END
#initials
list(U=c
(2000,2000,2000,2000,2000,2000,2000,2000,2000,2000,2000,2000,2000,2000,2
000,2000,2000,2000), q=c
(0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5),
S=0.5, r=0.5,tau2=1))
However, I would rather put a prior on the total population estimate (Popsum)
and assume the weekly proportion (p) follows a multinomial distribution based
on Pop[1:18] and use Dirichlet priors for p. However, N (Popsum) for a
multinomial in WinBUGS must be fixed. Therefore I tried adapting the code from
the users manual "learning about the parameters of a Dirichlet"and the
archives without success. When running the second model (code below), I
get the error message - "value of order binomial c[1] must be greater than
zero" despite changes in priors and initial values. On the surface this seems
relatively straightforward but I am not getting it. I would appreciate any
suggestions to this code? Thanks in advance for your time!
Dan Rawding
Model:{
#Vague Priors
for( i in 1 : N ) {
#U[i]~dunif(T[i],100000)
q[i] ~dbeta(1,1)
}
S ~dbeta(1,1)
r ~dbeta(1,1)
tau2 ~ dgamma(0.001,0.001)
#Estimate survival (S) from exponetial decay function
for (i in 1:NN){
log_mu[i] <- x[i]*log((1-r)*S)+log(r/(1-r))
ln_mM[i]~dnorm(log_mu[i],tau2)
}
#Use Binomial to estimate weekly abundance U[i]
for( i in 1 : N ) {
R[i] ~dbin(q[i],T[i])
C[i] ~dbin(q[i],U[i])
#Estimate X[i] carcasses that persisted from previous weeks
X[i] <- U[i] * S
#Estimate Weekly Escapement
Pop[i] <- U[i] - X[i] + m[i]
}
#Sum of weekly estimates
#Popsum <- sum(Pop[])
#Code to change to put priors on a single population
min <- sum(T[]) # min = all marked individuals
Popsum~dunif(min,50000) # prior for total population
for ( i in 1:18){
U[i] <- p[i]/Popsum
p[i] ~dgamma(alpha[i],1)
alpha[i] ~dgamma(1, 0.001)
}
}
# data is same as above
#initials
list(Popsum=15000, q=c
(0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5),
S=0.5, r=0.5,tau2=1))
-------------------------------------------------------------------
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
|
