Supplement 1. WinBUGS code for the statistical model that estimates vital rates, effects of time since fire, year and population effects.

File List
Supplement1Code.txt

Description
The code is commented throughout (following pound signs, ###).  Here we make note of indexing and notation issues, peculiarities of the data, and BUGS conventions that may be confusing.  Notation is consistent between the manuscript and BUGS code below, with two exceptions: (1) transition probabilities (ap,t,i,j) described in the manuscript are here denoted pi[p,t,i,j] and (2) the notation for random effects follows the parameterization described in Appendix A, Section 1.4, “Hierarchical Centering”.

The seedling survival submodel, a logistic regression of the number of seedlings surviving from the time of detection (at a quarterly census) to the time of the annual census in September, is split into several loops because there were certain years and populations when either 1) no seedlings emerged or 2) fire killed seedlings. In the first case, it is impossible to run a logistic regression, because the denominator (i.e., the number of seedlings that emerged) is zero. In the second case, we wish to exclude the direct effect of fire (100% mortality) on seedlings, in order to detect the subsequent effect of time-since-fire. Thus what would otherwise be a single loop:

for (p in 1:P){

    for (t in 2:20){

               Nsdlg[p,t] ~ dbin(s[p,t], Ngerm[p,t])

               logit(s[p,t]) <- b1.s*tsf[p,t] + yr.s[p,t]

  }

}

takes up many lines of code in several loops. The five study populations are indexed 1:5 (corresponding to population codes 2, 4, 10, 12, and 19 listed in Table 2).

In order to constrain the values taken on by the b2 parameters to +1 or -1, we use three lines of code. For example, b2.sd, the sign of the correlation between the effect of year variation on the number of seeds per branch and the model-wide year effect, is constrained as follows:

  b2.sd <- 2.0*I.b2.sd - 1.0

  I.b2.sd ~ dbern(pi.b2.sd)

  pi.b2.sd ~ dunif(0,1)

That is, a sample from a Bernoulli (coin flip) process is assigned to the intermediate parameter I.b2.sd. Thus I.b2.sd takes values 0 or 1, which is then transformed to values +1 or -1 by multiplying I.b2.sd by 2 and subtracting 1. The Bernoulli process has probability pi.b2.sd, which is assigned a uniform prior distribution (0,1).

            In BUGS, a normal distribution (such as the normal priors assigned for b0 and b1 terms, year effects, etc.) is defined in terms of mean and precision. Precision (denoted τ by convention) is the inverse of variance. To translate between precision and standard deviation (σ), we use the power function (pow in BUGS):

tau.b0.sd <- pow(sd.b0.sd, -2)

That is, τ = σ-2 or τ = 1/σ2.

          Finally, note that in the multinomial regression submodel, we used an indicator variable (ind[p,t]) to force BUGS to ignore (1) missing data or (2) the year that a population burned. Censusing of population 19 did not begin until 1990 (see Table 1), thus the indicator excludes the first two years of (missing) data from this population. The year after fire (tsf=1) is another case of missing data: the only plants in the population are newly emerged seedlings. Only in the following year (tsf=2) do those seedlings begin to transition into classes described by the multinomial model. We excluded years that a population burned (tsf=0) because we wish to exclude the direct effect of fire (100% mortality), as above in the seedling survival regression.

文件列表
Supplement1Code.txt

描述
本代码全程附带注释（以井号###开头）。下文将对索引与符号规则、数据集的特殊属性，以及可能造成混淆的BUGS软件使用约定进行说明。本文手稿与下述BUGS代码的符号体系基本一致，但存在两处例外：其一，手稿中提及的转移概率（transition probabilities）在代码中记为$pi[p,t,i,j]$；其二，随机效应（random effects）的符号遵循附录A第1.4节“分层中心化（Hierarchical Centering）”中的参数化方式。

幼苗存活子模型为逻辑回归（logistic regression）模型，用于拟合**季度普查时被检测到的幼苗**至**当年9月年度普查时**的存活数量。由于部分年份和种群存在两种特殊情况，该模型被拆分为多个循环：1）未出苗；2）火灾导致幼苗全部死亡。第一种情况中，由于分母（即出苗幼苗总数）为零，无法执行逻辑回归；第二种情况中，为了探究火灾后间隔时间的后续效应，我们需要排除火灾直接导致的全致死效应。因此，原本的单循环代码：

for (p in 1:P){
for (t in 2:20){
Nsdlg[p,t] ~ dbin(s[p,t], Ngerm[p,t])
logit(s[p,t]) <- b1.s*tsf[p,t] + yr.s[p,t]
}
}

在代码中被拆分为多个循环，占据了大量代码行。本研究的5个调查种群以索引1:5表示，分别对应表2中列出的种群代码2、4、10、12和19。

为了将$b_2$参数的取值约束为+1或-1，我们使用了三行代码。例如，针对“每年变异对每分支种子数量的效应与全模型年度效应之间的相关性符号”参数$b_{2.sd}$，其约束代码如下：

b2.sd <- 2.0*I.b2.sd - 1.0
I.b2.sd ~ dbern(pi.b2.sd)
pi.b2.sd ~ dunif(0,1)

具体而言，我们将伯努利过程（Bernoulli process）的抽样结果赋值给中间参数$I.b2.sd$，因此$I.b2.sd$的取值为0或1；随后通过将$I.b2.sd$乘以2再减1，将其转换为+1或-1。该伯努利过程的概率为$pi.b2.sd$，我们为其指定了均匀先验分布$U(0,1)$。

在BUGS软件中，正态分布（normal distribution）（例如为$b_0$、$b_1$项、年度效应等指定的正态先验）以均值和精度（precision）的形式定义。按照惯例，精度（记为$ au$）为方差的倒数。若要在精度与标准差（standard deviation, $sigma$）之间进行转换，可使用BUGS中的幂函数（pow）：

tau.b0.sd <- pow(sd.b0.sd, -2)

即$ au = sigma^{-2}$ 或 $ au = 1/sigma^2$。

最后需要说明的是，在多项回归（multinomial regression）子模型中，我们使用了指示变量（indicator variable）$ind[p,t]$，以实现两种功能：1）忽略缺失数据（missing data）；2）忽略种群发生火灾的年度。种群19的普查始于1990年（见表1），因此该指示变量会排除该种群最初两年的缺失数据。火灾发生后的次年（$tsf=1$）属于另一种缺失数据情况：此时种群中仅存在新出苗的幼苗，直到次年（$tsf=2$），这些幼苗才会开始转换为多项回归模型所定义的状态类别。与幼苗存活回归中的处理逻辑一致，我们排除了种群发生火灾的年份（$tsf=0$），以排除火灾直接导致的全致死效应。