model
{

	###### Hierarchical GAM model with additional random year-effects
	####  original priors with small improvements to parameterization to speed convergence
	####  uses original Crainiceanu et al basis function
	#### changed parameter names to better match model description
	#### smooth = GAM-based smooth


	for( k in 1 : ncounts )
	{
		Eloglambda[k] <- obs[strat[k],obser[k]] + eta*firstyr[k] + strata[strat[k]] + smooth[strat[k],year[k]]
	 	loglambda[k] ~ dnorm(Eloglambda[k], taunoise)
	 	log(lambda[k]) <- loglambda[k]
		count[k] ~ dpois(lambda[k])


	}



	### fixed effect priors

  taunoise ~ dgamma(0.001,0.001)
	sdnoise <- 1 / pow(taunoise, 0.5)
	eta ~ dnorm( 0.0,0.01)
	STRATA ~ dnorm( 0.0,0.01)
	taustrata ~ dgamma(0.001,0.0001)
	sdstrata <- 1/pow(taustrata,0.5)

  sdobs <- 1/pow(tauobs, 0.5)
  tauobs ~ dgamma(0.001,0.0001)

	#### stratum-level effects  ######
	for( i in 1 : nstrata )
	{
		#### observer effects  ######

		for( o in 1 : nobservers[i] )
		{

			obs[i,o] ~ dnorm(0.0, tauobs)
		}

		#### end observer effects  ######


		### stratum-level priors


		strata[i] ~ dnorm(STRATA,taustrata)

	}# end s strata loop and stratum-level effects




	###########COMPUTING GAMs for smooths##############
	# Following Crainiceanu, C. M., Ruppert, D. & Wand, M. P. (2005). Bayesian Analysis for Penalized Spline Regression Using WinBUGS. Journal of Statistical Softare, 14 (14), 1-24.
	# X.basis is data computed in R

	tauX~dgamma(1.0E-2,1.0E-4) #alternate prior, original from Cainiceanu et al. second gamma parameter == 0.0001 << (abs(mean(B.X[]))^2)/2, mean(B.X[]) ~ 0.2
	sdX <- 1/(pow(tauX,0.5)) # ~ dunif(0,5)
	taubeta ~ dgamma(2,0.2) # prior on precision of gam coefficients(
	sdbeta <- 1/(pow(taubeta,0.5))

	for(k in 1:nknots)
	{
		# Computation of GAM components
		B.X[k] ~ dnorm(0,tauX)

		for(i in 1:nstrata)
		{
			#beta.X.p[i,k] ~ dnorm(0,taubeta) #alternative non-centered parameterization
      #beta.X[i,k] <- B.X[k]+beta.X.p[i,k]
      beta.X[i,k] ~ dnorm(B.X[k],taubeta)


			for ( t in ymin : ymax )
			{
				X.part[i,k,t] <- beta.X[i,k]*(X.basis[t,k])
			}#t
		}#i
	}#k

		for(i in 1:nstrata)
		{
		for (t in ymin : ymax )
	{
			smooth[i,t] <- sum(X.part[i,1:nknots,t])
		}#t
	}#i


	#-------------------------------------------------#



	#### summary statistics  ######
	sdn_ret <- 0.5*sdnoise*sdnoise
	sdobs_ret <- 0.5*sdobs*sdobs

	for( i in 1 : nstrata )
	{

		for( t in ymin : ymax )
		{
			for(o in 1 : nobservers[i])
			{
				no[i,t,o] <- exp(strata[i]+smooth[i,t] + obs[i,o] + sdn_ret)
			}

			mn[i,t] <- mean(no[i,t,(1 : nobservers[i])])
			n[i,t] <- nonzeroweight[i]*(mn[i,t])
			n2[i,t] <- nonzeroweight[i]*exp(strata[i]+smooth[i,t] + sdn_ret + sdobs_ret) #n2 is an alternative way of calculating n

		}
	}


	#-------------------------------------------------#
}
