I am trying to model the number of seedlings in meter-square plots in a forest. Let old.true[i,j] and new.true[i,j] be the true numbers of old and new seedlings in plot i and year j. Each year a botanist goes to each plot looking for seedlings and counts old.count[i,j] and new.count[i,j]. The model is (not strict BUGS notation)old.count[i,j] ~ binomial ( old.count.p, old.true[i,j] )
new.count[i,j] ~ binomial ( new.count.p, new.true[i,j] )Each year some seedlings survive and others die. The numbers of survivors are modelled asold.surv[i,j] ~ binomial ( old.surv.p, old.true[i,j] )
new.surv[i,j] ~ binomial ( new.surv.p, new.true[i,j] )

 

Those that survive become the old seedlings in the subsequent year.

old.true[i,j+1] <- old.surv[i,j] + new.surv[i,j]

The true number of new seedlings each year is random

new.true[i,j] ~ categorical with values from 0 to 19.

The BUGS problem arises because the second parameter of these binomial distributions may be 0.

Here is a hack to work around the problem. I will first generate new.true[i,j], then new.surv[i,j].# generate new.true[i,j] #

new.index[i,j] ~ dcat(prior.new[]); # choose a number 1-20
new.true[i,j] <- new.index[i,j]-1;
# when new.true[i,j]=0 we can’t use it as the n in a binomial.
# instead, we create a phony binomial with n=1 and p=0.
# New.useless[i,j,] is a vector containing one 1 and the rest 0’s.
# It is used to generate new.n[i,j], which is the n in the phony binomial.
# Because new.n is generated using dcat, BUGS recognizes that it
# is positive.

 

 

new.useless[i,j,1] <- equals(0,new.true[i,j]) + equals(1,new.true[i,j]);

for ( k in 2:(max.new-1) ) {
new.useless[i,j,k] <- equals(k,new.true[i,j]);

}

new.n[i,j] ~ dcat( new.useless[i,j,] );

 

# new.surv.p[i,j] is the phony p.

new.surv.p[i,j]
= 0,
= new.surv.prob,

if new.true[i,j] == 0
if new.true[i,j] > 0

new.surv.p[i,j] <- new.surv.prob * ( 1 – equals( 0, new.true[i,j] ));

# Now we’re ready to generate

new.surv[i,j] ~ dbin( new.surv.p[i,j], new.n[i,j] )

model seedling;

const
n.plots = 60,
n.years = 5,
max.new = 20,
max.old = 20;

var
prior.old[max.old],
prior.new[max.new],

# unknowns
old.find.prob,
new.find.prob,

old.surv.prob,
new.surv.prob,

old.true[n.plots,n.years],
old.index[n.plots],
new.true[n.plots,n.years],
new.index[n.plots,n.years],

new.surv[n.plots,n.years-1],
old.surv[n.plots,n.years-1],

# data
new.count[n.plots,n.years],
old.count[n.plots,n.years],

# prior for no. of old seedlings in year 1
# prior for no. of new seedlings each year

# probabilities of finding seedlings

# survival probabilities

# number of old seedlings in plot i, year j
# index for no. old seedlings in plot i, year 1
# number of new seedlings in plot i, year j
# index for no. new seedlings in plot i, year j

# number of new sdlngs that survive to year j+1
# number of old sdlngs that survive to year j+1

# number of new sdlngs counted in plot i, year j
# number of old sdlngs counted in plot i, year j

# some phony binomial parameters to handle the case when the number of
# trials is 0
new.count.p[n.plots,n.years],
new.surv.p[n.plots,n.years-1],
new.n[n.plots,n.years],
new.useless[n.plots,n.years,max.new-1],
old.count.p[n.plots,n.years],
old.surv.p[n.plots,n.years-1],
old.n[n.plots,n.years],
old.useless[n.plots,n.years,max.old-1],

# logical constraints on n and o imposed by the model
min.new.surv[n.plots,n.years-1],
max.new.surv[n.plots,n.years-1],
min.new.index[n.plots,n.years],
max.new.index[n.plots,n.years],
min.old.surv[n.plots,n.years-1],
max.old.surv[n.plots,n.years-1],
min.old.index[n.plots],
max.old.index[n.plots];

data in “bugs.dat”;

inits in “bugs.init”;

{
old.find.prob ~ dunif(0,1);
new.find.prob ~ dunif(0,1);
old.surv.prob ~ dunif(0,1);
new.surv.prob ~ dunif(0,1);

##################
# the number of new seedlings in plot i, year j
##################
for ( i in 1:max.new ) {
prior.new[i] <- 1/max.new;
}

for ( i in 1:n.plots ) {
for ( j in 1:n.years-1 ) {
min.new.index[i,j] <- max( new.count[i,j], new.surv[i,j] ) + 1;
max.new.index[i,j] <- max.new;
new.index[i,j] ~ dcat(prior.new[]) I(min.new.index[i,j], max.new.index[i,j]);
new.true[i,j] <- new.index[i,j]-1;
}

min.new.index[i,n.years] <- new.count[i,n.years] + 1;
max.new.index[i,n.years] <- max.new;
new.index[i,n.years] ~ dcat(prior.new[]) I(min.new.index[i,n.years], max.new.index[i,n.years]);
new.true[i,n.years] <- new.index[i,n.years]-1;

# this bit handles binomial distributions when new.true[i,j]
equals 0
for ( j in 1:n.years ) {
new.useless[i,j,1] <- equals(0,new.true[i,j]) + equals(1,new.true[i,j]);
for ( k in 2:(max.new-1) ) {
new.useless[i,j,k] <- equals(k,new.true[i,j]);
}
new.n[i,j] ~ dcat( new.useless[i,j,] );
}

}

##################
# the number of old seedlings in plot i, year 1
##################
for ( i in 1:max.old ) {
prior.old[i] <- 1/max.old;
}
for ( i in 1:n.plots ) {
min.old.index[i] <- max( old.count[i,1], old.surv[i,1] ) + 1;
max.old.index[i] <- max.old;
old.index[i] ~ dcat(prior.old[]) I(min.old.index[i],
max.old.index[i] );
old.true[i,1] <- old.index[i]-1;
for ( j in 2:(n.years) ) {
old.true[i,j] <- old.surv[i,j-1]+new.surv[i,j-1];
}

# this bit handles binomial distributions when old.true[i,j]
equals 0
for ( j in 1:n.years ) {
old.useless[i,j,1] <- equals(0,old.true[i,j]) +
equals(1,old.true[i,j]);
for ( k in 2:(max.old-1) ) {
old.useless[i,j,k] <- equals(k,old.true[i,j]);
}
old.n[i,j] ~ dcat( old.useless[i,j,] );
}
}

######################
# the number of new seedlings that survive in plot i, year 1
######################
for ( i in 1:n.plots ) {
for ( j in 1:(n.years-1) ) {
new.surv.p[i,j] <- new.surv.prob * ( 1 – equals( 0, new.true[i,j]));
new.surv[i,j] ~ dbin( new.surv.p[i,j], new.n[i,j] )
I( min.new.surv[i,j], max.new.surv[i,j] );
}
}

######################
# the number of old seedlings that survive in plot i, year 1
######################
for ( i in 1:n.plots ) {
for ( j in 1:(n.years-1) ) {
old.surv.p[i,j] <- old.surv.prob * ( 1 – equals( 0, old.true[i,j]));
old.surv[i,j] ~ dbin( old.surv.p[i,j], old.n[i,j] )
I( min.old.surv[i,j], max.old.surv[i,j] );
}
}

######################
# the number of new seedlings counted in plot i, year 1
######################
for ( i in 1:n.plots ) {
for ( j in 1:n.years ) {
new.count.p[i,j] <- new.find.prob * ( 1 – equals( 0, new.true[i,j]));
new.count[i,j] ~ dbin( new.count.p[i,j], new.n[i,j] );
}
}

######################
# the number of old seedlings counted in plot i, year 1
######################
for ( i in 1:n.plots ) {
for ( j in 1:n.years ) {
old.count.p[i,j] <- old.find.prob * ( 1 – equals( 0, old.true[i,j]));
old.count[i,j] ~ dbin( old.count.p[i,j], old.n[i,j] );
}
}

# logical constraints imposed by the model
for ( i in 1:n.plots ) {
for ( j in 1:(n.years-2) ) {
min.new.surv[i,j] <- max( old.count[i,j+1] – old.surv[i,j],
max( 0, old.surv[i,j+1]-old.surv[i,j]));
max.new.surv[i,j] <- max.old – old.surv[i,j];

min.old.surv[i,j] <- max( old.count[i,j+1] – new.surv[i,j],
max( 0, old.surv[i,j+1]-new.surv[i,j]));
max.old.surv[i,j] <- max.old – new.surv[i,j];
}
min.new.surv[i,n.years-1] <- max( 0, old.count[i,n.years] –
old.surv[i,n.years-1] );
max.new.surv[i,n.years-1] <- max.old – old.surv[i,n.years-1];
min.old.surv[i,n.years-1] <- max( 0, old.count[i,n.years] –
new.surv[i,n.years-1] );
max.old.surv[i,n.years-1] <- max.old – new.surv[i,n.years-1];
}