《时间序列干扰分析非重复大尺度实验.docx》由会员分享，可在线阅读，更多相关《时间序列干扰分析非重复大尺度实验.docx（19页珍藏版）》请在三一办公上搜索。
1、第9章: 时间序列干扰分析：非重复大尺度实验Paul W. RasmussenDannis M. HeiseyErik V. NordheimThomas M. Frost91 非重复的研究一些重要的生态学问题，尤其是那些在大尺度或独特尺度下的生态学问题，重复和随机几乎不可能 (Schindler 1987; Frost 等 1988; Carpenter 等 1995)。例如，假如我们想知道湖水酸化是如何影响龟甲轮虫属种群, 那么我们主要的工作就是将单独的小湖泊进行酸化，因此，这样的操作也许只能在一个湖中进行。即使是有酸化前的基础数据，这样的实验也不可能重复地进行，因此不能够采用传统统计方法
2、来进行分析，例如方差分析（第4章）。一种可替代的方法就是采用一些生物学或物理学模型，这些有趣的系统模型允许进行重复。例如，我们可以在湖中建一些小而相同的围栏并将它们中的一些进行酸化。尽管这（在实施适当的情况下）能够使我们对这些小实验单元的差异进行有效的统计分析，但是，这种模型能否在湖泊生态系统中进行有效的生态学普遍应用仍有置疑（Schindler 1987）。对于大尺度现象，小尺度模型的实验应用于大尺度研究仍然不能令人信服（尽管他们可能提供有价值的补充信息）。 我们将在本章检验一些非重复研究类型是如何利用以时间序列数据发展起来的技术进行分析。时间序列是取自相同的实验单元在不同时间进行的重复测量
3、或采集子样本。时间序列分析技术包括能确定一个序列平均水平上非随机变化是否已经在以前特定时间内发生的方法。这种分析的结果能够帮助确定那些特定时间之前发生的其他变化或操作是否可能导致了观察系列中的变化。我们描述的这些方法利用了许多可用于大尺度研究中的长期数据（还见于Jassby和Powell 1990; Stow等 1980）。以前的一些作者已经提供了许多的评价这些研究的技术，包括从图示的方法（Bormann和Likens 1979）到更加复杂的统计分析 (Matson和Carpenter 1990)。这些技术都强调了时间序列数据，并通常需要对处理和参照系统中处理前后量测数据序列进行比较（例如，S
4、tewart-Oaten et al.1986; Carpenter et al. 1989）。尽管有很多其他的时间序列方法也可以使用，但我们将着重讨论时间序列方法在ARIMA模型技术的应用。我们首先考虑时间系列数据如何影响传统统计分析方法的应用和解释（第9.2节）。第9.3节中举出一些非重复时间序列设计的例子。第9.4节介绍一些时间系列分析的重要概念和思想。第9.5节描述干扰分析，这是一个时间序列分析的扩展，对检验人工处理或自然干扰的时间序列的影响是有用的。第9.6节举例说明干扰分析在整个湖水酸化实验中的应用。（如果你不是一个水生生态学家，那么你可以在你的脑子中用你喜欢的生物体和生态系统分别
5、来代替“龟甲轮虫”或“湖泊”，在我们的实例中，原理和潜在的应用都是通用的。）我们将在第9.7节讨论一些有普遍性的议题。92 重复和实验中的误差我们首先简要地回顾一下湖水酸化的例子。对于检验湖水酸化对龟甲轮虫密度影响的传统实验设计如下所示： 1 确认我们想要将进行推论应用的所有湖泊。2 从这个湖泊总体中随机选择四个湖。 3 从这四个湖中随机选择两个进行酸化。4 将其它两个湖作为参照，或控制湖泊。假设已有对四个湖泊的处理前后龟甲轮虫密度的估计；对每一个湖泊，我们能够计算处理前后之间的变化（差异）从而可以进一步分析。然后再考虑到典型湖泊变异性后，继续用传统的分析方法来检验酸化湖泊中的变化是否始终与参
6、照湖泊中的变化不同。这种典型湖泊变异性，或试验误差，仅仅能够利用每一组内已有重复湖泊进行估计。（要想评估试验误差，每一组最少两个湖泊，尽管每组仅有两个湖泊的设计通常得到相当粗略的实验误差估计Carpenter1989）。传统方法的一个重要的特征是我们能够详细给出当酸化不存在的情况下却错误地认定有酸化影响的概率（类错误的概率 译者注）。尽管将更多的湖泊进行酸化是很困难的，但从单独酸化处理的湖泊得到处理前后的长期数据确是可行的。使用这样的序列，根据典型湖泊的可变性判断，我们能够检验与是否在处理时数据序列发生了一些不寻常变化相关问题，这可与全部时间内数据序列发生正常变化进行比较做出判断。在这里检验的
7、假设要比传统的假设要弱。这种假设检验仅仅能够回答在处理时是否发生变化这类的问题；而不能解决这种变化是由于处理还是由于一些其他偶然的事件引起这样的问题（Frost等 1988; Carpenter等 1998）。无重复情况下由于处理而导致变化的案例基本上是类似于统计学问题的生态学问题；这些争论通常能够被那些证实了的数据所支持，例如围栏试验。在下一节，我们要讨论时间系列的设计，这种设计在某种程度上防范了对假变化的认可。 确定某一时间内时间序列的正常变化要比确定一个重复实验中的实验误差更困难。在一个随机重复实验中，假定实验单元是互相独立的，这样就能够极大地将统计模型简单化。另一方面，时间序列的量测通
8、常互相依赖或自相关；即将来能（至少部分地）从过去来预测。时间上的系统变化依赖于其自相关结构，并且时间序列分析大都集中在对这些结构模型的构建和评估上。我们提醒读者要防止将时间上的亚抽样（subsampling）视为真正的重复并在分析非重复的试验中当作是有重复的。在一些例子中，可能已经证明一些传统的检验如t-test能够应用到时间序列上。如果分析表明量测数据间不是自相关，那么这样的检验就能够被采用（Stewart-Oaten et al. 1986）。尽管如此，一定要记住这种被检验的假设必须在实施处理时对检验水平有所改变，而这种变化可能并不是由于处理而产生的。要想有效地实施生态学实验，就必须对这些
9、问题有一个的基本的理解；这方面的讨论可参看Hurlburt (1984) 和Stewart-Oaten 等 (1992)。照例我们不应该鼓吹非重复的设计，但是它们有些时候是必要的。尽管时间序列有很具体的方法，在解释结果时仍需谨慎。没有重复，我们就不能够保证我们对用于判断处理效果的基本误差有充足的了解。然而，有几个方法可以增强在这方面解释的置信度。这包括1）考虑自然的可变性的信息，这些信息是被操作系统类型所特有的 2）在大尺度的操作范围内融入小尺度的实验 3）发展评估过程的仿真、机理模型（Frost 等 1988）。93 非重复时间序列设计假定一个生态学家在一个湖泊中监控龟甲轮虫密度5年，每两周
10、观测一次，然后将湖水酸化并继续监控该湖泊5年。在这个例子中，对单个的实验单元（湖泊）随着时间的观察是有顺序的（时间序列）。我们将这样一个干扰内的单独序列称为前后设计。（酸化就是干扰）有时可能并不会像预期的那样，自然界会对某个系统施加干扰。这可称作是自然实验。在这种情况下，可能没有其它选择，不得不采用简单的前后设计。这类设计最倾向于出现会被错误地认为是干预影响的巧合“跳跃”（数据的突然变化译者注）或趋势。当处理是在观察者的控制下时，就有可能改善前后设计。一个简单的改进是设立多重干扰，或在处理和控制条件之间来回变动。如果每一个干扰都得到一个一致的响应，那么观察到的反应仅仅是巧合的可能性就会减少。采
11、用配对单元的设计也是有用的。一个单元在研究的过程中接受了一个处理（干扰），另一单元则作为参照或基础值而不施以处理。处理前后在这两个单元内同时重复取样。Stewart-Oaten等(1986) 已经讨论过这种设计，而更早的讨论见于Eberhardt (1976)，虽然这种方法已经在这两个讨论之前被野外生物学家们应用（见Hunt (1976)）。我们依Stewart-Oaten等(1986)将这种设计称为前后控制影响（BACI）设计。在一个BACI设计中，参照单元提供一个与处理单元进行比较的标准。这就有助于确定在处理单元中的变化是否是由于处理本身，长期的环境趋势，自然的变化，或其他一些原因所致。在
12、使用BACI设计时应考虑两个重要因素，即实验单元，时机的选择以及确定每一单元的样本数量。处理单元和对照单元应该在所研究系统中具有代表性，从而在这些系统中，观察到处理的影响能代表整个系统，并且它们相互之间物理和生物特性应该是相似的。常见的做法是样品采集间隔是均匀的，采样时间间隔依赖于所研究种群或现象变化的速率（Frost等1988）。例如，种群变化迅速的昆虫或浮游动物需要一年采样多次，而变化慢的脊椎动物或植物种群通常一年仅采样一次。采样的频率依赖于所研究的系统和所提出的具体问题。每一个评估时期的持续时间（处理前后）是由自然周期、随机时间波动幅度，和所研究的生物体的生命周期决定的。当随机波动很大、
13、循环周期更长或研究的生物体是长寿时，就需要更长一些的采样间隔时间。我们将采用一个多重干扰的BACI设计的例子来详细展示一个真正的非重复设计的一些分析和解释的细节。这个例子关注的是小石湖酸化的影响，小石湖是位于威斯康辛北部的一个小贫养渗出湖（Brezonik等1996）。从1983年开始，这个湖已经成为研究渗出湖对逐渐酸化反应的一个生态系统水平试验的地点（Watres和Frost 1989）。小石湖被不透水乙烯基塑料布分成处理和对照两部分。在经过一个观察的起始期后（1983年8月到1985年4月），处理部分逐渐用硫酸将其酸化成三个pH值目标水平，5.6 、5.1、 4.7，每一个pH值水平保持两
14、年。而对照部分在整个实验过程中保持平均pH值 6.1水平。详细的湖泊学特征可以参看Brezonik等 (1986)。我们在这里的分析将集中在随着酸化，龟甲轮虫（Keratella taurocephala）种群将会发生明显的变化。制定一些好的非重复设计需要一些创造力。要分析好结果数据也同样如此。对大多数有真正重复的研究来说分析方法（例如方差分析）相当简单。然而，大多数无重复试验的分析就不那么清晰，且常趋向于主观。如果时间序列的数据长度足够的话（至少50个观测资料），那么就有可能采用具有严格理论基础的统计技术。我们将在9.4章节 和9.5章节讨论一些这样的技术。这些统计模型通常的目标是发现一个简
15、化的模型、一个能够很好的处理数据的简单模型。简单的说，我们想要一个生物学和物理学都合理的模型，且该模型能够利用较少的参数来进行评估。我们将对单个干扰的时间序列采用恰当的分析技术。在接下来的两部分章节也将看到，正确地采用这些技术需要掌握某一定的基本技能并且要非常小心。除非你在这个领域已经很熟练了，否则我们建议在进行分析时请教一下统计学家。94 时间序列时间序列分析包括大量统计技术用于分析具有时间次序序列依赖性的观察结果。我们将集中于讨论以自回归综合移动均值模型（AutoRegressive Integrated Moving Average Models-ARIMA）发展的时间序列技术的一个子集
16、。对于ARIMA模型的经典参考文献请见Box和Jenkins (1976)。其他的参考文献有MaCleary和Hay (1980), Cryer(1986) ,Diggle (1990), 和Wei (1990)。ARIMA模型很有价值，因为他们能够仅仅利用很少的参数来描述广泛的过程（由随机模型产生的事件次序）。此外，用于识别，估计和检查ARIMA模型适应性的方法已经有很好的研究。这些模型在模拟离散、(通常)均匀时空间隔观察结果的时间序列时是很合适的。ARIMA模型包括两个基本类型的模型：自回归（AR）和移动均值模型（MA）。AR模型与时间t的观察结果与早期观察结果一起进行回归。首先，我们介绍
17、一些符号含义。设yt是初始时间序列，也就是t时间的龟甲轮虫密度。用以均值为中心的序列来论述较为方便，也就是说，Zt=yt-u; Zt是t时间密度对长期平均值的离差。则用于中心龟甲轮虫密度的p=2阶AR模型（标示为（）有以下形式:Zt = 1Zt-1 + 2Zt-2 +t（9.1）这里的 是系数（如同回归系数），t是在t时间的随机误差（通常假定与0均值和2方差不相关）。这个模型所表述的是目前龟甲轮虫种群密度是以前两次取样的密度加上随机误差的线性方程。MA模型将时间t观察的随机误差与现在和过去的随机误差t联系起来。例如， q=2阶的MA模型可表示为：Zt =t1t-1 +2t-2 （9.2）这里是
18、系数。更多通用的模型可包括AR 和MA，所谓的ARMA模型（注意，ARIMA模型包括ARMA模型，以下同）。拟合ARMA模型的作用就是用自回归和移动均值来描述所有的序列依赖性从而使残差或估计误差项看上去像无相关的随机误差或白噪音。ARMA过程要求在观察的整个时间内，无论什么样的时间间隔，都要有相同的均值、方差和自相关格局。 这是对所说的稳态(stationarity)基本特性的一种直觉的描述（更为严格的定义请见引用文献）。尽管缺乏真正的重复，稳态允许有效的统计推论和估计。就像重复是建立在过程中一样，稳态过程展示了不同的时间间隔下相同的行为。因为观察结果之间的自相关依赖于它们之间时间段的数目，所
19、以观察结果通常一定是有相等的空间间隔。然而对于一些生物学过程，这是否有必要还不清楚。例如，在浮游动物丰度之间的两周自相关系数在种群变化慢的冬天就要大于种群变化快的夏天。因为在自然界中许多的观察过程是非稳态的，因此找到一种修正观察时间序列的方式就十分必要，以便于修正后的序列是稳态的。如此，有很少参数的ARMA模型就能够用于修正后的数据了。要做到这一点有两种方法。一是将一些确定性的方程引入模型，例如整个时间内的线性趋势、已知时间里加一步或周期性的方程，例如正余弦，来表示季节性行为。另外的一种方法是针对序列进行差分，也就是说，计算新的观察结果如Zt-Zt-1(第一个差异)或者Zt-Zt-12（与周期
20、12的季节性变化）。差分是一个较通用的方法，因为它能够同时表示确定性的和随机的趋势。在另一方面，当确定性的趋势有特殊意义时，通过差分来拟合这些趋势要比简单的处理掉它们可能更合理。由于它的较大通用性，我们将集中讨论通过差分来得到稳态。描述一些过程可能需要AR和MA二者的参数以及差分。这就产生了有pAR参数、qMA参数和d差分ARIMA模型（标示为ARIMA(p,d,q)）。例如，ARIMA(1,1,1)模型有如下形式：Xt =1Xt-1+t1t-1 （9.3）同前述，式中xt = zt zt-1和t是不相关的误差。Box和Jenkins (1976)为确定ARIMA模型来自数据的恰当形式发展了一
21、些方法。他们的方法要求计算样本的自相关系数和相关的来自数据的方程，并利用这些方程的相关性质来确定可能的ARIMA模型。每一个样本在时间距k的自相关系数rk（即相距时间段观测值间的相关系数）在时间距为 1,2 被计算出来以获得样本的自相关方程（ACF）。一个相关的方程，样本的偏自相关方程（PACF）表示当调整为在中等时间距状态下的自相关时时间距为时样本观察结果间的自相关关系。此模型确定过程包括确定用于拟合ACF、PACF结果的ARIMA模型的子集并将原始序列标绘出来。原始序列可能会显示出长期趋势或季节性行为，因而意味着非稳态。对于非稳态序列的ACF也会非常缓慢地下降到零。如果序列出现非稳态，那么
22、它应该在进一步检验ACF和PACF之前差分，直至它稳定为止。理论上ACF 和PACF对任何特定ARIMA模型都给出可鉴别的信号。在ARIMA模型中，样本ACF 和PACF 是从希望能够辨认出这种信号的数据（可以是差分出来的数据）中来估计的，并因此确定恰当的ARIMA参数化方式。对于纯MA模型的理论ACF仅在低时间距的情况下才有大值，而PACF下降更为缓慢。具自相关的非零时间距数目给出MA模型的阶q。对于纯AR模型来说则正相反，在那里只有ACF下降更为缓慢时，PACF在低时间距的情况下才有大值。在这种情况下，非零偏自相关系数的数目确定AR模型的阶数p。通常低阶（1或2）AR或MA模型对观察到的序
23、列有好的拟合，但可能要求模型同时有AR 和MA二项。这个模型的拟和程度将通过检验残差来实现。残差样本ACF 和PACF在任何时间距都不会有大值，也不显示任何明显格局。Box 和Jenkins（1970）提出了一个对缺乏拟合度的全面检验，后来该检验已经被精炼了(Ljung和Box 1978）。95 干扰模型干扰分析扩展了ARIMA建模方法以研究已知事件或干扰对时间序列的影响。干扰或处理的反应可有不同的形式。两个最简单的形式是：一个是干扰后形成在新水平上的永久跳跃，另一个是干扰时暂时的波动，或峰值。例如，在ARIMA模型中一个步骤变化能够表示为：Zt = St + Nt （9.4）式中St=0表示
24、干扰之前，St=1表示干扰时和干扰后，是系数。Nt 是一个ARIMA模型。在这个公式里，如果原始数列不是稳态的，Zt可以使用差分出来的序列。要模拟一个有峰值或波动的反应，令St在干扰时为1，否则为0。Box和Tiao (1975)讨论了一些更复杂的模型，如线性或非线性增加到新水平。McCleary和Hay (1980)给出一些说明的例子。干扰和反应之间的时间距也能够通过用不同的时间距拟合模型来检测。一般没有直接的方法可以通过观察时间序列本身来确定干扰模型ARIMA形式。通常一个可能的ARIMA模型是从一个整体序列或者是单独地从干扰前后的序列得到。对于干扰的反应形式可以通过检验没有干扰的ARIM
25、A模型的残差，或通过关于期望反应的理论知识引出。通常一些可能的模型必须通过检验并从中选择最好的。如果干扰形式未知，McCleary 和 Hay (1980)建议首先应试一下峰值模型，然后逐渐地增加模型，最后是步骤模型。Box和 Jenkins (1976)已经强调了以迭代方法来建立模型，这样做的目的是为了得到一个简单且充分拟合数据的模型，同时这个模型能够具有合理的科学解释。按照这个迭代方法建立模型的思路，可以产生一系列的决定，从而以不同的方法对不同数据集进行分析。尽管我们能够给出干扰分析的一般步骤的框架，但是我们不能确切地描述适应所有数据集的方法。尽管每一个分析方法都是其独特性，但下面的步骤在
26、进行干扰分析时是有用的：1 划分时间数列。2 如果必要的话，对原始序列进行转换。3 决定序列是否应该差分。4 检验转换（很有可能）和差分的序列样本ACF 和PACF以确定可能的模型。5 对模型进行不断的迭代拟合，并评估他们的拟合程度直到获得一个好的模型。我们强调的是依赖于数据集进行分析的方法。接下来的例子不应该被视为一个干扰分析的药方，而应该作为通用过程的一个例子。在生态学和环境监测中应用到的其它干扰分析的例子已经在Pallesen 等(1985)，Madenjian等(1986), van Latestejin和Lambeck (1986), Noakes (1986), Bautista等
27、 (1992), Rudstam等(1993)的几本著作中有所讨论。96例子961数据我们将分析威斯康辛州小石湖龟甲轮虫丰富度数据(每升动物数量)。采样方法在Gonzalez等. (1990)的文章中已有描述。在湖水没有结冰时每两周采集一次数据，当湖水结冰后每五周采集一次。在头一年的监测期间这个安排有略微的变化，在所有年份中会不可避免的得到相同采样时间表时，我们放弃额外的观察结果或将一个观察结果用于两个连续的取样时间（我们放弃了所有八个观察值，把六个观察值用两次）。接下来我们将讨论观察结果的时间距影响，同时在9.6.5 节讨论获得等时间距的方法。本例的数据集包括每年19个观察资料（总共106个
28、数据），所有的观察在每年差不多相同的日子进行。我们将检验两个干扰的影响：pH值在1985年4月29日从6.1降到5.6，在1987年4月27日从5.6降到5.1。我们没有采用在1989年5月9日第三次干扰后收集到的数据。任何数据分析的早期步骤之一是应该仔细的检查数据的图表。所有两部分丰富度序列的图表表明整个时期内对照部分的龟甲轮虫丰富度发生了很少的改变，但是在干扰之后处理部分的龟甲轮虫丰富度和变异水平都有所增加。962 导出序列首先要确定的是分析两个分开的序列还是分析从两个原始序列导出的一个单独序列。这两个方法各有其优点。单一导出序列在分析上可能简单一些：这不仅是因为只有一个序列，而且因为导出
29、序列可以比原始数列有较少的连续自相关和较不显著的季节性行为。另外，导出序列的分析可以得到一个对干扰影响的单独直接检验，而处理序列和对照序列的分开分析则需要综合在一起考虑来评估干扰影响。另一方面，导出序列更远离观察数据，并且我们也许希望模拟原始序列本身的连续依赖性、季节性行为和干扰影响。如果两个序列的时间格局差异很大，那就有必要分开来分析两个序列。在许多案例中，这两种方法都是正确的。这里我们都将加以讨论。导出序列有两种可能的形式：1）每一取样日处理和对照的观察数据间差异的序列2）每一取样日处理和对照的观察数据比率的序列。在生态学文献中关于这个题目的绝大多数讨论已经得出结论，即对数转换比率序列(等
30、同于两对数序列之差)最能够满足统计检验的假定(Stewart-Oaten等 1986, Carpenter等 1989, Eberhart和Thomas 1991)。这个观点部分考虑到了对导致反应改变的因素是叠加的还是相乘的评估。因为影响种群的因素经常具有相乘的影响，比率序列最适合于丰富度数据分析。通常对两种导出序列最好都要进行检验。对龟甲轮虫数据来说，随着每次干扰后序列均值和变异水平逐渐增加，两种导出序列显示出相似的格局。变异的增加在差异序列中更为显著，我们选择比率序列进行分析是因为我们可以找到一种可使该序列方差更接近常量的转换。注意我们对每个值都加1，以避免被零除（注意我们为了避免被零除，
31、对每个值都加1）。963 转换如同许多标准统计模型一样，ARIMA模型的一个基本假定就是观察数据方差恒定。这种假定最常见的违例是观察数据的方差与均值之间有相关关系。当原始序列方差与均值某次幂成比例时，我们将用一个简单方法来确定一种稳定方差的转换(Box等 1978)。Poole（1978）和Jenkins(1979)描述了该方法在时间序列分析中的使用。Jenkins(1979)建议根据季节的长度将序列分成若干从4到12个观测时间段的子集（我们采用6个月的间隔）。然后计算每一个子集的均值和方差，对数标准差对对数均值的回归。回归的斜率（b）值表明恰当的转换（见附录9.1的SAS语句）。方差与均值某
32、次幂成比例的假定决定了转换的形式：z = y1-b(1-b)0（9.5） = log yb = 1对比率序列而言，对数标准差对对数均值的回归斜率是1.33。1.33的估计有一些与之具来的误差；传统观念上我们不愿采用准确的估计而愿意采用与更易于解释的转换相对应的近似估计。因此，我们可能考虑平方根倒数转换（原数值平方根的倒数；斜率为1.5时适用）或者图 9.1威斯康辛州小石湖酸化处理过程中龟甲轮虫（Keratella taurocephala)的时间序列。（A）原尺度处理和对照部分龟甲轮虫丰度序列，头/升。（B）同样序列以10 为底的对数尺度。（C）处理序列与对照序列比率以10为底的对数尺度。对数
33、转换（斜率为1.0时适用）。尽管平方根倒数尺度对方差进行稳定性处理要稍好一点，但是我们在两种尺度上进行分析发现了相似的结果。由于大家对对数转换更熟悉一些而且较容易解释，我们将用对数比率序列分析（我们用以10为底的对数）。964 ARIMA模型在log10尺度下两部分(图9.1B)的丰富度序列图展现了与丰富度序列(图9.1)非常不同的结果。在原始尺度上均值和变异都有很大增加的处理湖泊序列，在对数尺度里则增加很小。对照部分序列显示了在对数尺度里的缓慢增加趋势。对数比率序列（图9.1 C）显示在第一次干扰后好像水平有所增加，而且在第二次干扰后也许变异也有增加。当进行统计分析时这些主观的印象会保留在脑
34、海里。鉴别ARIMA模型的方法要求序列稳态。如果干扰改变了一个序列的水平，那么这个序列就不再是稳态了，标准模型识别程序就可能很难应用。决定干扰序列ARIMA形式的一个方法是使模型拟合两次干扰之间的序列部分。在讨论此对数比率序列过程的结果之前，我们先展示现有最长稳态序列，即对数对照序列，的模型识别过程。ARIMA 模型鉴别的第一步是检验序列的样本ACF和 PACF（附录9.2A）。样本ACF和PACF通常表示为图表（见表9.1的SAS结果打出页），在表中每个自相关或部分自相关是从零延伸，具有两个标准误（SE）限（大致95%的置信限）的一条线。对数对照序列样本ACF在时间距为1,2,3时具大值，而
35、其它地方的值比二倍SE区间要小。样本PACF在时间距为1时具有单一大值。由于PACF在时间距为1时切断，而ACF逐渐下降到零，这种格局意味着一个AR(1)过程。而AR(2)过程的PACF在时间距1和时间距2具有大值。ACF和PACF两个样本在时间距18或19几乎达到2SE极限，这意味着可能有一个弱年周期。一旦一个模型形式被暂定，下一步就是使该模型去拟合数据序列（附录9.2 B)并检验残差。在此例中，来自拟合AR(1)模型的残差具有随机不相关误差（白杂音）的特点：在样本ACF 和PACF 中没有大值或强的格局，残差是随机的和不相关的零假设检验不被拒绝无论任何时候ESTIMATE语句被用于去拟合一
36、个模型时，SAS程序ARIMA都会运行这个检验（ SAS Institute Inc. 1988年）。AR(1)模型产生了均值估计值1.44(SE=0.08)和自回归系数() 估计值0.65（SE=0.07），以及残差估计值0.084（在使用SAS ESTIMATE语句时，这些估计值和标准误出现在ARIMA打出的结果页中）在两次干扰之间对数比率序列的三个节段执行相同的过程表明这三个部分都能够很好的被AR(1)模型拟合。这就引出了下面的参数估计值（SE的值在圆括号里）：节段方差1: pH 6.1-0.10 (0.05)0.28 (0.18)0.042: pH 5.6 0.37 (0.09)0.4
37、7 (0.15)0.093: pH 5.1 0.54 (0.19)0.62 (0.13)0.21表9.1 出自附录7.2程序的自相关方程（ACF）和偏相关方程（PACF）的输出表。尽管这三个节段具有相同的AR(1)形式，但参数估计值和方差可能有差异。因为这些序列太短（31到38个观察数据之间），这些估计值可能比较差，但是这些差异也意味着当解释干扰模型时应该谨慎从事。三个节段均值间的差异为干扰的影响提供了估计值（见下面的讨论）。要注意的是，如果对照和处理序列开始时有相同均值，那么对数比率序列的估计均值在pH 6.1时应该是0。前述的分析表明我们能够通过假定一个AR(1)形式来开始用干扰模型拟合所
38、有序列。我们也必须详细说明干扰的形式。湖泊酸化是逐步进行的，因此，它们可以由变量来表示，干扰之前取值为0，干扰后为1。这些变量就像在回归中的假变量（dummy variables）（Draper和Smith 1981）。我们构建了两个假变量，每个对应一次干扰。第一个是在1985年4月29日至1987年4月27日之间取值为1，其他的则为0。在此参数化过程中，均值的参数用来估计序列的初始水平，干扰参数估计用来干扰之后水平与初始水平的离差。尽管干扰参数的解释会有所不同，但能够被采用的任何一个假变量都会导致相同的预测值。当用SAS过程的ARIMA程序ESTIMATE来拟合模型时，假变量可视为输入变量。
39、我们使干扰模型拟合对数比率序列时，也分别拟合对数处理部分序列和对数对照部分序列（附录9.3）。在所有案例中，尽管对于对数处理序列来说也有一些证据表明其具有年周期，但简单的AR(1)模型能够很好的进行拟合。我们还使模型能够拟合更复杂的干扰反应，但是增加的参数并没有显著区别于零。三个序列的结果如下表示如下（SE的值在每一估计值后圆括号里；1和2是对第一和第二干扰参数的系数）：序列方差比率-0.06(0.13) 0.57(0.08) 0.40(0.17) 0.60(0.18) 0.12处理 1.40(0.12) 0.56(0.08) 0.23(0.15) 0.69(0.16)0.09对照 1.46(
40、0.12) 0.60(0.08) -0.16(0.16) 0.09(0.16)0.08对于参数估计，我们能够用每一个估计的SE去除该估计值，就像线性回归一样构建与t统计量相似的统计。近似的检验能够以标准正态分布为基础，p-值能够通过利用统计表或软件包的Z分数来决定。所有三个序列的模型具有相似的自回归系数（），近似0.6值意味着所有三个序列均具强连续相关性。干扰系数提供了对自初始基础时期水平的变化估计。因而对于对数比率序列来说，0.4（）是从初始-0.06水平到第一次干扰之后的水平变化。并且0.6 (）是从初始水平到第二次干扰后的水平变化。这两个系数都是显著不同于0（Z=2.35,p=0.02;
41、 and Z=3.33, p=0.001）。在第二次干扰时的变化能够由2-1计算得到这个差异。对数比率序列均值水平最大的变化是第一次干扰，而在对数处理序列则是第二次干扰。要注意的是，处理和对照序列干扰参数之差近似等于比率序列的干扰参数即，0.23-（-0.16）=0.39和0.69-0.09=0.60。对于对数比率序列，第一次干扰时的较大变化是因为对数处理序列的增加和对数对照序列的减少，而第二个干扰时的对数比率序列的较小变化是由于对数对照和处理序列都增加的缘故一般情况下不应预料这种单个序列和其导出或转换的序列间的相加性，出现在对数序列中是由于log(x/y) = log(x) log(y)造成
42、的。当转换到原丰富度尺度时，干扰参数更容易解释。对数比率序列的初始预测均值水平是log10(t/r)=-0.06，意味着t=0.87r(t是处理部分值，r是对照部分值)。因而在第一次干扰前，处理丰富度是对照丰富度的0.87倍。在第一次干扰后，均值水平-0.06+0.4=0.34意味着t=2.19r,并且在第二次干扰后，均值水平变为-0.06+0.6=0.54，意味着t=3.47r。第一次干扰时对数序列0.23均值水平变化与100.231.70倍龟甲轮虫丰富度一致。第二次干扰的变化为初始丰富度的100.694.9倍。965 观察资料的间距我们以观察数据有相同间距为前提来分析我们的序列，而实际上夏
43、冬间的时间距是不同的。要研究观察数据的时间距是如何影响我们的分析，我们从对数比率序列构建了两个有相等时间距的观察数据序列，并用干扰模型拟合它们。用来构建有相等时间距观察数据序列的一些方法如下：1 内插缺失的观察数据。2 将数据整合而形成一个较长的取样间隔（Wei 1990）。3 采用统计模型来添补缺失的观察数据 (Jones 1980) 。我们用方法1通过在有5周采样间距的冬季观察数据间插数值来建立一个每年有24个观察数据的序列。我们也采用方法2通过平均每四周所有的观察数据来建立一个每年13个观察数据的序列。这两个建立的序列都能够充分地拟合AR(1)模型，且干扰影响的估计值与那些来自原始序列模
44、型的估计值非常相似。因为构建这些系列的过程包括滤波，所以这两个构建的序列的方差很小。自回归参数的大小随观察频率的增加而增加。这里重要的一点是干扰参数估计值受观察数据的时间距影响并不很大。因此，我们对于不等时间距没有歪曲干扰参数的估计感到满意。我们选择描述那些最符合实际观察数据集分析。966 与其他分析程序的比较如果我们忽视在整个时间内观察数据间的连续依赖性，我们就能够用t检验来比较干扰前后的数据。重新回忆一下，t检验以及相似的传统统计过程的基本假定是观察数据间相互独立。我们分析的龟甲轮虫序列在一个时间单位的观察数据有约0.57的自相关系数，且自相关系数随着观察数据时间距的延长成指数降低。当观察
45、数据正相关时，t-统计量就会增大，并且增大的程度随自相关系数的增加而增加(Box等 1978; Stewart-Oaten等 1986)。对于对数比率序列，用于比较干扰前后的干扰系数而计算得到的Z统计量只是通过标准t检验计算所得结果的1/3到1/2。在这个案例中，所有的检验都是显著的，但是在干扰具有较小影响的条件下，t检验就能导致错误的结论。随机化干扰分析（RIA）已被认为是分析实验设计数据的一个供选方法，正如我们在本章已经讨论过的(Carpenter等 1989)。随机化检验包括将观察统计量与通过计算那个观测值得到的分布进行比较，这个统计量是由原始观测数据的每个大量随机置换排列得来的（第7章
46、）。尽管随机化检验没有要求正态分布的假设，他们仍然能够受到观察数据间由于缺乏独立性而产生的负影响（第14章）。尽管Stewart-Oaten et al. (1992) 质疑这些评估的通用性，但是经验评估表明即使当观察数据有中度自相关RIA仍然能够很好的运转(Carpenter等 1989)。采用RIA方法对小石湖龟甲轮虫数据进行分析得到的结论与我们发现的结论是平行的(Gonzalez等 1990)。尽管当只有一个对照单位时这个方法并不可行，但Schindler等(1985)仍然采用了传统统计方法将单个控制系统与许多参照系统进行比较。当标准的监测能够从一些对照系统中得到数据时这也许是有用的（U
47、nderwood 1994）。967生态学结论根据我们的干扰分析，我们得出这样的结论：在小石湖中酸化处理的部分龟甲轮虫种群呈非随机化变化。然而，问题是这些变化是否可以归因于酸化或是否这些变化也反映了与pH操作一致的自然种群变化。小石湖酸化的增加显示了对pH处理的一种反应，龟甲轮虫种群的生态学信息对这一结论提供了强有力的支持。纵览各种湖泊样品，龟甲轮虫显示出一种对低pH值环境的强烈亲和力（MacIssac等 1987）。先前关于龟甲轮虫分布的信息导致下列预测，即在任何的酸化之前，龟甲轮虫的种群将在较低的pH水平下增加（Brezonik等 1986）。对小石湖的处理部分的龟甲轮虫对环境反应的详细机
48、制分析揭示了龟甲轮虫在更酸化的环境下种群增加是与它的捕食者种群减少相关（Gonzalez和Frost 1994）。97讨论干扰分析提供了其它分析相似时间序列数据的分析方法所不具有的三个优点。首先，即使不同观察数据间存在自相关，它也能提供对干扰影响的好的估计值和标准误。第二，干扰分析要求对数据的自相关结构有明晰的模仿，这也可提供额外的关于潜在过程的信息。第三，干扰分析所固有的迭代模型建立方法本身是一种了解数据特点的途径，否则这些数据特点可能很难被发现。与那些简单执行标准检验来计算一个p-值不同，干扰分析的确是一种很不同的方法。干扰分析的不利因素是它要求具有长序列观察结果（通常50或更多）来确定连续依赖格局，且在数据分析上花费大量时间，包括熟悉方法和执行每一步分析。我们没有涉及到干
