在一个农作物新品种的培育中，要对该品种和与之竞争的其它品 种采集性能数据。这些性能数据包括与指定农作物有关的多种农艺学 特性的测量值；例如，对于Zea mays，要测量其粮食产量、粮食湿度 和作物倒伏。
评价一个新的农作物品种(下文中称为“品种”)潜在的商业 价值时，要将其农艺学性能与其它品种的农艺学性能进行对比。其它 对照品种包括，培育该品种的公司中已商业化和未商业化的品种，以 及竞争对手公司的商业化品种。注意，对现有的商业化品种也会进行 这种同样类型的评价，以确定它们是应当保留在市场上还是应当被培 育中的、更新的品种替代。
新品种和对照品种的农艺学性能数据来自多个试验点。这些试验 点通常广泛分布于试验所包括品种的适应区域。覆盖这些试验点的适 应区域往往很大，其范围达数百平方英里。例如，一个新的Zea mays 培育品种的试验区可能从衣阿华州西部到密执安州东部，从威斯康星 州中部到伊利诺伊州南部。
由于试验项目中的差异，对于一个指定的品种和它的竞争者，数 据往往相当“不平衡”，指的是并非指定组中的所有品种都出现在 所有的试验点上。考虑只有一对品种的试验数据，也就是新品种和单 一的竞争者，某些试验点两个品种都有，而其它点只有这两个品种之 一。
对这些性能数据进行分析，以便确定与对照品种相比，新品种具 有足够大性能优势的地理范围，从而证明它应当引入这些地区的市场。 在理想情况下，在整个试验区域中，与所有对照品种相比，正在培育 的品种都具有显著的性能优势。然而在某些情况下，一个品种可能仅 在一个地区具有性能优势，但是它仍然可以在该地区中满足相当大的 市场需求。所以，不仅要在整个试验区域将指定品种与其它对照品种 相比确定性能特征，而且也要在各个地区内这样做，这一点很重要。
不同地理位置或地区两个品种相对性能的差异来自所谓的“环 境影响的基因型”(Sprague&Eberhart，1976)。环境影响的基因型 是由不同品种对环境条件的反应差异造成的。这些环境条件可能包括 如昼长、温度、土壤湿度和病虫害压力。注意，术语“环境”可以 指指定一年的不同位置，或者指定位置的不同年度，或者位置与年度 的某种组合。
Bradley等人(1988)介绍了品种性能评价的传统统计分析方法。 这些传统方法通常是基于“位置匹配的”数据，也就是，对于一个 指定品种相对于一个对照品种，在分析中使用的数据仅仅来自两个品 种共同种植的试验点。采用一次配对的t试验，对两个品种没有性能 差异的假设进行检验。此外，为了获得指定地区相对性能的有关结论， 分析中只使用该地区内试验点的数据。
对于位置匹配的数据采用t试验的传统分析效率很低，至少有五 个理由。第一，它没有使用全部数据；它使用的数据仅仅来自两个品 种共同种植的试验点。第二，为了进行地区性能对比，它不使用所关 注地区以外的邻近区域的数据。第三，它不利用与所关注之性能特征 有关的共变量(如灌溉或土系)来帮助解释和预测这些差异。第四， 在分析模型中它只使用年内的数据；如果一个模型使用多年度的数据， 可以获得更稳健的结论。第五，它是基于一个不正确的假设：一个试 验点的观测结果独立于其它点的结果。
有两个概括的理由，说明以上列举的低效为什么限制了广阔区域 中的农作物评价。第一是数据的使用效率。要作出结论时，实验人员 自然希望使用全部数据。在当前的情况下，这包括一个品种在没有种 植其它品种之处的数据，包括与所关注地区类似区域的数据，包括有 关共变量的数据，还包括跨年度的数据。统计方法应当争取充分利用 所有可用的信息。传统分析受到限制的第二个理由是它基于数据独立 的经典假设。品种性能数据几乎毫无例外地不满足这种独立性假设， 使统计的结论无效，当数据并不真正表明差异时，往往使人推断出存 在品种差异。
现有文献综述：
现有文献的简要综述也强调现有技术中潜在的不足，本发明争取 解决这些问题。
对农场地块的实验应用空间统计的领域中，关系最密切的论文之 一、可能是最值得关注的工作是由Zimmerman和Harville(1991)完 成的。作者认为观测结果是一个随机田的表现，引入了所谓的随机田 线性模型(RFLM)。在这个模型中，趋势是由一个平均结构所模拟 的，而小尺度的相关性通过空间自相关结构模拟。参数估计是通过一 种似然方法进行。通过实际数据分析，作者试图表明他们的模型优于 最近邻域分析(NNA)方法。注意，他们的研究专门针对小区域估计 的情形，空间相关性的范围限于一个试验点。
在空间预测中使用共变量的情况下，另一篇值得关注的论文是由 Gotway和Hartford(1996)完成的。对粮食产量作为一个共变量的 数据，通过应用协克里金法预测土壤硝酸盐含量，作者提出在空间预 测中使用辅助或次要变量。将他们的方法应用于一个试验点的数据后， 他们表明了其方法的效益超过更传统的外表漂移方法。他们研究的范 围也是仅限于点上的预测。
涉及多地点产量试验的近期论文之一是由Cullis等人(1998)完 成的。在这篇论文中，作者提出了一个多环境初代品种试验中空间分 析的方法。这个方法使用最优线性无偏预测(BLUP)求取基因型效应 和环境影响的基因型效应，使用REML求取空间参数和方差分量。不 过，该方法对每个试验单独模拟协方差结构，也就是不考虑试验之间 的相关性。
在广阔的陆地区域将地质统计学应用于土壤化学研究的领域中， Yost，Uehara和Fox(1982a，1982b)是农业科学研究的先驱。在两 篇相继发表的论文中，他们论述了在夏威夷岛上进行土壤化学空间预 测的研究结果。不过，他们的研究限于只使用克里金方法，而不考虑 任何协克里金的方式。此外，在他们的研究中，除了使用选定的土壤 化学变量以外，没有使用其它类型的变量。在土壤科学领域值得一提 的另一项类似却相对更近期的工作是由Ovalles和Collins(1988)完 成的。在整个佛罗里达西北部大约覆盖380公里×100公里的区域中， 作者使用通用克里金方法研究选定土壤性质的空间变化，记载的自相 关范围大约为40km。在他们的研究中没有试图进行任何空间估计或 协克里金。
涉及农作物产量空间预测的论文之一是由Bhatti，Mulla和 Frazier(1991)完成的。作者使用一块面积大约为655米×366米的 实验田，研究小麦产量与土壤有机质和土壤磷含量的关系。预测产量 时使用了克里金和协克里金。然而，研究限于单一的实验田，所以不 包含任何模拟大尺度趋势的方式，而该趋势在广阔区域的试验中往往 存在。农作物产量空间分析领域的另一篇论文是由Brownie，Bowman 和Burton(1993)完成的。在这篇论文中，作者对比了三种可选的空 间方法：趋势分析、最近邻域Papadakis分析和相关误差分析，以研 究玉米(Zea mays)和大豆(Glycine max)产量数据的空间变化。与 农作物产量的其它现有研究的情况相同，他们的研究也是仅在点上的， 也就是不考虑各实验点之间的空间相关性。事实上，作者在其论文的 结论中也注意到，虽然分析中结合了来自多点的数据，文献中还没有 农作物产量数据的跨点空间分析。
有关农作物产量空间分析的其它论文在文献中也是存在的。 Bhatti，Mulla，Kooehler和Gurmani(1991)使用了半变量图来识别 农作物产量的空间自相关性。他们通过研究应用NNA前后的半变量 图，显示了NNA消除空间可变性的效果。他们研究的范围也是仅限于 点上分析，空间自相关的最大范围大约是20米。此外，他们的研究中 不包括空间估计或预测的方式。Wu等人(1997)就谷类培育试验中 消除产量数据的空间变化，对比了变量误差一阶差分(FD-EV)法 (Besag和Kempton，1986)与更传统的Papadakis最近邻域法和经典 的随机化完整区块(RCB)分析。不过，他们的方法不需要预先指定 的趋势和空间自相关结构模型。另外，他们的研究仅限于点上空间变 化。
Stroup，Baenziger和Mulitze(1994)使用了若干育种区的数据， 对于有效地消除空间可变性产生的干扰，就处理效果比较，对比了传 统的RCB分析、NNA和随机场线性模型分析(Zimmerman and Harville，1991)。自然，他们的空间相关性研究限于每个育种区之内 的数据。
Gotway和Stroup(1997)的一篇论文的独特之处在于，作者把 广义线性模型理论拓展到包括离散和分类的空间变量的空间估计和预 测。他们将其拓展应用于两个数据集，一个有关作物虫害，另一个有 关杂草计量。不过，与所有其它研究的情况相同，他们研究的范围及 其应用限于每个实验点之内的数据。
在以上几段中，回顾了与本发明有关领域中现有的文献。这些领 域中研究工作的总体主要可以划分为两个范畴：(1)地质统计学和空 间统计学应用于土壤科学领域的大区域和小区域试验，(2)地质统计 学和空间统计学应用于农田内小区域试验的农作物反应分析。现有文 献缺乏在大区域试验情况下农作物反应分折(如Zea mays的粮食产量) 领域的研究，其中空间相关性的范围要扩展到超出单个的实验点。
在农作物品种的性能评价中，根本问题是在多种环境下，也就是 覆盖多个试验点的广阔地理范围内得出结论。在文献综述中也应当注 意，现有文献的任何研究中都没有涉及空间估计和空间预测这二者的 使用。相反，正如下面将要详细讨论的，本发明涉及解决品种评价问 题的一种新方法，它不仅是基于多环境、大区域试验，它还具有两个 部分来回答两个明显的问题：(1)一个估计部分，回答的问题是一个 品种的长期平均性能或两个品种的长期平均性能差异，(2)一个预测 部分，回答的问题是在若干指定的时间(年度)点上和/或若干指定的 空间(地理)点上，一个品种的性能或两个品种的性能差异。通过通 用克里金方法和通用区块克里金方法，本发明背后的方法考虑了大尺 度的趋势，而且通过协克里金方法容易使用共变量。现有文献中的工 作没有把以上的全部功能结合在一个统一的方法中，以便研究农作物 品种的性能。
以上列举的问题并非面面俱到，对于目前已知的作物分析技术中 有损于有效性的许多问题，只是涉及了一部分。其它值得关注的问题 也可能存在；然而，上述问题应当足以表明现有技术中采用的方法总 而言之还是不能令人满意。

本发明的若干实施例采用一种称为线性混合模型的统计模型和地 质统计方法，由广阔区域的试验数据评价某个农作物品种的性能。评 价农作物性能可以通过测量商业上的重要特性，比如(但不限于)产 量、粮食湿度或作物倒伏。混合模型能够使用共变量，比如年度、土 壤类型、灌溉等等，作为固定效应，而且能够使用随机效应，比如试 验位置，这样做有助于解释农作物的性能。固定效应和随机效应不能 解释的剩余变化采用地质统计方法模拟。地质统计模型考虑了数据中 的空间自相关性，能够获得有效的置信区间，以评价估计和预测中的 不确定性。
本发明的实施例具有两个明显的部分：空间估计和空间预测。估 计部分可用于：
(1)估计一个品种在一个指定空间位置的平均性能(点估计)， 然后使用这些点估计结果构造
(a)一个品种在广阔区域(如多国或多州的集合)的平均性 能曲面(跨越若干年度)，以及
(b)两个品种的性能在广阔区域(如多国或多州采集)的平 均差异曲面(跨越若干年度)，
(2)估计一个品种在一个规定地理区域比如一个市场区(性能的 区块估计)跨越若干年度的平均性能，以及
(3)估计两个品种之间在一个规定地理区域比如一个市场区(相 对性能的区块估计)跨越若干年度的平均性能差异。
预测部分可用于：
(1)预测一个品种在一个指定年度一个指定空间位置的平均性能 (点预测)，然后使用这些点预测结果构造
(a)一个品种在广阔区域(如多国或多州的集合)每一年度 的性能预测曲面，
(b)两个品种在广阔区域(如多国或多州的集合)每一年度 的性能差异预测曲面，
(2)预测一个品种在一个规定地理区域比如一个市场区(区块性 能预测)每一年度的平均性能，以及
(3)预测两个品种之间在一个规定地理区域比如一个市场区(区 块相对性能预测)每一年度的平均性能差异。
决定选用新培育的、与其它品种相比具有最优性能的品种时，需 要如上所述的品种性能估计和预测。对于一个候选品种是否进入商业 化状态需要作出决定，如果进入，将它置于适当的地区，使之与竞争 对手相比表现良好。
从营销的角度看，需要进行一个品种相对性能的定量评价有两个 明显的理由，首先，为了向市场引入一个新品种而将该品种推进到商 业化状态，其次，为了决定由一个相对性能表现更好的新品种取代一 个现有的商业化品种。
以上的决定所依据的最重要的准则之一，是一个品种的期望长期 性能优势。评价品种在市场中的商业价值时需要这一点，其中度量一 个品种在市场中的价值要跨越它在市场中的整个寿命。通过评价一个 品种在一个指定位置上或一个指定地理区域中的长期相对性能，本发 明的估计部分提供了这个问题的答案。
在一个指定地点上或一个指定地理区域中，本发明的预测部分可 以评价在一个指定年度，一个品种的性能以及品种之间的性能差异。 确定一个品种跨越空间和/或时间情况下性能或相对性能的一致性时， 需要这种特定空间和/或特定时间的性能评价。对于一个商业化品种来 说，性能的一致性——在作物培育术语中称为“品种稳定性”—— 是一个极为必要的属性。
在广阔区域的试验中，本发明优于传统方法的长处包括以下几点。 第一，它并不要求数据是位置匹配的(数据集内的每个位置都包含两 个品种的观测结果)，而是使用全部试验点的数据，只要该点至少具 有所研究的品种之一。第二，为了进行地区性能评价，不限制它只使 用本地区内的数据，而是也使用周围地区的试验点上的数据。第三， 它易于加入共变量的信息，这些信息与所关注的主要特性有关(例如 当产量是所关注的主要特性时，有关土系的数据)。第四，它使用的 模型加入年度作为一个模型因素，因此可以容纳多年度的数据。第五， 它并不假设不同试验点的观测结果具有独立性，而是利用不同试验点 之间的空间相关性以提供改善的统计结论。
本发明的一个实施例中的步骤可归纳如下：
1.构建农作物在广阔区域中性能的一个数据库，包括试验点的空 间坐标、每个试验点上品种性能的测量值和每个试验点上有关共变量 的数据。
2.选择两个品种进行对比，例如，一个正准备推进到商业化状态 的品种“头领”和一个已经在市场中的品种“竞争者”。
3.通过目测以及采用基于帽子矩阵的统计试验和Cook距离方 法，从两个品种的数据中排除异常值。
4.为线性混合模型选择一个空间协方差模型。
5.估计空间协方差参数以及线性混合模型中其它固定效应和随 机效应的参数。
6.使用估计的参数
(a)估计跨年度平均值或按年预测，
(b)计算在地理区域中连续曲面或区块平均值，
(c)计算每个品种性能特征或性能特征的差异，以及
(d)求取估计结果或预测结果的标准差。
7.使用估计结果和预测结果及其标准差评价品种性能，例如评价 头领品种的相对性能，以决定将其推进到商业化状态。
一方面，本发明是一种在广阔区域中评价农作物品种性能的方法， 它使用的线性混合模型包含地质统计部分，还包括固定效应、随机效 应和共变量的参数。“广阔区域”仅仅表明一个试验品种假定的区 域。例如，“广阔区域”可能是指一个试验品种的适应区域，它可 能包括许多试验点。“农作物品种”是指培育一个指定的植物种类 以及业内周知的任何其它用途。构建一个广阔区域的数据库，其中包 括一个或多个农作物品种的试验点空间坐标以及一个或多个农作物品 种的性能特征值。“数据库”是指数据的任何集合。例如，一个“数 据库”可能指的是可搜索数据的一个电子形式的集合。“性能特征 值”是指与特定品种相关联的任何受关注的农艺学特征。比如“性 能特征值”可能指的是业内周知特征的任何数字。例如，它可能指的 是比如Zea mays的粮食产量、粮食湿度和作物倒伏。根据广阔区域数 据库中的数据拟合线性混合模型，从而估计固定效应、随机效应和共 变量的参数。在一个或多个指定空间位置中的每一个，使用参数估计 结果，估计农作物品种的长期期望性能。“估计长期期望性能”意 为在多个试验期间(如年度)估计一个品种农艺学特征的一个数值， 以及业内周知的该词组的正常用途。例如，它可能指的是对于该品种 从产生直到它退出市场的整个时间周期，估计一个特征的一个数值， 不过它也可以理解为该时间周期可随估计过程的需要而改变，而且它 可能比从产生到退出的时间短得多。
在其它方面，数据库进一步可能也包括共变量数据。估计长期期 望性能可能包括估计该农作物品种与另一个农作物品种之间的长期期 望性能差异。估计这些参数的过程可能包括约束最大似然方法。估计 长期期望性能的过程可能包括广义最小平方方法。本方法可能也包括 在参数估计之前，使用一种倍率或Cook距离的方法排除数据。本方 法可能也包括计算与长期期望性能相关联的一个标准差。本方法可能 也包括形成长期期望性能的一个输出。该输出可能包括文本输出。该 输出可能包括图形输出。该图形输出可能包括一个等值线图，表示长 期期望性能的连续曲面。
另一方面，本发明是一种在广阔区域中评价农作物品种性能的方 法，它使用的线性混合模型包含地质统计部分，还包括固定效应、随 机效应和共变量的参数。构建一个广阔区域的数据库，其中包括一个 或多个农作物品种的试验点空间坐标、试验点所处的地理区域以及一 个或多个农作物品种的性能特征值。根据广阔区域数据库中的数据拟 合线性混合模型，从而估计固定效应、随机效应和共变量的参数。在 一个或多个指定地理区域中的每一个，使用参数估计结果，估计农作 物品种的长期期望性能。
另一方面，本发明是一种在广阔区域中评价农作物品种性能的方 法，它使用的线性混合模型包含地质统计部分，还包括固定效应、随 机效应和共变量的参数。构建一个广阔区域的数据库，其中包括一个 或多个农作物品种的试验点空间坐标以及一个或多个农作物品种的性 能特征值。根据广阔区域数据库中的数据拟合线性混合模型，从而估 计固定效应、随机效应和共变量的参数。在一个或多个指定空间位置 中的每一个，在一个或多个指定期间中的每一个，使用参数估计结果， 预测农作物品种的平均性能。“预测平均性能”仅仅表明预测一个 品种的一个农艺学特征的一个数值以及业内周知的该词组的正常用 途。
在其它方面，估计这些参数的过程可能包括约束最大似然方法。 数据库可能也包括共变量数据。共变量数据可能包括一个反应变量。 预测平均性能的过程可能包括通用协克里金方法。共变量数据可能包 括固定效应。预测平均性能的过程可能包括通用克里金方法。预测平 均性能的过程可能包括预测该农作物品种与另一个农作物品种之间的 平均性能差异。本方法可能也包括在参数估计之前，使用一种倍率或 Cook距离的方法排除数据。本方法可能也包括计算与预测平均性能相 关联的一个标准差。本方法可能也包括形成预测平均性能的一个输出。
另一方面，本发明是一种在广阔区域中评价农作物品种性能的方 法，它使用的线性混合模型包含地质统计部分，还包括固定效应、随 机效应和共变量的参数。构建一个广阔区域的数据库，其中包括一个 或多个农作物品种的试验点空间坐标、试验点所处的地理区域以及一 个或多个农作物品种的性能特征值。根据广阔区域数据库中的数据拟 合线性混合模型，从而估计固定效应、随机效应和共变量的参数。在 一个或多个指定地理区域中的每一个，在一个或多个指定期间中的每 一个，使用参数估计结果，预测农作物品种的平均性能。
另一方面，本发明是一种杂交品种培育方法，培育一个杂交品种， 获得该杂交品种和对照杂交品种的性能数据。采用广义最小平方法使 一个三次多项式曲面拟合到每个杂交品种的性能数据上，使用球面变 量图模拟残差方差。新的杂交品种与对照杂交品种的性能进行对比。
另一方面，本发明是一种系统，包括一台计算机和一种程序。该 程序在该计算机上执行并包括以下程序代码：采用广义最小平方法使 一个三次多项式曲面拟合到每个杂交品种的性能数据上；使用球面变 量图模拟残差方差；以及新的杂交品种与对照杂交品种的性能进行对 比。
本公开文件的优点应当理解为，预测、区块预测、估计和点估计 能够以任何数目的各种排列进行组合，以获得有价值的性能评价。所 有这些组合都不脱离本发明的范围。
附图简要说明
本发明采用举例说明，但不受附图的限制。本专利的文件包含至 少一张彩色图。专利和商标事务所在收到请求和必需的费用之后，将 会提供本专利的拷贝，包括彩色附图。
图1是一个等值线图，显示1996年两个商业化杂交品种—— DK512和DK527——产量差异的预测曲面。图中显示了三条产量差异 等值线(-5、0和+5蒲式耳/英亩)和四条显著级别等值线(产量差 异±5％和±20％)，也显示了两个杂交品种的试验点以及市场区的边 界和名称。
图2是一个等值区域图，显示1996年在若干地理区域中的每一个 区域内，两个商业化杂交品种——DK512和DK527——产量差异的区 块预测。等值区域图上的图例既表明了产量差异的级别(-3、0和+ 3蒲式耳/英亩)，又表明了差异是否达到5％的显著级别。
图3是一个等值区域图，显示1996年在若干地理区域中的每一个 区域内，两个商业化杂交品种——DK512和DK527——产量差异的简 单均值。等值区域图上的图例既表明了产量差异的级别(-3、0和+ 3蒲式耳/英亩)，又表明了差异是否达到5％的显著级别。
图4、图5、图6表示的信息分别与图1、图2、图3相同，只是 时间为1997年。
图7、图8、图9表示的信息分别与图1、图2、图3相同，只是 时间为1998年。
图10是一个等值线图，显示1996年、1997年和1998年两个商业 化杂交品种——DK512和DK527——产量差异的跨年度估计曲面。图 中显示了三条产量差异等值线(-5、0和+5蒲式耳/英亩)和四条显 著级别等值线(产量差异±5％和±20％)，也显示了两个杂交品种的 试验点以及市场区的边界和名称。
图11是一个等值区域图，显示1996年、1997年和1998年在若干 地理区域中的每一个区域内，两个商业化杂交品种——DK512和 DK527——产量差异的跨年度区块估计。等值区域图上的图例既表明 了产量差异的级别(-3、0和+3蒲式耳/英亩)，又表明了差异是否 达到5％的显著级别。
图12是一个等值区域图，显示1996年、1997年和1998年在若干 地理区域中的每一个区域内，两个商业化杂交品种——DK512和 DK527——产量差异的跨年度简单均值。等值区域图上的图例既表明 了产量差异的级别(-3、0和+3蒲式耳/英亩)，又表明了差异是否 达到5％的显著级别。
示范实施例的说明
现在将要详细介绍在此公开之发明的优选实施例，它的实例在附 图中展示。
本发明涉及以上在背景技术中论述的传统统计分析方法中的所有 五个限制。第一，它并不要求位置匹配的数据，而是使用全部试验点 的数据，只要该点至少具有试验中的品种之一。第二，为了进行地区 性能对比，不限制它只使用本地区内的数据，而是也使用周围地区的 试验点上的数据。第三，它易于使用共变量。第四，它使用的模型加 入年度作为一个模型因素，因此可以容纳多年度的数据。第五，它并 不假设不同试验点的观测结果具有独立性，而是利用不同试验点之间 的空间相关性以提供改善的统计结论。
期望品种性能数据具有空间相关性是相当自然的，因为品种具有 对环境条件的敏感反应。然而，为了对品种性能数据尤其是产量数据 中空间相关的程度作出定量评价，对来自Monsanto(DEKALB)玉米研 究机构的Zea mays产量数据进行了一项研究。分析了三年中(1994- 1996)包括许多Zea mays品种的历史产量数据。分析结果表明，大约 80％品种的数据具有很强的空间自相关性，其证明是拟合的协方差模 型中很大的正范围(Bhatti，Mulla，Kooehler and Gurmani，1991)。 这一结果受到一个1994年和1995年Zea mays产量数据混合模型分析 的支持，那是一项似然比试验，分别对每个品种试验空间自相关的程 度。平均说来，75％的品种在其数据中表现出显著的自相关性。
本发明使用一种统计模型，它既能进行空间估计，又能进行空间 预测，这使它在品种性能评价中具有超过传统方法的某些明显的优点。 空间估计能够评价品种的长期平均性能，也能够评价不同品种的长期 平均性能差异。这种信息的价值在于确定一个指定的品种相对于其它 竞争品种是否具有某种长期的性能优势。可以在特定的位置确定，也 可以在更大的区域中确定。这一点对于评价该品种在市场上的商业价 值是必需的，因为度量一个品种的价值是要跨越它在市场上的整个寿 命。空间预测能够在指定的时间点(年度)和/或指定的空间点(地理 位置)或聚集在更大区域的空间点上，评价品种的性能和性能差异。 在确定一个品种跨越空间和/或时间情况下性能的一致性时，这种特定 空间和/或特定时间的性能评价是有用的。对于一个商业化品种来说， 性能的一致性——在作物培育术语中称为品种稳定性——是一个极为 必要的属性。采用品种性能的混合模型分析，可以解决这种稳定性问 题(Littell等人，1996，250页)。
本发明在对比品种性能的方法中引入了地质统计模型。地质统计 模型使用从一种地貌的不同位置获取的样本点，并由这些样本点产生 (插值出)一个连续曲面。这些样本点为某种现象的测量结果，比如 空气中的污染物质含量或地质岩芯的矿物级别。据信地质统计学的起 源有数个来源。使用协方差的插值——也称为最优线性无偏预测 (BLUP)——来自Wold(1938)、Kolmogorov(1941)、Wiener (1949)、Gandin(1959)和Goldberger(1962)。使用变量图的插 值归功于Gandin(1959，1963)和Matheron(1962，1969)。Cressie (1990)对地质统计学的起源给出了更多细节。地质统计学最初应用 于矿业领域，随后应用的领域诸如气象学和环境科学。近年来它已经 应用于农业有关的领域，比如土壤性质和农作物产量的模拟。
混合效应模型：
“混合效应模型”形成了理解本发明的背景，所以应当进行更 为详细的讨论，以向读者提供与在此讨论的实施例有关的背景信息。 考虑以下普通模型：
y＝Xθ+Dr+ε，                [1]
式中y为反应向量，θ为未知固定效应参数的向量，r为未知随机 效应参数的向量，X和D为已知结构的矩阵，ε为未知随机误差分量 的向量。因为这种模型具有固定和随机两种效应参数，所以称为混合 效应模型，或者更简单地称为“混合模型”。假设r和ε为正态分 布，其参数为
$E\left(\begin{array}{c}r\\ ϵ\end{array}\right)=\left(\begin{array}{c}0\\ 0\end{array}\right),---\left[2\right]$
以及
$\mathrm{var}\left(\begin{array}{c}r\\ ϵ\end{array}\right)=\left(\begin{array}{cc}G& 0\\ 0& \Sigma \end{array}\right).---\left[3\right]$
在以下段落中，将要说明地质统计模型与这种混合模型的关系。
地质统计学中的混合模型：
地质统计学假设的一个空间线性模型由下式给定
Z(s)＝μ(s)+ε(s)，              [4]
式中向量s包含空间坐标(如经纬度)，表示随机变量Z(s)的空 间位置。[4]式由两部分组成：一个确定部分μ(s)和一个随机部分ε(s)。 对于两个不同的点s和u，随机变量ε(s)和ε(u)是空间相关的。因为 这种相关性存在于同一随机过程ε(·)内的各项之间，在术语中称为 “自相关”。为了在统计模拟中进行复制，假设这种空间自相关仅仅 依赖于分离s和u的距离和方向，而与它们的确切位置无关。这种性 质在术语中称为“平稳性”。
在地质统计学中最普通的模型是[4]式的一种特殊情况，确定曲面 μ(s)不随位置s而变化。这种模型由下式给定
Z(s)＝μ+ε(s).                [5]
[4]式和[5]式的模型都可以用于预测，术语“预测”对应于随机 量的推断(Cressie，1993，15页)。这里，推断一词是指随机量的预 测以及与该预测相关联的不确定性。在[4]式的模型中，确定曲面μ(s) 随位置变化，对于某个位置s0，Z(s0)的预测在术语中称为“通用克 里金”。在[5]式的模型中，模型的确定部分是一个常数，不随位置变 化，Z(s0)的预测在术语中称为“普通克里金”。
[4]式和[5]式的模型也可以用于估计，术语“估计”是指固定或 确定量的推断(Cressie，1993，13页)。这里，推断是指固定量的估 计以及与该估计相关联的不确定性。注意，在[4]式的模型中估计指的 是确定μ(s)，而在[5]式的模型中估计指的是确定μ。
以上关于预测和估计的介绍涉及在若干点上的推断；也就是在特 定点上某些数量的预测或估计。预测和估计也可以应用于空间聚集或 称为区域的级别；例如在整个州。在这种情况下，在某个区域A中， 我们或者预测随机变量Z(s)的平均值，或者估计确定趋势曲面μ(s)的 平均值。这两种运算在术语中分别称为“区块预测”Z(A)≡∫AZ(s)d s/|A|和“区块估计”μ(A)≡∫Aμ(s)ds/|A|，式中|A|≡∫A1ds为A的 面积。
[4]式的模型可以写为一个普通线性模型，
z＝Xθ+ε                       [6]
表示一个品种所有随机变量的联合期望。这里，随机误差的方差 -协方差矩阵给定为var(ε)＝∑。注意，[6]式的模型可视为[1]式模型 的一种特殊情况。它可以推广为容纳多变量的情况(例如两个品种)， 写为
$\left(\begin{array}{c}{z}_1\\ z_2\end{array}\right)=\left(\begin{array}{cc}{x}_1& 0_1\\ 0_2& x_2\end{array}\right)\left(\begin{array}{c}{\theta }_1\\ {\theta }_2\end{array}\right)+\left(\begin{array}{c}{ϵ}_1\\ {ϵ}_2\end{array}\right),---\left[7\right]$
式中
$\mathrm{var}\left(ϵ\right)=\mathrm{var}\left(\begin{array}{c}{ϵ}_1\\ {ϵ}_2\end{array}\right)\equiv \Sigma =\left(\begin{array}{cc}{\Sigma }_{11}& {\Sigma }_{12}\\ {\Sigma }_{12}& {\Sigma }_{22}\end{array}\right),---\left[8\right]$
它仍然是[1]式模型的一种特殊情况。例如，若是第i个品种有ni 个观测结果，i＝1，2，那么[7]式中的向量zi和εi为ni维的向量，i＝1，2。 矩阵0i，i＝1，2是元素为零的矩阵。0i和Xi的行数是ni，i＝1，2，而这两 个矩阵的列数和向量θi，i＝1，2，的长度取决于模型中使用的参数数目。
在[7]式的模型中，对第i个品种预测Zi(s0)，i＝1，2时，这个过程 称为简单情况下的协克里金
$\mathrm{X\theta }=\left(\begin{array}{cc}{1}_{n_1}& 0_{n_1}\\ 0_{n_2}& 1_{n_2}\end{array}\right)\left(\begin{array}{c}{\mu }_1\\ {\mu }_2\end{array}\right),---\left[9\right]$
式中1ni和0ni分别为ni个1和0的向量，i＝1，2。对于比[7]式的模 型更一般的情况，这个过程称为通用协克里金。当联合预测多个品种 时——由向量z(s0)表示，这个过程称为多变量空间预测(Ver Hoef和 Cressie，1993)。
协方差的地质统计模型
到目前为止，我们还没有指定随机误差协方差矩阵var(ε)＝∑的 结构。确实，使本方法成为“地质统计学”的，是使用空间坐标指 定矩阵∑的结构。令一个位置的空间坐标包含在2×1的向量s中。令 第i个品种(如第一个品种)的第j个观测结果位于s＝(xs，ys)，并令第 k个品种(如第二个品种)的第m个观测结果位于u＝(xu，yu)。那么第 i个品种的第j个观测结果与第k个品种的第m个观测结果之间的协 方差可以由一个地质统计协方差模型来模拟。这种协方差模型的一个 实例是球面模型(Cressie，1993，61页)，由下式给出
$\mathrm{cov}(Z_{\mathrm{ij}},Z_{\mathrm{km}})=\left\{\begin{array}{ccc}{\phi }_{\mathrm{ik}}[1-(\frac{3\left|\right|h\left|\right|}{{2\rho }_{\mathrm{ik}}}-\frac{{\left|\right|h\left|\right|}^3}{{2\rho }_{\mathrm{ik}}^3})]& \mathrm{for}& \left|\right|h\left|\right|\le {\rho }_{\mathrm{ik}},\\ 0& \mathrm{for}& \left|\right|h\left|\right|>{\rho }_{\mathrm{ik}},\end{array}---\left[10\right]\right.$
式中h＝u-s，‖h‖是向量h的长度。本质上‖h‖是从s到u的欧 几里得距离，等于 $\sqrt{{(x_s-x_u)}^2+{(y_s-y_u)}^2}.$ 在地质统计学中，参数ik 称为门槛，而ρik称为射程。如果i＝k，[10]式就定义了自相关，否则 它就定义了互相关。从物理上说，门槛参数表示方差，而射程表示一 个距离，超过后相关性下降至零。
必须注意确保方差-协方差矩阵var(ε)＝∑是正定的。在单一品 种的情况下，使用[10]式保证了∑是正定的。例如，注意到对于自相关， 每个i的φii和ρii都大于零，足以保证正定性。然而对于互相关，要 保证正定性，对φii、φkk和φik就有所限制(Goovaerts，1977，108 页)。传统方法中使用了一类称为共区域化法的模型，以保证矩阵∑是 正定的。Ver Hoef和Barry(1998)已经给出了模拟互相关的其它方式。 注意到本实施例中∑12＝∑21＝0，当使用[10]式时，∑11和∑22中每个i的 φii和ρii都大于零，显而易见[8]式中的方差-协方差矩阵是正定的。
现在我们要介绍一类共区域化法的模型，作为本发明和本实例的 细节。对于一个指定的整数m，m＝1，2，……，M，让我们开始为第i个 变量构建一个过程如下，i＝1，……，K，
$Y_{\mathrm{im}}\left(s\right)=\sqrt{1-r_{\mathrm{im}}^2}{W}_{\mathrm{im}}\left(s\right)+r_{\mathrm{im}}{W}_{0m}\left(s\right),---\left[11\right]$
式中-1≤rim≤1，对于m＝1，2，……，M，Wim是另一个空间过程，对 于所有i≠k(i，k＝0，1，……，K)，具有以下协方差性质：
cov[Wim(s)，Wkm(s+h)]＝0，
以及对于i＝0，1，……，K，
fm(h；γm)≡cov[Wim(s)，Wim(s+h)]，
在这里，函数fm(h；γm)是某种无尺度(fm(0；γm)＝1)的模型，作 为协方差模型是有效的，也就是所有可能的协方差矩阵都是正定的。 例如，按照[10]式定义、但是没有尺度参数φik的球面模型，作为协方 差模型是有效的，如下式给出

式中γm＝ρm。对于fm(h；γm)，另一个显而易见的例子是无尺度 的“熔核效应”，对于所有h，都有fm(h；γm)＝1。熔核效应通常是 由于误差和/或可变性的采样尺度小于观测数据集中最小的试验点距 离造成的。大多数地质统计模型对于某种其它的模型都包括一种熔核 效应，比如以上[12]式中的球面模型，当h＝0时，这种效应造成原点 处的不连续。
使用[11]式和fm(h；γm)的定义，我们有
$\mathrm{cov}[Y_{\mathrm{im}}\left(s\right),Y_{\mathrm{km}}\left(s\right)]=\left\{\begin{array}{cc}{f}_m(h;{\gamma }_m)& \mathrm{ifi}=k,\\ r_{\mathrm{im}}{r}_{\mathrm{km}}{f}_m(h;{\gamma }_m)& \mathrm{ifi}\ne k.\end{array},---\left[13\right]\right.$
下面，对于i＝1，2，……，K，我们通过定义以下求和，构建共区域化 法的线性模型，
${ϵ}_i\left(s\right)=\underset{m=1}{\overset{M}{\Sigma }}{\phi }_{\mathrm{im}}{Y}_{\mathrm{im}}\left(s\right),---\left[14\right]$
注意，[14]式的构建提供了一个有效的空间随机过程。换句话说， 对于i＝1，2，……，K；m＝1，2，……，M，只要分量过程Yim(s)存在，对于 i＝1，2，……，K，过程εi(s)也就存在。还要注意，对于m＝1，2，……，M， 过程Yim(s)可能具有不同的基本协方差模型(例如球面的、指数的等 等)。
[14]式的构建为构建εi(s)提供了以下协方差模型(共区域化法的 模型)：
$C_{\mathrm{ik}}(h;\gamma )\equiv \mathrm{cov}[{ϵ}_i\left(s\right),{ϵ}_k\left(s\right)]=\left\{\begin{array}{cc}\underset{m=1}{\overset{M}{\Sigma }}{\phi }_{\mathrm{im}}^2{f}_m(h;{\gamma }_m)& \mathrm{ifi}=k,\\ \underset{m=1}{\overset{M}{\Sigma }}{r}_{\mathrm{im}}{r}_{\mathrm{im}}{\phi }_{\mathrm{im}}{\phi }_{\mathrm{km}}{f}_m(h;{\gamma }_m)& \mathrm{ifi}\ne k.\end{array}\right.---\left[15\right]$
注意，当K＝2时，参数r1m和r2m并不能各自确定，但是它们的 积是可以确定的，我们简单记为rm≡r1m r2m。现在只剩下估计以上[15] 式协方差模型中的参数，以及固定效应θ。我们在下一节讨论这些问 题。
约束最大似然：
对于[1]式中论述的混合模型，[6]式中的地质统计模型是一种特殊 情况，必须使用数据来估计模型的参数。一般说来，参数可以分为两 部分：[1]式中给出的固定效应θ，以及[8]式给出的∑中包含的协方差 参数。
让我们把∑对参数向量的依赖记为∑()。例如，可能是[12]式 中球面协方差模型的参数，或者[15]式中共区域化法模型中的参数。 一个普通的估计方法是使用约束(也称为残差)最大似然(REML) 方法估计θ和。
我们现在简要介绍REML估计方法，它是由Patterson和 Thompson(1971，1974)研究的。假设数据服从一个多变量正态分布 N(Xθ，∑())，传统的最大似然估计方法能够用于估计θ和，只要 在θ和的参数空间中，使双倍负对数似然，
L(θ，)＝Nlog(2π)+log|∑()|+(y-Xθ)′∑()-1(y-Xθ)，[16]
同时对于θ和达到最小。在的条件下，θ的一个解由下式给出

它是一个众所周知的、θ的广义最小平方估计量。把[17]式代入[16] 式，协方差参数的似然函数变为

众所周知，使[18]式中的L()最小化而得出的协方差参数的最大 似然估计是有偏的(Mardia和Marshall，1984)。REML估计是最大 似然估计的一种改进，它校正了的最大似然估计结果的偏置。对于 最小化


就获得了REML估计量，式中d为X的秩。[19]式需要一个数值 解法。关于空间协方差结构的更多REML细节可以参见Cressie(1993， 91页)。[19]式最小化的计算细节可以参见Harville(1977)，Zimmerman (1989)和Wolfinger爹人(1994)。虽然研究REML时假设数据服从一个 多变量正态分布，Heyde(1994)和Cressie和Lahiri(1996)说明了[19] 式也是一个估计方程(Godambe，1991)；所以即使数据并不服从一个多 变量正态分布，的REML估计也是无偏的。
参考了以上讨论，现在就可以详细讨论本发明的实施例了。
模型：
对于单一品种或者品种对比，为了在集合的点级别或某种区块级 别进行估计和预测，研究了一种统计模型。在品种培育的过程中，品 种的对比是最重要的目的，几乎不可能发生实验者一次只考虑一个品 种或者拟合单品种模型的情况。注意，如果对杂交品种各自的数据拟 合单品种模型，就不可能进行品种对比的统计试验。相反，使用一个 双品种或多品种的模型是最普通的情况，而且在品种对比的同时，也 能够得出对单一品种的推断(估计结果和预测结果)。所以，详细考 虑这种一般模型。此处考虑的模型是一种统计模型，既包含固定效应， 又包含随机效应，称为混合效应模型。该模型为
y＝Xθ+Dr+z+ε.                     [20]
[20]式模型中的向量y包含了所关注的主要变量(如粮食产量)， 可能还有其它伴随变量。作为这些伴随变量的一个例子，考虑生长季 节中的平均降水量。提供这些降水量数据的气象台，处于测量所关注 变量的试验点所在大区域。向y引入这些其它变量的观测结果的理由， 是使用这些变量与预测和估计中所关注的变量之间的协方差。在预测 的情况下，读者可以参考前面有关协克里金的讨论。
在一个实施例中，y可能包含恰好两个品种中一个变量的反应， 后面将要如此假设。注意，这个实施例仍然表示向量y的一般形式， 因为它能够使用两个品种的该变量反应之间的协方差。在这个实施例 中，向量y写为y＝(y1′，y2′)′，式中 $y_i^{\prime }=(Y_{i1},Y_{i2},...,Y_{{\mathrm{in}}_i}),i=\mathrm{1,2},$ 包括第一个品种的n1个观测结果和第二个品种的n2个观测结果。换句 话说，一共有来自两个品种的N≡nl+n2个观测结果。
[20]式模型中的矩阵X是固定效应的设计矩阵，
$X=\left(\begin{array}{cc}{X}_1& 0\\ 0& X_2\end{array}\right),---\left[21\right]$
式中
$X_i=\left(\begin{array}{cc}{v}_{i1}^{\prime }& x_{i1}^{\prime }\\ v_{i2}^{\prime }& x_{i2}^{\prime }\\ ·& ·\\ ·& ·\\ ·& ·\\ v_{m_i}^{\prime }& x_{{\mathrm{in}}_i}^{\prime }\end{array}\right),i=\mathrm{1,2},---\left[22\right]$
i＝1，2。在[22]式中，列向量Vij，i＝1，2，j＝1，2，……，ni包含用于模拟趋 势曲面的固定效应。这些向量包括第j个观测结果的空间坐标函数， j＝1，2，……，ni。例如，第j个观测结果的x坐标可能是是纬度，而y坐 标可能是该观测结果的经度。也可以使用其它坐标系统，比如通用横 墨卡托(UTM)。如果x坐标和y坐标记为xj和yj，i＝1，2，那么对于第 i个品种的第j个观测结果，Vij的实施例包含一个三次多项式内的项，
$v_{\mathrm{if}}={(1,x_j,y_j,x_j^2,y_j^2,x_j,y_j,x_j^3,y_j^3,x_j{y}_j^2,x_j^2{y}_j)}^{\prime }.---\left[23\right]$
在[22]式中，对于第i个品种的第j个观测结果，列向量 xij，i＝1，2，j＝1，2，……，ni包含其它固定的共变量效应。这些共变量可能包 括某些因素，比如有没有灌溉、年度的影响等等。所以固定效应参数 向量θ可以写为
$\theta =\left(\begin{array}{c}{\alpha }_1\\ {\beta }_1\\ {\alpha }_2\\ {\beta }_2\end{array}\right),---\left[24\right]$
式中子向量αi是一个ai×1的固定效应向量，用于多项式vij′ai
形成的、反映第i个品种的空间趋势曲面，βi是一个bi×1的 固定效应向量，反映第i个品种的共变量，i＝1，2。令θ的维数为d×1， 也就是d＝a1+b1+a2+b2。在[23]式给出的实施例中，对于i＝1，2，αi是 一个10×1的向量。
在[20]式的模型中，r、z和ε是随机效应的向量。向量r包含位置 的随机效应。如果有nc个位置两个品种都有观测结果，那么就只有 N-nc个不同的位置，r就变成(N-nc)×1的向量。假设E(r)＝0以及
$\mathrm{var}\left(r\right)={\sigma }_0^2{I}_r,---\left[25\right]$
式中Ir为(N-nc)×(N-nc)的单位矩阵。在[20]式的模型中，D为一 个N×(N-nc)的关联矩阵，反映位置的随机效应。关联矩阵D的行中 所有元素都置为0，每行只有一个1表明一个特定位置与该观测结果 的关联(例如，参见Searle et al.，1992，138-139页)。此外，一列中 出现两个1表明在一个位置两个品种都有。
在[20]式的模型中，z表示一个空间随机过程。向量z可能包含两 个品种的空间相关随机变量；z＝(z1′，z2′)′式中 $Z_i^{\prime }=(Z_{i1},Z_{i2},...,Z_{{\mathrm{in}}_i}).$ ，i＝1，2。假设E(z)＝0以及
$\mathrm{cov}\left(z\right)=\left(\begin{array}{cc}{\Sigma }_{11}& {\Sigma }_{12}\\ {\Sigma }_{21}& {\Sigma }_{22}\end{array}\right),---\left[26\right]$
具有某种空间自相关和互相关的变量图/协方差模式。令一个位置 的空间坐标包含在2×1的向量s中。令第i个品种的第j个观测结果 位于s＝(xs，ys)，并令第k个品种的第m个观测结果位于u＝(xu，yu)。那 么第i个品种的第j个观测结果与第k个品种的第m个观测结果之间 协方差的一个实施例是一个球面协方差模型(Cressie，1993，61页)， 由下式给出
$\mathrm{cov}(Z_{\mathrm{ij}},Z_{\mathrm{km}})=\left\{\begin{array}{ccc}{\phi }_{\mathrm{ik}}[1-(\frac{3\left|\right|h\left|\right|}{{2\rho }_{\mathrm{ik}}}-\frac{{\left|\right|h\left|\right|}^3}{{2\rho }_{\mathrm{ik}}^3})]& \mathrm{for}& \left|\right|h\left|\right|\le {\rho }_{\mathrm{ik}},\\ 0& \mathrm{for}& \left|\right|h\left|\right|>{\rho }_{\mathrm{ik}},\end{array}---\left[27\right]\right.$
式中h＝u-s，‖h‖是向量h的长度。本质上‖h‖是从s到u的欧 几里得距离， $\sqrt{{(x_s-x_u)}^2+{(y_s-y_u)}^2.}$ 在地质统计学中，参数 φik称为门槛，而ρik称为射程。如果i＝k，[10]式就定义了自相关， 否则它就定义了互相关。
必须注意确保方差-协方差矩阵[26]是正定的。传统方法中使用了 一类称为共区域化法的模型，以保证这个矩阵是正定的，但是这些方 法相当受限。注意到对于自相关，每个i的φii和ρii都大于零，足以 保证正定性。然而对于互相关，要保证正定性，对φii、φkk和φik就 有所限制(Goovaerts，1977，108页)。Ver Hoef and Barry(1998) 已经给出了模拟互相关的其它方式。注意到一个实施例中∑12＝∑21＝0， 当使用[10]式时，∑11和∑22中每个i的φii和ρii都大于零，显而易见 [26]式中的方差-协方差矩阵是正定的。
在[20]式的模型中，ε是一个包含独立随机误差的向量。向量ε可 能包含两个品种的随机变量；ε＝(ε1′，ε2′)′式中 ${ϵ}_i^{\prime }=[{ϵ}_{i1},]{ϵ}_{i2},...,{ϵ}_{{\mathrm{in}}_i]}$ ，i＝1，2。假设E(ε)＝0以及
$\mathrm{cov}\left(ϵ\right)\left(\begin{array}{cc}{\sigma }_1^2{I}_1& 0\\ 0& {\sigma }_2^2{I}_2\end{array}\right),---\left[28\right]$
式中Ii为ni×ni的单位矩阵，i＝1，2。还假设r、z和ε相互独立。
根据[20]、[25]、[26]和[28]式，以及r、z和ε相互独立的事实， 可知反应向量y的方差-协方差矩阵由下式给出
$\Sigma ={\sigma }_0^2{\mathrm{DD}}^{\prime }+\left(\begin{array}{cc}{\Sigma }_{11}& {\Sigma }_{12}\\ {\Sigma }_{21}& {\Sigma }_{22}\end{array}\right)+\left(\begin{array}{cc}{\sigma }_1^2{I}_1& 0\\ 0& {\sigma }_2^2{I}_2\end{array}\right).---\left[29\right]$
注意，矩阵DD’仅包含0和1，1在对角线上。在DD’中偏离对角 线的位置出现1表明在某些位置两个品种一起进行了试验。
干扰诊断：
在拟合[20]式的模型以便得出推断之前，很重要的是定位有干扰的 数据点以及评价它们对模型的冲击。在进行模型拟合之前，最好排除 这些数据点。注意，对于数据的空间分布进行一次目测，即可排除一 个点或一簇点，因为在空间上脱离其它观测结果。在我们的过程中， 除了对数据的目测，我们还使用干扰的以下两种目标测度，以甄别异 常或有干扰的观测结果。第一种干扰测度称为倍率。对于混合效应模 型最一般的情况，Chrestensen et al.(1992)定义了所谓的“帽子”矩 阵，
H＝X(X′∑-1X)-1X′                   [30]
作为倍率的一种测度。如果hjj＞2t/N，那么第j个数据点就被视为 一个高倍率点，式中hjj为H的第j个对角线元素，t为H的迹，N为 数据点的总数目。在一个实施例中，当∑是一个单位矩阵时，帽子矩 阵退化为以下形式(Hoaglin and Welsh，1978)：
H＝X(X′X)-1X′.                            [31]
在这种情况下，如果hjj＞2d/N，那么第j个数据点就被视为一个高 倍率点(Montgomery and Peck，1992，182页)，正如前面的定义，式中 d为θ中参数的数目。在我们的过程中使用的第二种干扰测度为对于 第j个观测结果定义的Cook距离Dj(Cook，1977)。对于混合效应模型 最一般的情况，这个距离由下式给出，
$D_j=\frac{{({\hat{\theta }}_{\left(j\right)}-\hat{\theta })}^{\prime }{X}^{\prime }{\Sigma }^{-1}X({\hat{\theta }}_{\left(j\right)}-\hat{\theta })}{d},---\left[32\right]$
式中 为排除第j个观测结果后θ的一个估计结果， ∑为估计的方差-协方差矩阵。Di＞1的观测结果通常视为有干扰的 (Montgomery and Peck，1992，182页)。在一个实施例中，当∑是一个 单位矩阵时，为了甄别有干扰的观测结果，采用以下表达式计算Cook 距离：
$D_j=\frac{{({\hat{\theta }}_{\left(j\right)}-\hat{\theta })}^{\prime }{X}^{\prime }X({\hat{\theta }}_{\left(j\right)}-\hat{\theta })}{d{\hat{\sigma }}^2},---\left[33\right]$
如前面一样，Dj＞1的观测结果视为有干扰的。在我们的过程中为 了进行干扰诊断，我们使用[31]式和[33]式，以及一个没有任何共变量 效应的模型。
为了筛选出有干扰的观测结果，既使用H又使用Dj的理由如下。 H的计算完全根据数据点的空间位置(而不是数值)。所以H对空间 的异常敏感，可以用于筛选空间异常。相反，Cook距离Dj则是部分 根据空间位置，部分根据观测值。所以，Dj可以甄别的数据点具有异 常的高低值但是不必是空间异常。
混合模型参数的估计：
介绍了[20]式给出的数据模型之后，下一步应当涉及参数估计的问 题。一般说来，这些参数可以划分为两个范畴：[24]式给出的固定效应 θ和[29]式给出的∑中包含的协方差参数。让我们把∑对于参数 的依赖记为∑()。在一个优选实施 例中，约束最大似然(REML)方法用于估计参数θ和。采用REML 估计出∑之后，利用广义最小平方估计θ，由Searle(1971，87页)给出 如下
$\hat{\theta }={\left(X^{\prime }{\Sigma }^{-1}X\right)}^{-1}{X}^{\prime }{\Sigma }^{-1}y,---\left[34\right]$
正如前面的定义，式中y为全部选定数据点的向量。
参数估计结果用于两个不同的过程：估计和预测。在两个过程中， 所关注的都是或者在单一的位置，或者是一个地理区域的平均。前者 称为点估计(预测)，后者称为区块估计(预测)。注意，点估计结 果也用于构建估计趋势曲面，也就是对于覆盖着数据的网格，在每一 个网格点上进行点估计。同样，对于覆盖着数据空间范围的网格，在 所有网格点上进行点预测，就可以构建预测曲面。
在以下两个小节中，我们将要考虑通过拟合[20]式的模型以进行估 计和预测。注意，在各个区域进行区块估计或区块预测之前，使用所 有选定的数据点(原始数据集减去有干扰的观测结果)拟合[20]式的 模型只进行了一次。由这个模型估计出的参数然后用于各个点和区域 均值的估计和预测。
区块和点估计：
假定目标是估计固定空间效应趋势曲面的均值。对于一个位置 s＝(x，y)’和第i个品种，i＝1，2，让我们定义一个空间坐标的向量值函数
$w_i\left(s\right)={(w_{i1}\left(s\right),w_{i2}\left(s\right),......w_{i{a}_i}\left(s\right))}^{\prime }.---\left[35\right]$
注意，在[23]式给出的实施例中，我们有a1＝a2＝10和 wi(s)＝(1，x，y，x2，y2，xy，x3，y3，xy2，x2y)′。所以，在某个位置s的趋势曲面值由 wi′(s)αi给出。不过，往往受到关注的是某个地理区域A的趋势曲面 平均值。如上所述，估计A上的这个平均值称为区块估计，对于i＝1，2， 区块估计结果由下式给出
${\mu }_i\left(A\right)=x_{i0}^{\prime }{\beta }_i+\frac{1}{\left|A\right|}{\int }_A{w}_i^{\prime }\left(s\right){\alpha }_i\mathrm{ds},---\left[36\right]$
式中|A|为区域A的面积，对于i＝1，2，向量x10包含一组对应于区 域A的共变量数值。假设共变量xi0与趋势曲面没有相互作用；指定 xi0的数值仅仅在相应的区域内使趋势曲面升降而不影响其形状。
在[20]式中模型的一个实施例中，因为某些理由，比如没有共变量 数据或缺少一个适当的共变量，[22]式中的共变量效应xij可以完全从 模型中消除。在这种情况下，对应于共变量的xi0′Bi，项就会 从[36]式中消失。注意，从[36]式的积分中提出xi0′βi项有助 于区分共变量与趋势曲面分量。这些共变量可能是空间的，也可能是 非空间的，可能是观测的，也可能是用户指定的。它们的数值可能是 连续的，也可能是离散的。这些共变量类型的某些实例是有规律的。 作为离散非空间共变量的一个实例，考虑年度效应，它不是空间可积 的。当一个共变量是空间的，并且数值是连续的，向量xi0中的相应分 量就表示一个积分平均结果，与[36]式中的积分相似。对位于s的第i 个品种，一个空间共变量的数值记为gi(s)，并令此共变量对应于[22] 式中向量xij的第k个分量。那么xi0的第k个分量会是 $\frac{1}{\left|A\right|}{\int }_A{g}_i\left(u\right)\mathrm{du},$ ，式中gi(u)假设是黎曼可积的。作为空间共变量观测值的一 个实例，考虑土系。在所有点上此共变量都必须已知。幸运的是所有 主要的农业区域都有土系地图，并对任何区域(例如县)都按多种土 系进行了分类。这意味着在所有可能的空间位置，都能获得作为一个 共变量的土系数据。作为用户指定的空间共变量的一个实例，考虑灌 溉，它可以作为二进制值函数(1≡灌溉的，2≡非灌溉的)对待。例 如，用户可能关注在100％灌溉情况下的一个区域中，对两个品种之 间的差异进行估计(预测)。如果作为一个共变量，灌溉对应于向量 xij的第l个分量，那么这就意味着xi0的第l个分量指定的数值为1表 示灌溉条件。
在一个实施例中有恰好两个品种的一个反应变量，对于区域A中 两个品种的区块平均，让我们将其差异定义为μ0(A)。那么μ0(A)由下 式给出，
${\mu }_0\left(A\right)\equiv {\mu }_1\left(A\right)-{\mu }_2\left(A\right)=x_{10}^{\prime }{\beta }_1-x_{20}^{\prime }{\beta }_2+\frac{1}{\left|A\right|}{\int }_A[w_1^{\prime }\left(s\right){\alpha }_1-w_2^{\prime }\left(s\right){\alpha }_2]\mathrm{ds}.---\left[37\right]$
现在，让我们定义一个集合SA＝{p：sp∈A}，其中的点{sp}或者是位 于覆盖在A上的规则网格上，或者是在A之内随机选取的，SA中有 nA个元素。让我们再定义
${\lambda }_1\left(s_p\right)\equiv \left(\begin{array}{c}{w}_1\left(s_p\right)\\ x_{10}\\ 0\\ 0\end{array}\right),{\lambda }_2\left(s_p\right)\equiv \left(\begin{array}{c}0\\ 0\\ w_2\left(s_p\right)\\ x_{20}\end{array}\right),\mathrm{and}{\lambda }_0\left(s_p\right)\equiv \left(\begin{array}{c}{w}_1\left(s_p\right)\\ x_{10}\\ -w_2\left(s_p\right)\\ -x_{20}\end{array}\right).---\left[38\right]$
那么，近似[36]式和[37]式中的积分，并使用[24]式中θ的定义， 我们有
${\mu }_i\left(A\right)\approx x_{10}^{\prime }{\beta }_i+\frac{1}{n_A}\underset{{p\in S}_A}{\Sigma }{w}_i^{\prime }\left(s_p\right){\alpha }_1=\frac{1}{n_A}\underset{{p\in S}_A}{\Sigma }{\left[{\lambda }_i\left(s_p\right)\right]}^{\prime }\theta ,---\left[39\right]$
i＝1，2，以及
${\mu }_0\left(A\right)\approx x_{10}^{\prime }{\beta }_1-x_{20}^{\prime }{\beta }_2+\frac{1}{n_A}\underset{{p\in S}_A}{\Sigma }[w_1^{\prime }\left(s_p\right){\alpha }_1-w_2^{\prime }\left(s_p\right){\alpha }_2]=\frac{1}{n_A}\underset{{p\in S}_A}{\Sigma }{\left[{\lambda }_0\left(s_p\right)\right]}^{\prime }\theta .---\left[40\right]$
在[39]式和[40]式中代入[34]式中θ的广义最小平方估计量，我们 就得到了μi(A)的估计量
${\hat{\mu }}_1\left(A\right)\approx \frac{1}{n_A}\underset{{p\in S}_A}{\Sigma }{\left[{\lambda }_1\left(s_p\right)\right]}^{\prime }{\left(X^{\prime }{\Sigma }^{-1}X\right)}^{-1}{X}^{\prime }{\Sigma }^{-1}y,---\left[41\right]$
i＝0，1，2。 的方差为
$\mathrm{var}\left({\hat{\mu }}_1\left(A\right)\right)\approx \frac{1}{n_A^2}\underset{{p\in S}_A}{\Sigma }\underset{{q\in S}_A}{\Sigma }{\left[{\lambda }_i\left(s_p\right)\right]}^{\prime }{\left(X^{\prime }{\Sigma }^{-1}X\right)}^{-1}\left[{\lambda }_i\left(s_q\right)\right],---\left[42\right]$
其标准差为 $\mathrm{se}\left({\hat{\mu }}_i\left(A\right)\right)=\sqrt{\mathrm{var}\left[{\hat{\mu }}_i\left(A\right)\right]}$ i＝0，1，2。
现在让我们考虑点估计的情况。目的是估计以下函数，它是[36] 式的空间情况，其中区域A现在退化为一个位置(点)。在s0的点估 计为
μi(s0)＝xi0′βi+wi′(s0)αi.          [43]
正如区块估计时的情况，共变量xi0的数值仅仅通过使点估计结果 或增或减，在[43]式中引入了一个附加的效应。然而，由于这种附加效 应的幅度和方向在不同的点可能有变化，模型中有了这些共变量之后， 采用点估计结果构建的整个趋势曲面的形状可能发生实质性的变化。
作为[41]式中给出的区块估计量的一种特殊情况，μi(s0)的估计量 可以由下式获得
${\hat{\mu }}_i\left(s_0\right)={\left[{\lambda }_i\left(s_0\right)\right]}^{\prime }{\left(X^{\prime }{\Sigma }^{-1}X\right)}^{-1}{X}^{\prime }{\Sigma }^{-1}y,---\left[44\right]$
i＝0，1，2。类似于[42]式中的区块估计结果，可以获得 标的估计方差，由下式给出，
$\mathrm{var}\left[{\hat{\mu }}_i\left(s_0\right)\right]={\left[{\lambda }_i\left(s_0\right)\right]}^{\prime }{\left(X^{\prime }{\Sigma }^{-1}X\right)}^{-1}\left[{\lambda }_i\left(s_0\right)\right].---\left[45\right]$
点估计结果的标准差为 $\mathrm{se}\left[{\hat{\mu }}_i\left(s_0\right)\right]=\sqrt{\mathrm{var}\left[{\hat{\mu }}_i\left(s_0\right)\right]}$ i＝0，1，2。
通用区块和点协克里金：
让我们首先考虑区块预测，它也称为区块协克里金，正如在本发 明的背景技术中的讨论。这里的目标是在某个指定的地理区域，估计 反应曲面的平均值，它是以下三个分量的和：固定的空间效应的趋势 曲面、固定的协方差效应和自相关的随机曲面。定义随机函数Zi(s)， 其中集合{Zi(s)；s∈A}；i＝1，2形成了一个随机曲面，可以进行 Cressie(1993，107页)定义的积分。令Z0(s)≡Z1(s)-Z2(s)为一个随机函 数，其中{Z0(s)；s∈A}形成了一个随机曲面，如果Z1(s)和Z2(s)可积， 它就可积。现在，对于指定的、xi0，i＝1，2给出的共变量数值集合，让 我们定义函数
$S_i\left(A\right)\equiv x_{i0}^{\prime }{\beta }_i+\frac{1}{\left|A\right|}{\int }_A[w_i^{\prime }\left(s\right){\alpha }_i+Z_i\left(s\right)]\mathrm{ds},---\left[46\right]$
i＝1，2，为第i个品种在区域A上平均实现值。正如区块估计时的 情况(参见[36]式后的讨论)，[46]式中的xi0′βi项可以容纳 所有类型的共变量，并表示对应于共变量xi0的一个平均的结果。
这里的目的是根据数据的某种线性组合，比如说ωi′y 预测区块均值S1(A)、S2(A)和S0(A)≡S1(A)-S2(A)。向量ωi，i＝0，1，2 的估计是由一种称为最优线性无偏预测的方法完成的，预测量是数据 的一个线性函数，是无偏的，并且在所有线性无偏预测量中具有最小 的均方预测误差。
让我们首先假设[20]式给出的线性模型。为了预测区块均值Si(A)， i＝0，1，2，通过最小化
E[Si(A)-ωi′y]2，                 [47]
可以获得最优线性无偏预测量(BLUP)，此时需要无偏条件
E[Si(A)]＝E(ωi′y)，              [48]
i＝0，1，2。
现在，使用[20]式的模型和[46]式，i＝1，2时，[48]式可以写为
$\frac{1}{\left|A\right|}{\int }_A{w}_i^{\prime }\left(s\right){\alpha }_i\mathrm{ds}+x_{0i}^{\prime }{\beta }_i={\omega }_i^{\prime }\mathrm{X\theta },---\left[49\right]$
而i＝0时，[48]式可以写为
$\frac{1}{\left|A\right|}{\int }_A{w}_1^{\prime }\left(s\right){\alpha }_1\mathrm{ds}+x_{01}^{\prime }{\beta }_1-\frac{1}{\left|A\right|}{\int }_A{w}_2^{\prime }\left(s\right){\alpha }_2\mathrm{ds}-x_{02}^{\prime }{\beta }_2={\omega }_0^{\prime }\mathrm{X\theta }.---\left[50\right]$
一般说来，使用[35]式中定义的向量值函数wi(s)，i＝0，1，2，[49]式和 [50]式可以结合成单一的方程
ui′θ＝ωi′Xθ，                         [51]
i＝0，1，2，式中
$u_1\equiv \left(\begin{array}{c}\frac{1}{\left|A\right|}{\int }_A{w}_{11}\left(s\right)\mathrm{ds}\\ .\\ .\\ .\\ \frac{1}{\left|A\right|}{\int }_A{w}_{{1a}_1}\left(s\right)\mathrm{ds}\\ x_{01}\\ 0\\ .\\ .\\ .\\ 0\\ 0\end{array}\right),u_2\equiv \left(\begin{array}{c}0\\ .\\ .\\ .\\ 0\\ 0\\ \frac{1}{\left|A\right|}{\int }_A{w}_{21}\left(s\right)\mathrm{ds}\\ .\\ .\\ .\\ \frac{1}{\left|A\right|}{\int }_A{w}_{{2\mathrm{a\alpha }}_2}\left(s\right)\mathrm{ds}\\ x_{02}\end{array}\right),\mathrm{and}{u}_0\equiv \left(\begin{array}{c}\frac{1}{\left|A\right|}{\int }_A{w}_{11}\left(s\right)\mathrm{ds}\\ .\\ .\\ .\\ \frac{1}{\left|A\right|}{\int }_A{w}_{{1a}_1}\left(s\right)\mathrm{ds}\\ x_{01}\\ -\frac{1}{\left|A\right|}{\int }_A{w}_{21}\left(s\right)\mathrm{ds}\\ .\\ .\\ .\\ -\frac{1}{\left|A\right|}{\int }_A{w}_{{2a}_2}\left(s\right)\mathrm{ds}\\ -x_{02}\end{array}\right).---\left[52\right]$
为了满足无偏的条件，对于每一个可能的θ，[51]式都必须成立， 因此有
ui′＝ωi′X，                             [53]
i＝0，1，2。
[52]式中的向量可以近似为
$u_1\approx \left(\begin{array}{c}\frac{1}{n_A}\underset{{p\in S}_A}{\Sigma }{w}_{11}\left(s_p\right)\\ .\\ .\\ .\\ \frac{1}{n_A}\underset{{p\in S}_A}{\Sigma }{w}_{a11}\left(s_p\right)\\ x_{01}\\ 0\\ .\\ .\\ .\\ 0\\ 0\end{array}\right),u_2\approx \begin{array}{c}\left(\begin{array}{c}0\\ .\\ .\\ .\\ 0\\ 0\\ \frac{1}{n_A}\underset{{p\in S}_A}{\Sigma }{w}_{21}\\ .\\ .\\ .\\ \frac{1}{n_A}\underset{{p\in S}_A}{\Sigma }{w}_{2a_2}\left(s_p\right)\\ x_{02}\end{array}\left(s_p\right)\right)\end{array},\mathrm{and}{u}_0\approx \left(\begin{array}{c}\frac{1}{n_A}\underset{{p\in S}_A}{\Sigma }{w}_{11}\left(s_p\right)\\ .\\ .\\ .\\ \frac{1}{n_A}\underset{{p\in S}_A}{\Sigma }{w}_{{1\mathrm{a\alpha }}_1}\left(s_p\right)\\ x_{01}\\ -\frac{1}{n_A}\underset{{p\in S}_A}{\Sigma }{w}_{21}\left(s_p\right)\\ .\\ .\\ .\\ -\frac{1}{n_A}\underset{{p\in S}_A}{\Sigma }{w}_{{2\mathrm{a\alpha }}_2}\left(s_p\right)\\ -x_{02}\end{array}\right),---\left[54\right]$
式中SA已经在前面定义过。
现在，使用[20]式、[46]式和[52]式，均方预测误差[47]式可以写为
$E{[S_i\left(A\right)-{\omega }_i^{\prime }y]}^2=E{[u_i^{\prime }\theta +\frac{1}{\left|A\right|}{\int }_A{Z}_i\left(s\right)\mathrm{ds}-{\omega }_i^{\prime }(\mathrm{X\theta }+\mathrm{Dr}+z+ϵ)]}^2---\left[55\right]$
i＝0，1，2。
让我们定义N向量 $\delta =({\delta }_{11},{\delta }_{12},......,{\delta }_{1n_1},{\delta }_{21},{\delta }_{22},.....{\delta }_{2n_2})$ δ＝Dr+z+ε，并 令 $Z_i\left(A\right)\equiv \frac{1}{\left|A\right|}{\int }_A{Z}_i\left(s\right)\mathrm{ds},$ i＝0，1，2。然后，使用[51]式，[55]式可 以写为
E[Si(A)-ωi′y]2＝E[Zi(A)-ωi′δ]2，         [56]
i＝0，1，2。
求取[56]式中的数学期望需要求取cov[Zi(A)，δ]，它可以近似为
$\left(\begin{array}{c}\frac{1}{n_A}\underset{{p\in S}_A}{\Sigma }\mathrm{cov}[Z_i\left(s_p\right),{\delta }_{11}]\\ .\\ .\\ .\\ \frac{1}{n_A}\underset{{p\in S}_A}{\Sigma }\mathrm{cov}[Z_i\left(s_p\right),{\delta }_{{1n}_1}]\\ \frac{1}{n_A}\underset{{p\in S}_A}{\Sigma }\mathrm{cov}[Z_i\left(s_p\right),{\delta }_{21}]\\ .\\ .\\ .\\ \frac{1}{n_A}\underset{{p\in S}_A}{\Sigma }\mathrm{cov}[Z_i\left(s_p\right),{\delta }_{{2n}_2}]\end{array}\right),---\left[57\right]$
i＝0，1，2。
现在，由于cov(z，ε)＝0和cov(z，r)＝0，[57]式中的向量可以写为
$c_1\equiv \mathrm{cov}[Z_1\left(A\right),\delta ]\approx \left(\begin{array}{c}\frac{1}{n_A}\underset{{p\in S}_A}{\Sigma }\mathrm{cov}[Z_i\left(s_p\right),Z_1\left(s_1\right)]\\ .\\ .\\ .\\ \frac{1}{n_A}\underset{{p\in S}_A}{\Sigma }\mathrm{cov}[Z_i\left(s_p\right),Z_1\left(s_{n_1}\right)]\\ \frac{1}{n_A}\underset{{p\in S}_A}{\Sigma }\mathrm{cov}[Z_i\left(s_p\right),Z_2\left(s_1\right)]\\ .\\ .\\ .\\ \frac{1}{n_A}\underset{{p\in S}_A}{\Sigma }\mathrm{cov}[Z_i\left(s_p\right),Z_2\left(s_{n_2}\right)]\end{array}\right),---\left[58\right]$
i＝0，1，2
求取[56]式中的数学期望还需要求取Ci(A，A)≡var[Zi(A)]，它可以 近似为
$C_i(A,A)\approx \frac{1}{n_A^2}\underset{{p\in S}_A}{\Sigma }\underset{{q\in S}_A}{\Sigma }\mathrm{cov}[Z_i\left(s_p\right),Z_i\left(s_q\right)],---\left[59\right]$
i＝0，1，2。现在，由于 $Z_i\left(A\right)\equiv \frac{1}{\left|A\right|}{\int }_A{Z}_i\left(s\right)\mathrm{ds},$ 它可以近似为 $\frac{1}{n_A}\underset{{p\in S}_A}{\Sigma }{Z}_i\left(s_p\right),$ 我们注意到
$E\left[Z_i^2\left(A\right)\right]=E{\left[\frac{1}{n_A}\underset{{p\in S}_A}{\Sigma }{Z}_i\left(s_p\right)\right]}^2$
$=\frac{1}{n_A^2}\underset{{p\in S}_A}{\Sigma }\underset{{q\in S}_A}{\Sigma }E\left[Z_i\left(s_p\right)Z\left(s_q\right)\right]$
$\frac{1}{n_A^2}\underset{{p\in S}_A}{\Sigma }\underset{{q\in S}_A}{\Sigma }\mathrm{cov}[Z_i\left(s_p\right),Z_i\left(s_q\right)]$
$=C_i(A,A)$
因为对于任何p∈SA以及i＝0，1，2，E[Zi(sp)]＝0。所以，[56]式中的 均方预测误差可以写为
$E{[Z_i\left(A\right)-{\omega }_i^{\prime }\delta ]}^2=E[Z_i^2\left(A\right)-2Z_i\left(A\right){\omega }_i^{\prime }\delta +{\omega }_i^{\prime }\delta {\delta }^{\prime }{\omega }_i]$
$=E\left[Z_i^2\left(A\right)\right]-2{\omega }_i^{\prime }E\left[Z_i\left(A\right)\delta \right]+{\omega }_i^{\prime }E\left(\delta {\delta }^{\prime }\right){\omega }_i$
$=C_i(A,A)-2{\omega }_i^{\prime }{c}_i+{\omega }_i^{\prime }\Sigma {\omega }_i,$
[60]
i＝0，1，2。现在，为了满足[53]式中给出的无偏条件，增加一个拉格 朗日乘子mi，[60]式可以写为
L(ωi，mi)≡Ci(A，A)-2ωi′ci+ωi′∑ωi-2(ui′-ωi′X)mi，   [61]
i＝0，1，2。[61]式对ωi和mi求偏导，有
$\frac{\partial L({\omega }_i,m_i)}{\partial {\omega }_i^{\prime }}=-2c_i+2\Sigma {\omega }_i+2X{m}_i,---\left[62\right]$
以及
$\frac{\partial L({\omega }_i,m_i)}{\partial m_i^{\prime }}=-2(u_i-X^{\prime }{\omega }_i),---\left[63\right]$
i＝0，1，2。
[62]式和[63]式中的这两个偏导数置为0，然后联立求解，使[60] 式中的均方预测误差在[53]式中给出的无偏条件下最小化，并产生区 域A中品种性能(或性能差异)的最优线性无偏预测量(BLUP)如下。 区块预测(协克里金)方程为
$\left(\begin{array}{cc}\Sigma & X\\ X^{\prime }& 0\end{array}\right)\left(\begin{array}{c}{\omega }_i\\ m_i\end{array}\right)=\left(\begin{array}{c}{c}_i\\ u_i\end{array}\right),---\left[64\right]$
i＝0，1，2。对方程[64]求解ωi由下式给出，
    ∑ωi+Xmi＝ci
    ωi＝∑-1(ci-Xmi)，                          [65]
i＝0，1，2。所以，区域A中品种性能(或性能差异)的预测量由下 式给出，
${\hat{S}}_i\left(A\right)={\omega }_i^{\prime }y,---\left[66\right]$
i＝0，1，2。把从方程[65]解出的ωi代入[60]式的表达式中，就提供了 均方预测误差或称为预测方差的表达式如下，
$\mathrm{var}\left[{\hat{S}}_i\left(A\right)\right]=C_i(A,A)-2{\omega }_i^{\prime }{c}_i+{\omega }_i^{\prime }\Sigma \left[{\Sigma }^{-1}(c_i-X_{m_i})\right]$
$=C_i(A,A)-2{\omega }_i^{\prime }{c}_i+{\omega }_i^{\prime }{c}_i-{\omega }_i^{\prime }X{m}_i$
$=C_i(A,A)-{\omega }_i^{\prime }{c}_i-{\omega }_i^{\prime }X{m}_i$
$=C_i(A,A)-{\omega }_i^{\prime }(c_i+X{m}_i),$
[67]
i＝0，1，2。把从方程[65]解出的ωi代入从矩阵方程[64]式下半部获得 的方程X’ωi＝ui中，就提供了拉格朗日乘子mi的表达式如下，
X′ωi＝ui
X′∑-1(ci-Xmi)＝ui
X′-1ci-X′∑-1Xmi＝ui                [68]
mi＝(X′∑-1X)-1(X′∑-1ci-ui)，
i＝0，1，2。
我们现在考虑点预测或称为点协克里金的情况。注意，方程[46] 的一种特殊情况是区域A是单一的点，比如说s0，也就是， Si(s0)≡xi0′βi+wi′(S0)αi+Zi(s0)，i＝1，2。现在的目标是根据数据的某种线性组合， 比如说ωi*′y，在点s0进行点预测，也就是，预测S1(s0)、 S2(s0)和S0(s0)≡S1(s0)-S2(s0)。换句话说，点预测量由下式给出，
${\hat{S}}_i\left(s_0\right)={{\omega }_i^{*}}^{\prime }y,---\left[69\right]$
i＝0，1，2，式中协克里金权向量，ωi*需要确定。与方程[65] 类似，权向量ωi*由下式给出，
${\omega }_i^{*}={\Sigma }^{-1}(c_i^{*}-X{m}_i^{*}),---\left[70\right]$
式中向量ci*类似于为了[58]式中的区块预测而定义的 向量ci，拉格朗日乘子向量mi*类似于为了[68]式中的区块 预测而定义的向量mi。所以向量mi*由下 式给出，
$c_i^{*}\equiv \mathrm{cov}[Z_i\left(s_0\right),\delta ]\approx \left(\begin{array}{c}\mathrm{cov}[Z_i\left(s_0\right),Z_1\left(s_1\right)]\\ .\\ .\\ .\\ \mathrm{cov}[Z_i\left(s_0\right),Z_1\left(s_{n_1}\right)]\\ \mathrm{cov}[Z_i\left(s_0\right),Z_2\left(s_1\right)]\\ .\\ .\\ .\\ \mathrm{cov}[Z_i\left(s_0\right),Z_2\left(s_{n_2}\right)]\end{array}\right),---\left[71\right]$
以及
$m_i^{*}={\left(X^{\prime }{\Sigma }^{-1}X\right)}^{-1}(X^{\prime }{\Sigma }^{-1}{c}_i^{*}-u_i^{*}),---\left[72\right]$
i＝0，1，2。在方程[72]中，向量ui*，i＝0，1，2类似于为了 区块预测的情况而在方程[52]中给出的向量ui，i＝0，1，2，现在定义为
$u_1^{*}\equiv \left(\begin{array}{c}{w}_{11}\left(s\right)\\ .\\ .\\ .\\ w_{{1\mathrm{a\alpha }}_1}\left(s\right)\\ x_{01}\\ 0\\ .\\ .\\ .\\ 0\\ 0\end{array}\right),u_2^{*}\equiv \left(\begin{array}{c}0\\ .\\ .\\ .\\ 0\\ 0\\ w_{21}\left(s\right)\\ .\\ .\\ .\\ w_{{2\mathrm{a\alpha }}_2}\left(s\right)\\ x_{02}\end{array}\right),\mathrm{and}{u}_0^{*}\equiv \begin{array}{c}\left(\begin{array}{c}{w}_{11}\left(s\right)\\ .\\ .\\ .\\ w_{{1\alpha }_1}\left(s\right)\\ x_{01}\\ -w_{21}\left(s\right)\\ .\\ .\\ .\\ -w_{{2a}_2}\left(s\right)\\ -x_{02}\end{array}\right).---\left[73\right]\end{array}$
最后，正如在方程[67]中一样，获得点预测量 的预 测方差为
$\mathrm{var}\left[{\hat{S}}_i\left(s_0\right)\right]=C_i(s_0,s_0)-{{\omega }_i^{*}}^{\prime }(c_i^{*}+X{m}_i^{*}).---\left[74\right]$
以下包括的实例是为了说明本发明的特定实施例。本领域的技术 人员应当理解，在以下的实例中公开的技术是发明人发现的，在本发 明的实践中表现良好，因此可视为构成了其实践的特定模式。不过， 利用本公开文件的利益，本领域的技术人员将会理解，在这些公开的 特定实施例中可以进行许多改变，依然可以获得相象或相似的结果， 而不脱离本发明的实质和范围。阅读以下实例时应当参考图1-图12， 在以上的附图简要说明中已有小结。
实例一
两个Zea mays品种在广阔区域试验中的产量对比
我们首先说明本发明如何用于两个Zea mays品种在广阔区域的产 量对比。这两个品种称为DK512和DK527，使用从1996年到1998 年三年的广阔区域试验数据进行对比。
本过程的第一步包括使用帽子矩阵[31]式和Cook距离[33]式整理 数据。对位置也进行了目测，以排除任何异常的空间簇。整理之后， 对于DK512和DK527这两个杂交品种，跨年度的样本数目分别为 n1＝445，n2＝389。两个杂交品种的数据合并之后，样本数目为N＝834。 按年度的样本数目由下表给出。
表1   年度   杂交品种   样本数目   1996   DK512   141   113   1997   DK512   150   155   1998   DK512   154   121
本过程的下一步包括拟合[20]式中给出的混合模型。对于两个杂交 品种，[20]式中的固定效应Xθ都是由[23]式中给出，但是没有截距项 1。因此，[24]式中向量θ的β1和β2分量，每个的维数都是9×1。在 向量θ中，α1和α2分量是3×1的子向量，表示1996年到1998年三 年中每一年的截距参数。所以，[22]式中的Xii＝1，2为ni’12，[21]式中 的X为N’24矩阵。对于[26]式中的方差协方差矩阵，我们使用[15]式 中给出的共区域化法模型，而且M＝2，其中f1(h；γm)和f2(h；γm)都是 [12]式中给出的球面模型。[28]式中给出了对每个变量带有一种单独的 熔核效应的模型。使用的共区域化法模型带有熔核效应，以及如方程 [25]之前的段落所介绍的一种位置匹配的随机效应，所有数据的方差 -协方差模型由[29]式给出。
下一步是估计固定效应参数和协方差参数。使[19]式中给出的 REML方程最小化，可以估计方差一协方差参数。估计出的参数由下 表给出。
表2   参数   估计结果   参考方程   σ12   59.6   [28]，[29]   σ22   57.6   [28]，[29]   φ112   144.6   [15]，[29]   φ212   108.3   [15]，[29]   r1   0.941   [15]，[29]   ρ1   0.050   [12]，[15]，[29]   φ122   0.307   [15]，[29]   φ222   36.1   [15]，[29]   r2   0.999   [15]，[29]   ρ2   0.012   [12]，[15]，[29]   σ02   367.6   [25]，[29]
获得了参数估计结果之后，使用该模型进行两个杂交品种之间产 量差异的空间预测和空间估计。对于单独的地理位置(点预测和估计) 和地理区域(区块预测和估计)都要进行产量差异的空间预测和空间 估计。除了这些空间计算之外，还使用传统的简单均值法估计地理区 域的平均产量差异。(在这个实例中选用的地理区域表示推销商业化 品种所使用的区域。)
回顾一下，估计的目标是确定两个品种产量差异的长期平均。预 测的目标是在特定的年度和/或在特定的位置确定两个品种的平均产 量差异。因为这个理由，空间估计是以跨年度的方式进行的，而空间 预测是以按年度的方式进行的。为了提供与传统均值方法对比的基础， 对按年度和跨年度的分析，都进行了简单均值分析。(跨年度简单均 值方法是基于所有三年的汇总数据。)
使用表2中的参数估计结果和[38]式中的λ0(sp)，以及 x10＝x0＝(1/3，1/3，1/3)￠，从[41]式获得市场区的跨年度平均产量差异的 区块估计结果。(简单地设置x10＝x0＝(1/3，1/3，1/3)￠，在估计中为所有 三个年度提供了相同的加权系数。)使用[41]式获得区块估计结果时， 为每个市场区{sp}∈A产生一个致密、规则的点网格，其中A是市场区。 另外，从[41]式中估计出的数值，以其标准差——[42]式的平方根—— 划分时，假设具有标准正态分布；称这为z统计。如果z统计小于-1.96 或者大于1.96，对于两个品种之间无差异的虚假设，就假定其显著程 度在α＝0.05的级别。
使用表2中的参数估计结果和[54]式中的u0，从[66]式获得多个市 场区按年度平均产量差异的区块预测。在使用[54]式时，对于1996、 1997和1998年，x10和x20都分别设置为(1，0，0)￠、(0，1，0)￠和(0，0，1)￠。 对于每个市场区，在使用致密、规则的点网格方面，以及为了确定统 计显著程度而计算和使用z统计方面，区块预测过程都类似于区块估 计过程。在区块预测中计算z统计时，从[66]式中预测出的数值，以其 标准差——[67]式的平方根——划分时，假设具有标准正态分布。
下一步我们将考虑两个品种之间产量差异的点估计和点预测。点 估计和点预测都是在25×25的矩形网格上进行的，该网格覆盖着采集 每个杂交品种数据的区域。在25×25网格上的点估计结果和点预测结 果共同形成了产量差异的趋势曲面，它可以用于研究杂交品种试验的 地理区域中产量差异的变化。
在以上讨论的25×25网格上，跨年度的点估计使用[44]式计算产 量差异μ0(s0)。为了确定统计显著程度，已经发现，拒绝一个真实的、 两个品种之间无差异的虚假设的概率为
$2\{1-\Phi \left(\right|{\hat{\mu }}_0\left(s_0\right)/\sqrt{\mathrm{var}\left[{\hat{\mu }}_0\left(s_0\right)\right]}\left|\right)\},$
式中 由[45]式给出，Φ(·)为一个标准正态随机 变量的累积分布函数。
在以上讨论的25×25网格上，按年度的点预测使用[69]式计算产 量差异S0(s0)。为了确定统计显著程度，已经发现，拒绝一个真实的、 两个品种之间无差异的虚假设的概率为
$2\{1-\Phi \left(\right|{\hat{S}}_0\left(s_0\right)/\sqrt{\mathrm{var}\left[{\hat{S}}_0\left(s_0\right)\right]}\left|\right)\},$
式中 由[74]式给出，如同前面的情况，Φ(·)为 一个标准正态随机变量的累积分布函数。
以上分析的全部结果以图形方式展示在图1至图12中。这些图分 为两种类型，一种是等值线图，另一种是等值区域图。等值区域图是 一种二维彩色地图，表示一个变量数值在一个区域的变化。
在等值线图中，三个级别(-5、0和+5蒲式耳/英亩)的产量差异 和四个级别(产量差异±5％和±20％)的统计显著程度都用等值线表 示。在等值线图中还显示了两个杂交品种的试验点和市场区的边界和 名称。
在等值线图中，实等值线表示产量差异，虚等值线表示显著程度 级别。如同任何等值线图，各条线的关系表示同一变量(如产量差异) 的不同级别，用于确定曲面的大体形态。例如，对于产量差异，一个 跨越红色、蓝色和绿色实线的横切面(以此顺序)，就是在产量差异 曲面上向上移动，分别是从-5至0至+5蒲式耳/英亩。
对于产量差异和显著程度级别而言，等值线之间的关系用于在产 量差异和统计显著程度的多种级别，确定产量差异发生显著变化的地 理区域。例如，一条绿虚线围绕的区域落在一条绿实线围绕的区域中， 后者又落在一条蓝实线围绕的区域中，表明该区域具有大于+5蒲式耳 /英亩的、统计意义上显著的产量差异(在显著程度级别0.05)。
另一种类型的图形，也就是等值区域图，对于若干地理区域(市 场区)中的每一个，显示这两个杂交品种在其中的平均产量差异。等 值区域图上的图例既表明产量差异的级别(小于-5、-5和+5之间以及 大于+5蒲式耳/英亩)，又表明差异是否达到5％级别的统计显著程度。
对于空间预测方法(仅仅按年度，1996年到1998年)和空间估 计方法(仅仅跨年度)产生等值区域图和等值线图。等值线图不适合 简单均值方法，所以只产生了等值区域图(既有按年度也有跨年度)。 为了保持一致，如果某个市场区没有一个试验点种植了两个杂交品种， 该市场区的产量差异就不会展示在等值区域图上(回顾一下，简单均 值方法需要位置匹配的数据。)
为了对比传统的简单均值方法和空间方法，简单均值方法的按年 度地图可以与按年度的区块预测地图对比，简单均值方法的跨年度地 图可以与区块估计方法的跨年度地图对比。
实例二
有效性研究
为了表明本发明在农作物品种性能评价的应用条件下的功效，进 行了几项有效性研究。对于本发明的两个组件：估计和预测，进行了 分组的有效性研究。回顾以上讨论的实例，这两个组件中的每一个又 包括两个部分：(1)点估计(或预测)部分，用于构建带有等值线的估 计(或预测)曲面地图；(2)区块估计(或区块预测)部分，在实例中 通过映射每个区域的平均值，用于构建市场区性能地图。
我们以对预测的有效性研究开始。首先，我们讨论点预测或点克 里金的有效性研究。这项研究的基础是互证有效(Cressie，1993，101 页)，使用Monsanto玉米研究系统对于两个DEKALB玉米杂交品种， DK512和DK527，从1996年到1998年的实际数据。对于两个杂交品 种估计模型的预测质量，互证有效用作一种诊断检验。在互证有效的 过程中，基本的思路是逐个地每次删除一个数据，然后使用其它数据 预测已排除的数据。对于估计和预测的地质统计学方法，这是一种常 见的评价方式。在进行互证有效时，834个数据中的每一个，都是一 次排除一个，然后使用我们的过程预测已排除的数据。注意，这里我 们将要预测随机变量Yi(s0)；i＝1，2，它是从数据集中排除的观测值之一。 我们不想预测光滑值Si(s0)；i＝1，2，它具有与排除的熔核效应相关联的 误差项。我们使用了协克里金预测量，使用Cressie(1993，138页)介绍 的协方差。
在互证有效的研究中，计算了以下五种统计量。
1.均方预测误差(MSPE)，以量 ${[Y_i\left(s_0\right)-{\hat{Y}}_i\left(s_0\right)]}^2$ 的平均值计 算，其中 是已排除数据的预测值，而Yi(s0)是观测值。
2.平均估计预测方差(Cressie，1993，140页，3.2.52式)。
3.80％预测区间包含真实值次数的比例。
4.标准化偏置： $[Y_i\left(s_0\right)-{\hat{Y}}_i\left(s_0\right)]/\mathrm{se}\left({\hat{Y}}_i\left(s_0\right)\right),$ ，(Cressie， 1993，102页，2.6.15式)。
5.标准化MSPE，以量 ${[Y_i\left(s_0\right)-{\hat{Y}}_i\left(s_0\right)]}^2/\mathrm{var}\left({\hat{Y}}_i\left(s_0\right)\right).$ 的平均值的平方根计算。
在构建预测区间时，预测区间的下限计算为 ${\hat{Y}}_i\left(s_0\right)-1.28*\mathrm{se}\left({\hat{Y}}_i\left(s_0\right)\right),$ ${\hat{Y}}_1\left(s_0\right)+1.28*\mathrm{se}\left({\hat{Y}}_i\left(s_0\right)\right).,$ 上限计算为 这里 原文53页19行2是预测标准差，它计算为估计预测方差(Cressie，1993， 140页，3.2.52式)的平方根。
跨越所有三个年度归纳了以上的统计量。我们现在按照支持本发 明的需要和实际服从的规律，展示互证有效研究的结果。
a.MSPE(统计量1)和平均估计预测方差(统计量2)应当相互 接近。这两个统计量的互证有效研究结果分别为269.1和297.9。
b.预测区间(统计量3)应当接近名义(80％)覆盖。研究中的 平均覆盖为82.7％。
c.标准化偏置(统计量4)应当接近0。研究中的标准化偏置计算 结果为-0.021。
d.标准化MSPE(统计量5)应当接近1。研究中的标准化MSPE 为1.049。
我们可以看到，互证有效的结果非常接近理论上需要的值。换句 话说，该有效性结果支持以下说法：作为点预测工具，在[20]式中提 出的模型符合我们的期望。
对于本发明中区块克里金组件的有效性研究，也是基于玉米杂交 品种产量的实际数据分析。这种有效性研究的基础是把这个组件与传 统的方法对比。回顾一下，在传统的方法中，来自一个区域的所有位 置匹配的、并排种植的品种差异数据的简单均值，用于预测该区域的 均值。我们将称这种方法为均值方法。有效性研究包括的步骤如下：
1.选择在国内有充足重叠数据的一对玉米杂交品种。
2.对于这对杂交品种，搜集其全部位置匹配的、并排种植的差异 数据。这将被称为“原始种群数据”。
3.在原始种群数据的空间范围之内构建一个矩形。我们将称这个 矩形为“预测矩形”，我们将在其中构建平均预测值。
4.使用预测矩形内部的所有原始种群数据的均值表示该矩形的 “真实均值”。
5.随机选择原始种群数据的50％。这表示可用于分析的“样本 数据”。
6.使用落在预测矩形内部的样本数据的均值表示用均值方法预 测的、该矩形的均值。
7.使用(预测矩形内外的)全部样本数据，应用本发明的通用区 块克里金组件预测该预测矩形的均值。这表示区块预测的数值。
8.在步骤4的真实均值与步骤6和步骤7的两个预测均值之中的 每一个之间，计算平方差异。
9.对于几对杂交品种重复以上所有步骤。
10.对于所有的杂交品种对，使用步骤8的平方差异，对均值方 法和区块克里金方法计算均方预测误差(MSPE)。
在这项有效性研究中，使用了Monsanto(DEKALB)研究系统的 12对玉米杂交品种的数据。研究结果表明，区块克里金方法的MSPE 比均值方法的MSPE低36％。注意，MSPE是预测误差的一种直接度 量，所以，较低的MSPE对应于较高的预测精确度。
下一步，我们讨论对本发明的估计组件进行的有效性研究。这项 分析是基于一项广泛的仿真研究，它以一种独特的方式，既使用真实 数据，又使用仿真数据，正如我们下面的介绍。首先回顾一下，本发 明的估计组件是基于广义最小平方(GLS)过程，它是对更传统的普通最 小平方(OLS)过程的一种改进。与GLS过程不同，OLS忽略数据中的 空间自相关，仅仅模拟大尺度的趋势。这项有效性研究的计划，是要 对比三种截然不同的方法：(a)GLS方法、(b)OLS方法和(c)传统的均 值方法，在品种对比的情况下，后者仅仅包括有关的位置匹配的、并 排种植的两个品种的差异均值的计算。这项有效性研究的主要目的之 一，是要了解先从均值方法到OLS方法，再从OLS方法到GLS方法， 是否有增加的利益。
下面提供了这项有效性研究过程各步骤的介绍。
1.(从1994年、1995年或1996年中)选择一年，再选择相对成 熟程度(RM)接近的一对玉米杂交品种。不考虑位置匹配，搜集这两个 杂交品种的全部研究试验数据。
2.整理数据，排除空间异常值，并用一个矩形框围绕整理后的数 据。我们将称之为数据矩形。这将表示可用于分析的样本数据。
3.在数据矩形中构建另一个内部矩形，在东西方向和南北方向都 覆盖1/3到2/3的数据范围。我们将命名这个矩形为估计矩形。这就是 我们将要进行估计的矩形。
4.使用数据矩形内的全部观测结果，对每个杂交品种分别拟合一 个三次多项式趋势曲面。拟合这个模型是通过GLS完成的，先拟合一 个球面协方差模型(参见背景技术一节)到每个杂交品种的残差协方 差，然后使用两个品种共同的估计熔核效应。注意，在这两个杂交品 种之间，并未假设相互的协方差。
5.在估计矩形中对两个品种及其差异(也就是两个拟合趋势曲面 的差异)，计算趋势曲面均值的估计结果。我们将称它们为原始估计 结果。
6.对于欧几里得距离0.01和0.05(分别对应于大约40和200英 里)，计算每个杂交品种的自相关和两个杂交品种之间的互相关。
7.使用对应趋势曲面的估计参数和协方差参数以及步骤4的熔核 效应，在数据矩形的0.01网格上，对每个杂交品种模拟数据。
8.在x方向和y方向都采用随机的振幅和频率产生两种正弦波。 对于第一种正弦波，随机产生之频率的均值和标准差分别为100和5， 而对于第二种正弦波它们分别是17和3。同样，对于第一种正弦波， 随机产生之振幅的均值和标准差分别为4和2，而对于第二种正弦波 它们分别是3和1。(注意，在模拟的数据加入这些正弦波，意味着 引入某种不易被三次模型拟合的系统性变化。现在的数据结果等于三 次曲面加正弦波加自相关后的随机误差。)
9.在估计矩形中，使用(通过步骤4中的GLS拟合出的)三次 曲面参数和(步骤8中的)正弦波参数，计算两个杂交品种均值的“真 实”值及其差异。这些真实的均值将用于在以下步骤中计算估计的偏 置和置信区间的覆盖范围。
10.独立地对每个杂交品种随机地去除40％的数据点，在网格点 上产生两个杂交品种的人工不匹配。
11.在估计矩形内部，仅在两种模拟值都有的网格点上，使用模 拟值计算杂交品种的平均值及其差异。这将称为均值估计结果。
12.对于估计矩形内外的全部样本值(也就是使用整个数据矩 形)，通过拟合三次多项式(两个杂交品种为不同的多项式)，以及 假设两个杂交品种的残差相互独立但是残差的方差不同，使用OLS估 计两个杂交品种(在估计矩形上)的均值及其差异。这就提供了OLS 估计结果。
13.使用估计矩形内外的全部样本值(也就是使用整个数据矩 形)，通过拟合三次多项式(两个杂交品种的不同)和一个带有熔核 效应和残差互相关参数的球面协方差结构(两个杂交品种的不同)， 使用GLS估计两个杂交品种(在估计矩形上)的均值及其差异。这就 提供了GLS估计结果。
14.对于使用(i)均值估计结果、(ii)OLS估计结果和(iii)GLS估计 结果的真实均值及其相应的标准差构建80％置信区间。标注置信区间 是否覆盖步骤9中获得的真实均值。对以上三个估计结果中的每一个， 还要计算对于真实均值的偏置。
15.重复步骤7到步骤14一百次。这是模拟的数目。
16.使用所有100次模拟的估计结果，对于(a)均值的均值估计结 果、(b)均值的OLS估计结果和(c)均值的GLS估计结果计算标准差。 通过采用步骤14中获得的偏置(包括所有100次模拟)的平均值，还 要计算三个估计结果中每一个的平均偏置。
17.对于均值、OLS和GLS方法中的每一种，计算(在100次模 拟中)置信区间覆盖真实均值次数的百分比。
18.对于覆盖宽广的RM范围的27对杂交品种和对于1994年、 1995年和1996年三年的所有数据，重复以上的所有步骤。
注：如上所述，注意到这种有效性研究以一个独特的方式结合了 实际数据和模拟结果，首先使用每对玉米杂交品种的实际数据获得趋 势和相关参数，然后使用这些参数产生模拟研究所用的数据。注意， 在这项有效性研究中不能仅仅使用实际数据的理由是，这个过程需要 “真实”均值来计算覆盖范围的百分比和偏置，如步骤14中的介绍。
我们现在归纳这项有效性研究的结果。
1.以标准差度量提高了精度，从均值方法到GLS方法为37％。 这种提高在从均值到OLS和从OLS再到GLS之间等分。
2.GLS方法的80％置信区间覆盖真实均值次数的百分比(77％) 非常接近名义值，而对于均值和OLS方法则比80％低很多(均值为 46％，OLS方法为56％)。
3.对于均值和OLS方法，这项研究揭示了真实标准差的总体欠 估计，这意味着对无差异的虚假设有太多的假拒绝，也就是品种之间 确实没有差异时，结论却是品种差异存在。
结论：
我们在这里介绍过的所有这些有效性研究，共同表明了本发明在 农作物品种性能的空间估计和空间预测方面的功效。
对于1999年10月15日提交的、标题为“SYSTEM FOR PLANT PERFORMANCE ANALYSIS”、序列号为No.60/159,802的临时专 利申请，本申请书要求其优先权。这里专门引用上述公开文件的全文 作为参考，没有放弃权利要求。
参考文献
这里专门引用下列参考文献和说明书中提及的所有参考文献作为 参考。
Besag，J.，and Kempton，R.(1986)，Statistical analysis of field experiments using neighboring plots，Biometrics，42，231-251.
Bhatti，A.U.，Mulla，D.J.，and Frazier，B.E.(1991)，″Estimation of soil properties and wheat yields on complex eroded hills using geostatistics and thematic mapper images，″Remote Sens.Environ.，37，181-191.
Bhatti，A.U.，Mulla，D.J.，Koehler，F.E.，and Gurmani，A.H.(1991)，″Identifying and removing spatial correlation from yield experiments，″Soil Sci.Soc.Am.J.， 55，1523-1528.
Bradley，J.P.，Knittle，K.H.，and Troyer，A.F.(1988)，″Statistical methods in seed com product selection，″J.Prod.Agric.，1，34-38.
Brownie，C.，Bowman，D.T.，and Burton，J.W.(1993)，″Estimating spatial variation in analysis of data from yield trials：a comparison of methods，″Agron.J.，85， 1244-1253.
Cook，R.D.(1977)，“Detection of influential observation in linear regression，” Technometrics，19，15-18.
Cressie，N.(1990)，″The origins of kriging，″Mathematical Geology，22，239-252.
Cressie，N.(1993)，Statistics for spatial data，Revised Edition，John Wiley andSons， New York，USA.
Cressie，N.，and Lahiri，S.N.(1996)，“Asymptotics for REML estimation of spatial covariance parameters，”Journal of Statistical Planning and Inference，50， 327-341.
Christensen，R.，Pearson，L.M.，and Johnson，W.(1992)，“Case-deletion diagnostics for mixed models，”Technometrics，34，38-45.
Cullis，B.，Gogel，B.，Verbyla，A.，and Thompson，R.(1998)，“Spatial analysis of multi-environment early generation variety trials，”Biometrics，54，1-18.
Gandin，L.S.(1959)，“The problem of optimal interpolation，”Trudy GGO 99，67-75.
Gandin，L.S.(1963)，Objective Analysis of Meteorological Fields， Gidrometeorologicheskoe Izdatel′stvo(GIMIZ)，Leningrad(translated by Israel Program for Scientific Translations，Jerusalem，1965).
Godambe，V.P.(editor)(1991)，Estimatingfunctions，Clarendon Press，Oxford.
Goldberger，A.S.(1962)，“Best linear unbiased prediction in the generalized linear regression model，”Journal of the American Statistical Association，57，369- 375.
Goovaerts，P.(1997)，Geostatistics for Natural Resources Evaluation，Oxford University Press，New York.
Gotway，C.A.，and Hartford，A.H.(1996)，“Geostatistical methods for incorporating auxiliary information in the prediction of spatial variables，”Journal of Agrricultural Biological and Environmental Statistics，1，1，17-39.
Gotway，C.A.，and Stroup，W.W.(1997)，“A generalized linear model approach to spatial data anlysis and prediction，”Journal of Agricultural Biological and Environmental Statistics，2，2，157-178.
Harville，D.A.(1977)，“Maximum likelihood approaches to variance component estimation and to related problems”Journal of the American Stattistical Association，72，320-340.
Heyde，C.C.(1994)，“A quasi-likelihood approach to the REML estimating equations，”Statistics and Probability Letters，21，381-384.
Hoaglin，D.C.，and Welsch，R.E.(1978)，“The hat matrix in regression and ANOVA，”American Statistician，32，17-22.
Kolmogorov，A.N.(1941)，“Interpolation and extrapolation of stationary random sequences，”Isvestiia Akademii Nauk SSSR，Seriia Matematicheskiia，5，3-14. (Translation，1926，Memo RM-3090-PR，Rand Corp.，Santa Monica，CA).
Littell，R.C.，Milliken，G.A.，Stroup，W.W.，and Wolfinger，R.D.(1996)，SAS System for Linear Models，SAS Institute Inc.，Cary，North Carolina，USA.
Mardia，K.V.，and R.J.Marshall.(1984)，“Maximum likelihood estimation of models for residual covariance in spatial regression，”Biometrika，71，135-146.
Matheron，G.(1962)，Traitéde Géostatistique Appliquée，Tome I.Mémoires du Bureau deRecherches Géologlogiques et Minières，14，Editions Technip，Paris.
Matheron，G.(1969)，“Le Krigeage Universel，”Cahiers du Centre de Morphologie Mathematique，1，Fontainebleau，France.
Montgomery，D.C.，and Peck，E.A.(1992)，Introduction to Linear Regression Analysis Second Edition，Wiley，New York.
Ovalles，F.A.，and Collins，M.E.(1988).“Evaluation of soil variability in northwest Florida using geostatistics，”Soil Sci.Soc.Am.J.，52，1702-1708.
Patterson，H.D.，and R.Thompson.(1971)，“Recovery of interblock information when block sizes are unequal，”Biometrika，58，545-554.
Patterson，H.D.，and R.Thompson.(1974)，“Maximum likelihood estimation of components of variance，”Pages 197-207in Proceedings of the 8th International Biometric Conference，Biometric Society，Washington，DC， USA.
Searle，S.R.(1971)，Linear Models，Wiley，New York.
Searle，S.R.，Casella，G.，and McCulloch，C.E.(1992)，Variance Components，Wiley， New York.
Sprague，G.F.，and Eberhart，S.A.(1976)，“Corn breeding，”Corn and Corn Improvement，J.A.Dudley and G.F.Sprague(eds)，Iowa State Univ.Press.
Stroup，W.W.，Baenziger，P.S.，and Mulitze，D.K.，(1994)“Removing spatial variation from wheat yield trials：a comparison of methods，”Crop Science，86， 62-66.
Ver Hoef，J.M.，and Barry，R.D.(1998)，“Modeling crossvariograms for cokriging and multivariable spatial prediction，”Journal of Statistical Planning and Inference，69，275-294.
Ver Hoef，J.M.，and Cressie，N.(1993)，“Multivariable spatial prediction，” Mathematical Geology，25，219-240.
Weiner，N.(1949)，Extrapolation，Interpolation，and Smoothing of Stationary Time Series，MIT Press，Cambridge，MA.
Wold，H.(1938)，A Study in the Analysis of Stationary Time Series，Almqvist and Wiksells，Uppsala.
Wolfinger，R.D.，Tobias，R.D.，and Sall，J.(1994)，“Computing Gaussian likelihoods and their derivatives for general linear mixed models，”SL4M Journal on Scientifc Computing，15，1294-1310.
Wu，T.，Mather，D.E.，and Dutilleul，P.(1998)，“Application of geostatistical and neighbor analyses to data from plant breeding trials，”Crop Science，38，1545- 1553.
Yost，R.S.，Uehara，G.，and Fox，R.L.(1982a)，“Geostatistical analysis of soil chemical properties of large land areas.I.Semi-variograms，”Soil Sci.Soc.Am. J.，46，1028-1032.
Yost，R.S.，Uehara，G.，and Fox，R.L.(1982b)，“Geostatistical analysis of soil chemical properties of large land areas.II.Kriging，”SoilSci.Soc.Am.J.，46， 1033-1037.
Zimmerman，D.L.(1989)，“Computationally efficient restricted maximum likelihood estimation of generalized covariance functions，”Mathematical Geology，21， 655-672.
Zimmerman，D.L.，and Harville，D.A.(1991)，“A random field approach to the analysis of field-plot experiments and other spatial experiments，”Biometrics， 47，1，223-239.