Supporting Code

WinBUGS Script

model{
for(i in 1:totN){  
   
    # predict survival
    logit(Surv[i])<-(S[year[i],from[i]]+SQ[quad[i],from[i]])+
Sind[year[i],from[i]]*plant[i]+
S1[year[i],from[i]]*popS1[i]+
S2[year[i],from[i]]*popS2[i]+
S3[year[i],from[i]]*popS3[i]
    # convert to probability, preventing numerical overflow
    # and make sure bare ground is 0
    PSurv[i]<-max(0.000000000000,min(0.999999999999,Surv[i]))*
sdummy[from[i]]
    
    # add "self" abundance in neighborhood to conspp. abundance
    popTot[i,1]<-popS1[i]+plant[i]*dummy[from[i],1]
    popTot[i,2]<-popS2[i]+plant[i]*dummy[from[i],2]
    popTot[i,3]<-popS3[i]+plant[i]*dummy[from[i],3]
    
    for(k in 1:3){

      # predict colonization
      logit(Cest[i,k])<-(C[year[i],k]+CQ[quad[i],k])+
C1[year[i],k]*popTot[i,1]+
C2[year[i],k]*popTot[i,2]+C3[year[i],k]*popTot[i,3]
      # convert to probability and prevent numerical overflow
      PC[i,k]<-max(0.000000000000,min(0.999999999999,Cest[i,k]))

     }   # next k  
    
# combine survival and colonization probabilities to get 
#multinomial probabilities
p[i,1]<-PSurv[i]*dummy[1,from[i]]+(1-PSurv[i])*(PC[i,1]*
 (1-PC[i,2])*(1-PC[i,3])+PC[i,1]*PC[i,2]*(1-PC[i,3])/2+
 PC[i,1]*(1-PC[i,2])*PC[i,3]/2+PC[i,1]*PC[i,2]*PC[i,3]/3)
     p[i,2]<-PSurv[i]*dummy[2,from[i]]+(1-PSurv[i])*(PC[i,2]*
 (1-PC[i,1])*(1-PC[i,3])+PC[i,1]*PC[i,2]*1-PC[i,3])/2+
 (1-PC[i,1])*PC[i,2]*PC[i,3]/2+PC[i,1]*PC[i,2]*PC[i,3]/3)
     p[i,3]<-PSurv[i]*dummy[3,from[i]]+(1-PSurv[i])*(PC[i,3]*
 (1-PC[i,1])*(1-PC[i,2])+PC[i,1]*(1-PC[i,2])*PC[i,3]/2+
 (1-PC[i,1])*PC[i,2]*PC[i,3]/2+PC[i,1]*PC[i,2]*PC[i,3]/3)
     p[i,4]<-max(0,min(1,(1-sum(p[i,1:3]))))
    
    # compare predicted and observed
    y[i,1:4]~dmulti(p[i,1:4],1)
} # next i

# PRIORS
# set all surv. and coloniz. parameters for bare ground to 0
S.mu[4]<-0
S.tau[4]<-0
Sind.mu[4]<-0
Sind.tau[4]<-0
S1.mu[4]<-0
S1.tau[4]<-0
S2.mu[4]<-0
S2.tau[4]<-0
S3.mu[4]<-0
S3.tau[4]<-0
C.mu[4]<-0
C.tau[4]<-0
C1.mu[4]<-0
C1.tau[4]<-0
C2.mu[4]<-0
C2.tau[4]<-0
C3.mu[4]<-0
C3.tau[4]<-0
for(iyr in 1:nyrs){
   S[iyr,4]<-0
   Sind[iyr,4]<-0
   S1[iyr,4]<-0
   S2[iyr,4]<-0
   S3[iyr,4]<-0
   C[iyr,4]<-0
   C1[iyr,4]<-0
   C2[iyr,4]<-0
   C3[iyr,4]<-0
}
# quad effects
for(q in 1:nquads){
  SQ[q,4]<-0
  CQ[q,4]<-0
} # next q

for(k in 1:3){
   # update means and variances
   S.mu[k]~dnorm(0,0.001)
   S.tau[k]~dgamma(1,100)
   Sind.mu[k]~dnorm(0,0.0002)
   Sind.tau[k]~dgamma(1,100)
   S1.mu[k]~dnorm(0,0.0002)
   S1.tau[k]~dgamma(1,100)
   S2.mu[k]~dnorm(0,0.0002)
   S2.tau[k]~dgamma(1,100)
   S3.mu[k]~dnorm(0,0.0002)
   S3.tau[k]~dgamma(1,100)
   C.mu[k]~dnorm(0,0.001)
   C.tau[k]~dgamma(1,100)
   C1.mu[k]~dnorm(0,0.0002)
   C1.tau[k]~dgamma(1,100)
   C2.mu[k]~dnorm(0,0.0002)
   C2.tau[k]~dgamma(1,100)
   C3.mu[k]~dnorm(0,0.0002)
   C3.tau[k]~dgamma(1,100)
   for(iyr in 1:nyrs){
      # update year-specific effects
      S[iyr,k]~dnorm(S.mu[k],S.tau[k])
      Sind[iyr,k]~dnorm(Sind.mu[k],Sind.tau[k])
      S1[iyr,k]~dnorm(S1.mu[k],S1.tau[k])
      S2[iyr,k]~dnorm(S2.mu[k],S2.tau[k])
      S3[iyr,k]~dnorm(S3.mu[k],S3.tau[k])
      C[iyr,k]~dnorm(C.mu[k],C.tau[k])
      C1[iyr,k]~dnorm(C1.mu[k],C1.tau[k])
      C2[iyr,k]~dnorm(C2.mu[k],C2.tau[k])
      C3[iyr,k]~dnorm(C3.mu[k],C3.tau[k])
   } # next j
   # update quadrat effects   
   for(q in 1:(nquads-1)){
      SQ[q,k]~dnorm(0,0.001)
      CQ[q,k]~dnorm(0,0.001)
   } # next q   
   SQ[4,k]<-0
   CQ[4,k]<-0
} # next k

} # end model