总结

我们提出了一类半参数函数回归模型来描述向量值协变量对响应曲线样本的影响。每个观察到的曲线都被视为一个随机过程的实现,由一个总体平均函数和随机分量组成。有限维协变量通过包括未知光滑连接和方差函数的单指标模型影响本征函数展开的随机分量。单指标模型的参数分量是通过拟核估计方程估计的,其中链接函数和方差函数是非参数估计的。我们得到了几个基本的渐近结果。通过对1000只地中海雌果蝇(medflies)产卵曲线数据集的分析,说明了所提出的函数回归模型。

1.简介

由于这类数据的收集越来越普遍,因此需要对由曲线或类似无限维对象组成的数据进行分析的方法。目前正在收集内分泌学(Brumback和Rice,1998)和生长研究(Gasser.,1984; Gasser和Kneip,1995; Kneip和Gasser,1992). 在Ramsay和Silverman中可以找到被视为数据的函数所使用的不同方法和方法的概述(1997,2002).

我们在这里采用的观点是,每个观察到的曲线都是一个随机过程的实现,反映了曲线样本中包含的各个曲线的随机性质。对于曲线数据,该方法由卡斯特罗首创. (1986)并在莱斯和西尔弗曼的开创性工作中得到进一步发展(1991)他强调了平滑的重要性。对于许多应用,将协变量或相关变量之间的关系建模为预测因子和为每个受试者获得的响应曲线是特别有趣的。之前在卡普拉和米勒(Müller)研究过以时间偏移的形式纳入协变量效应(1997),其中一个协变量Z轴假设通过乘法方案影响时间刻度,t吨=ψ(Z轴)t吨,对于未知但平滑的函数ψ这可以从曲线样本中估计。将协变量纳入函数数据模型的其他最新方法包括方差函数分析方法的发展(Brumback和Rice,1998; 斯坦尼斯瓦利斯和李,1998)变系数回归模型(胡佛等。,1998). Shi在最近的工作中提出了一个相关的随机效应模型,其中将观测曲线假设为具有随机系数的回归样条曲线等。(1996)赖斯和吴(2001).

我们的出发点是随机过程的样本X(X)(t吨),= 1,…,n个,t吨T型,其中T型是一个间隔。实际上,流程X(X)通常记录在时间点网格上t吨ij公司, 1 ⩽j个,其中我们假设它们是等距的,或者如果不是等距的话,则在由设计密度生成的规则网格上。对于每个过程X(X),我们观察到协变量或相关变量的向量Z轴,Z轴∈ ℝ第页,第页⩾ 1. 观察结果{Z轴,X(X)假设(·)}是独立的同分布的。流程X(X)(·)假设受预测因素影响的响应曲线的作用Z轴假设总平均函数μ(·)和特征函数或主成分函数ρk(·),k=1,2,…,通过过程的自方差结构定义X(X)(请参见附录A详细信息)不依赖于协变量Z轴,而主成分(也称为主成分得分)和过程的条件平均值X(X)do.高维协变量向量的情况Z轴通过将协变量向量减少为单个索引来处理。然后将该单一指数嵌入到半参数拟似然回归(SPQR)模型中。我们使用具有未知链接和方差函数估计的拟似然法或QLUE方法(参见Chiou和Müller(1998))它提供了一个含有非参数估计分量的估计方程、渐近分布理论以及一个拟合参数和非参数分量的算法。我们通过生育率数据证明了所提出的半参数方法的有效性,其中曲线样本由1000只地中海雌果蝇(地中海果蝇)的每日产卵量组成;有关插图,请参见图1.

20个随机选择的平滑产卵曲线,用于总产卵量的四分之一中的每一个:(a)第一个四分之一;(b) 第二四分位数;(c) 第三四分位数;(d) 第四个四分位数
图1

为卵子总数的四个四分位数中的每一个随机选择20条平滑的产卵曲线:(a)第一个四分位数;(b) 第二四分位数;(c) 第三四分位数;(d) 第四个四分位数

论文组织如下:函数平滑随机效应模型曲线数据的协变量效应建模见第2节中讨论了模型分量的估计、平滑过程和渐近结果第3节将此方法应用于产蛋数据是第4节其中还包括对各种实际问题的讨论,例如特征函数数目和平滑参数的选择。其他技术材料,如特征函数方面的过程表示、SPQR估算方程的证明和细节,见附录A-C。

2.函数平滑随机效应模型

数据向量的主成分分析是一种常用的多元技术。曲线样本的分析可以以类似的方式进行。卡斯特罗. (1986)定义了最佳的概念K(K)-随机过程的多维线性模型,并表明该模型对应于第一个K(K)基于过程协方差核的特征函数的Karhunen–Loève展开的项。Rice和Silverman(1991)讨论了在曲线数据分析中实现这一概念的平滑技术的应用,并提出了一种选择平滑参数的留一曲线交叉验证技术。

预测变量Z轴∈ ℝ第页随过程观察X(X)(t吨),t吨T型,其中T型是一个间隔。我们假设Z轴影响随机分量A类k在以下一类函数平滑随机效应模型(模型1)中。

  • 给定平滑函数μ(t吨),t吨T型,随机变量A类k,k=1,2,…和本征函数ρk:T型→ℝ,
    X(X)(t吨)=μ(t吨)+k=1A类kρk(t吨),
    并假设观测到的随机曲线以协变量为条件,
    μZ轴(t吨)=E类{X(X)(t吨)Z轴}=μ(t吨)+k=1E类A类kZ轴ρk(t吨).

进行了以下假设。

  • (a)

    有一个整数K(K)<∞,这样特征值λk=E类A类k2(请参见方程(17)英寸附录A)满足λk<∞k= 1,…,K(K)、和λk=0,对于k>K(K)意味着只有第一个K(K)上述扩展条款。

  • (b)

    本征函数ρk(·)在中是正交的L(左)2([0,T型])和是两倍连续可微的。

  • (c)
    对于所有人k,使用Kronecker符号δ肯尼亚=1如果k=δ肯尼亚=0,如果k,
    E类A类k=0冠状病毒A类k,A类=λkδk.
    从假设(c)中,我们发现条件协方差(和设置=t吨,以及条件方差)
    覆盖(cov){X(X)(),X(X)(t吨)Z轴}=k,覆盖(cov)A类k,A类Z轴ρk()ρ(t吨).
    (1)
  • (d)
    随机效应的每个条件链接函数都是由数据形成的单个指数的平滑函数;即有一个参数向量βk∈ ℝ第页, ║βk║=1和一个二次连续可微函数αk:→这样
    αkβkz(z)=E类A类kZ轴=z(z),
    对于k= 1,2,….

这些单索引模型用作降维设备。

所提出的函数平滑随机效应模型是随机过程经典Karhunen–Loève表示的自然扩展(参见方程式(18)英寸附录A). 它包含了协变量效应的影响Z轴通过主成分得分的条件分布,即模型中的随机效应,通过平滑回归关系(d)。

链接功能αk在假设(d)中,它们与本征函数一起具有内在的意义ρk,确定流程依赖性的性质X(X)关于预测因子Z轴对于给定的本征函数,通过检查本征函数具有较大主成分的过程,可以找到对潜在主题问题的解释;这些“主成分得分”又由特征函数的形状决定。基于对k第个特征函数ρk,链接功能αk然后揭示了假设(d)和模型1作为协变量表示中该特征函数相对重要性的变化Z轴变化。

我们以三句话结束本节。首先,模型的条件协方差结构,如方程式(1)具有很强的适应性,同时保持非参数分量的维数较低,以避免与高维非参数平滑相关的“维数灾难”。

其次,如果域是T型= (t吨0)那么,一个单身汉X(X)只是一个真实的随机变量,并且只有一个本征函数:ρ1(t吨0) = 1. 在这种情况下,模型1崩溃为具有未知链接函数的半参数回归模型(·)和参数向量θ∈ ℝ第页这样的话E类{X(X)(t吨0)|Z轴=z(z)} =(θz(z)),这取决于(Z轴,X(X)(t吨0)). 这将我们的方法与具有未知链接和方差函数的拟似然回归模型联系起来(参见附录C). 从这个意义上讲,模型1提供了SPQR的自然扩展,适用于响应是函数的情况。

第三,考虑一个带有显式误差过程的修正模型,如下所示X(X)˜(t吨)=X(X)(t吨)+ε(t吨),其中错误处理ε(t吨)独立于X(X)Z轴并满足E类{ε(t吨)}=0和cov{ε(),ε(t吨)} =σ(,t吨),其中σ(·,·)是平方可积的。这意味着流程X(X)˜具有特征函数ρ˜k以及存在不相关的零均值随机变量A类˜k这样的话

A类˜kρ˜k(t吨)=A类kρk(t吨)+ε(t吨).

因此,基于随机效应,模型1仍然适用A类˜k和本征函数ρ˜k以及相应的修改解释。

3.模型组件的估计

3.1. 前期工作

为了进行估计和预测,我们假设K(K)模型1中具有非零特征值的特征函数是有限的。如果实际情况下K(K)组件,目标是第一个组件跨越的空间上的投影K(K)本征函数。我们还假设从现在开始T型流程的T型= [0,τ]以及离散化测量X(X)(t吨ij公司)进程数X(X)(·)可用,= 1,…,n个j个= 1,…,。我们提出了等间距设计的论点t吨ij公司=t吨j个,j个= 1,…,、和t吨j个由平滑的设计密度生成(例如,参见Müller(1984)),但是,如果测量次数不同,则可以通过添加预平滑步骤来覆盖不等距的情况以及测量次数因受试者而异的情况的扩展对于第个主题满足0<c(c)1/c(c)2<∞,= 1,…,n个,对于常量c(c)1,c(c)2>0.

由于平滑在所提出的方法中起着核心作用,我们简要讨论了一些相关问题。给定任意散点图数据(W公司,Y(Y))= 1,…,,具有基本回归函数(w个) =E类(Y(Y)|W公司=w个),我们定义平滑估计^()通过

^(w个)=S公司w个,b条,W公司,Y(Y)=1,,.

在这里,是散点图中的数据数量b条是平滑器的平滑参数S公司,在参数处求值w个.

对于我们的分析,只要满足一些基本的一致性属性并且平滑器在数据中是线性的,那么应用哪种特定的平滑方法并不重要。核估计、平滑样条或局部加权最小二乘的局部多项式拟合是可能的选择。我们在这里选择局部加权最小二乘平滑器,表示为S公司L(左),根据数据拟合本地线。正式定义是

S公司L(左)w个,b条,W公司,Y(Y)=1,,=参数最小值0最小值1=1K(K)w个W公司b条Y(Y)0+1w个W公司2.

这意味着参数中更平滑的值w个作为回归线在w个,其中该行适用于掉入窗口的数据[w个b条,w个+b条]通过加权最小二乘法。在这里,K(K)(·)⩾0是一个非负核函数,常见的选择是K(K)(x个)=(1−x个2)1(|x个|1)K(K)(x个)=经验(−x个2/2).

3.2. 估计总体均值和特征函数

总体平均函数的粗略估计μ可以通过平滑逐点平均值来获得μ˜()对于等间距设计,这些逐点平均值为

μ˜t吨j个=1n个=1n个X(X)t吨j个,
(2)

鉴于n个采样曲线X(X)1,X(X)2,…,X(X)n个。要在中获得平滑函数t吨,我们可以使用较小的带宽添加平滑步骤b条,

μ^(t吨)=S公司t吨,b条μ,t吨j个,μ˜t吨j个j个=1,,.
(3)

对于非等间距设计,可以首先执行预平滑步骤,目的是在同一时间点对所有曲线进行采样。

设∑为样本×随机过程的方差-协方差矩阵X(X),包含元素

v(v)第页=1n个=1n个X(X)t吨第页μ^t吨第页X(X)t吨μ^t吨,1第页,.
(4)

表示方式ρ˜k=ρ˜kt吨1,ρ˜kt吨2,,ρ˜kt吨对应于k矩阵∑的第个最大特征值。遵循卡普拉和米勒使用的方法(1997)这是赖斯和西尔弗曼提议的变体(1991),一个平稳的估计ρ^k()本征函数的ρk(·)只需对矢量进行平滑即可获得ρ˜k,

ρ^k(t吨)=S公司t吨,b条ρk,t吨j个,ρ˜kt吨j个j个=1,,,k=1,,K(K).
(5)

3.3. 估计随机影响

链路功能估计αk和参数向量βk基于估计的总体平均函数μ^和本征函数ρ^k.离散近似值A类k=〈ρk,X(X)μ〉(见表达式(19)英寸附录A)然后激发更新的随机效应估计

A类^k=Δj个=1X(X)t吨j个μ^t吨j个ρ^kt吨j个,k=1,,K(K),
(6)

式中Δ=t吨ij公司t吨(j个−1).更新估算α^kβ^kZ轴,用于Z轴∈ ℝ第页第页>1,我们采用

β^k,α^k()=昆士兰州Z轴,A类^k=1,,n个,
(7)

以Chiou和Müller的QLUE方法的形式使用SPQR(1998)总结如下附录C.对于协变量Z轴是一维的,我们可以得到α^k()通过简单的非参数回归步骤。在这种情况下,没有参数向量,我们可以回归A类^k非参数地Z轴,= 1,…,n个,以获得平滑函数估计

α^k(v(v))=S公司v(v),b条αk,Z轴,A类^k=1,,n个.
(8)

3.4. 预测和估计算法

对于预测,给定一个协变量向量Z轴如中所述第2节,我们替换未知函数μ, {αk}k= 1,…,K(K)和{ρk}k= 1,…,K(K)在具有估计值的模型中μ^,α^kρ^k分别是。预测μZ轴(t吨) =E类{X(X)(t吨)|Z轴}就是那个时候

μ^Z轴(t吨)=μ^(t吨)+k=1K(K)α^kβ^Z轴ρ^k(t吨).
(9)

模型组件的估计过程总结如下。

  • (a)

    获取估算值μ^(t吨)总平均函数的μ(t吨)由方程式(2)(3).

  • (b)

    获取估算值ρ^k(t吨)本征函数的ρk(t吨)由方程式(5).

  • (c)

    获取链接函数拟合α^k()和参数向量估计β^k通过方程式(7)QLUE方法。

平滑和QLUE步骤需要选择平滑参数。QLUE算法中使用的选择是基于拟似然的特殊特征;见Chiou和Müller(1998). 函数估计的平滑参数选择ρ^kμ^0需要其他方法。我们将在第4.2节和第4.4节中讨论与留一曲线技术相结合的实际选择。

3.5. 一些基本渐近结果

遵循佩祖利和西尔弗曼的论点(1993),可以获得特征函数估计的一致性属性ρ^kk=1,,K(K),在温和的规则条件下。对于函数ϕL(左)2(d)ν)在域上T型,套

φ=T型φ2(t吨)d日ν(t吨)1/2.

这里我们考虑模型1的简化版本,其中Z轴∈ ℝ1因此βk=1用于k= 1,…,K(K).

为了更平滑S公司在平滑步骤中,我们需要以下两个最小属性,几乎所有常用的平滑方法都满足这两个属性。

  • (a)
    给定独立同分布随机对(W公司,Y(Y))= 1,…,,来自具有回归函数的分布(w个) =E类(Y(Y)|W公司=w个),如果回归函数(·)是两倍连续可微的,概率分布函数(f)W公司W公司s在某一点上是连续的w个并满足(f)W公司(w个)>0,然后
    S公司w个,b条,W公司,Y(Y)=1,,(w个)=O(运行)第页τn个,
    对于序列τn个0,这取决于特定的平滑器和平滑参数选择。
  • (b)
    数据中的平滑度是线性的Y(Y),
    S公司w个,b条,W公司,Y(Y)=1,,==1G公司(w个)Y(Y)
    对于(可能是随机的)权重函数G公司,独立于Y(Y),并满足
    =1G公司2(w个)=11G公司(w个)0=O(运行)第页(1).

满足这些属性的平滑器示例包括局部线性拟合或核平滑器等。

定理1。如果平滑器S公司用于估计非参数模型组件的参数满足属性(a)和(b),并且协变量Z轴在任何给定点具有连续且为正的密度z(z),然后

μ^E类X(X)=O(运行)第页n个1/2,
(10)

如果另外ρ^kρk=O(运行)第页σn个,k=1,,K(K)对于序列σn个0,然后

|α^(z(z))α(z(z))|=O(运行)第页σn个+τn个+n个1/2+1.
(11)

这个结果的证据在附录B一致性要求(毫不奇怪)n个观察到的样本曲线以及测量次数每条曲线满足n个→∞ →∞.

我们注意到,对于方程式(11):

啜饮z(z)|α^(z(z))α(z(z))|=O(运行)第页νn个σn个+τn个+n个1/2+1,

如果我们要求的不是财产(b)

啜饮z(z)=1n个G公司2(z(z))=1n个1G公司(z(z))0=O(运行)第页νn个.

通过大多数平滑方法获得的典型值是νn个={日志(n个)/编号}−1/2,其中b条是等效带宽。因此,一致性要求另外νn个(σn个+τn个+n个−1/2+−1)→0.在某些假设下,Pezzulli和Silverman(1993)已经证明了ρ^kρk=第页(1)从而得到一致的结果|α^(z(z))α(z(z))|=第页(1).

利用Chiou和Müller中描述的QLUE的渐近一致性,这些结果可以推广到多元预测的情况(1998). 对于未知链路函数的QLUE估计的以下结果αk和参数向量βk为了简单起见,我们假设μρk是已知的,因此{A类伊克}在SPQR步骤中扮演通常的“响应”变量的角色。

定理2。假设总平均函数μ和本征函数ρk模型1中给出了。然后,在规则条件下,以下事件对每个k,k= 1,…,K(K),概率为1−δ,对于给定的任意小δ>0:存在QLUE估计β^k参数向量的βk,令人满意β^k=1、和α^k链接功能的αk,例如n个→∞,

啜饮z(z)|α^(z(z))α(z(z))|=第页(1),
(12)
n个1/2(f)β^k(f)βkD类N个第页0,(D类(f))βkΣ1(D类(f))βkT型,
(13)

哪里

Σ=n个1n个=1n个D类T型V(V)1D类

具有D类=αβkT型x个/βkT型,(f)(βk) = ((f)1(βk),…,(f)第页−1(βk))T型具有(f)j个(βk) =β千焦/║βk

(D类(f))βk=(f)βkβkj个1第页1,1j个第页,

定义为(第页−1)×第页矩阵。

我们注意到,由于中的假设(d)中的可识别性约束,回归参数估计的自由度以及相应的渐近协方差估计的秩减少了1第2节定理证明的细节,包括渐近协方差估计的调整,与Chiou和Müller中的类似(1998,2002)因此省略了。

4.medflies产卵曲线的应用

4.1. 背景

最近,人们对老龄化和定量生物人口学的研究越来越感兴趣,关注的焦点是生殖和老龄化模式之间的关系。进化是由生殖成功驱动的,衰老的进化和生殖之间的联系非常有趣,尤其是因为它被证明是难以捉摸的。解开这一联系可能有助于我们理解衰老过程。繁殖模式通常是从涉及大型实验队列的实验中推断出来的。大型群蜻蜓繁殖行为的研究(地中海实蝇)得出了一些有趣的发现。例如,正如Müller所报道的那样,营养和性别影响之间的相互作用导致蛋白质缺乏情况下女性和男性预期寿命差异的逆转. (1997)在Carey、Liedo、Müller、Wang和Vaupel中发现了两种不同的衰老和繁殖模式以及饮食(1998). 据推测,生殖活动对寿命有负面影响,Partridge和Harvey发明了“生殖成本”一词(1985). 苍蝇是研究繁殖和寿命的理想工具,因为可以通过记录每天的产卵量和存活率来饲养大量苍蝇。这里,我们将重点放在繁殖曲线中的模式之间的关系上,繁殖曲线被视为相关数据,苍蝇产卵总数和寿命是预测因子。产卵总数是衡量生殖成功的一个指标,研究生殖模式与总体生殖成功和寿命的关系是很有意义的。

为我们的分析提供数据的实验于1992年至1995年在墨西哥恰帕斯Metapa的medfly大规模饲养和消毒设施(Moscamed)进行,包括n个=1000只记录了每日产卵量的雌性蜉蝣。有关实验生物学特性的更多细节,请参阅凯里、利多、米勒、王和邱(1998).

X(X)(·)表示第次飞行和Z轴食肉动物的寿命和卵子总数。随机子样本产卵曲线显示在图1,显示总卵子数每四分之一的20条曲线。从1000只苍蝇的初始样本中,那些没有产卵的苍蝇被丢弃,因为这些曲线的值为0,留下的产卵数据为n个=分析中936只苍蝇。我们将分析限制在产蛋的前50天,因为我们发现在较高年龄段产蛋的变化很大,这将完全支配特征函数,并且作为凯里、利多、米勒、王和邱的早期分析(1998)结果表明,卵的总数与寿命的关系在51天时出现了显著的变化点。因此,产卵曲线X(X)(t吨),=1,…,936被视为随机过程的实现T型= [0,τ]与τ=50天。

由于许多苍蝇在50天之前死亡,因此有必要考虑是否会出现信息缺失,因为产卵的平均功能会受到信息缺失模式的影响。建立了产卵过程的近似参数模型,并将其用于米勒. (2001)定义死亡后的假定产卵随机过程,特别是剩余产卵潜能随机函数。虽然在寿命和死亡时的剩余产卵潜能之间建立了有趣的关联,但发现寿命和死亡后的估算剩余产卵过程之间的关联至多非常微弱,在本例中,为假设退出是非信息性的提供了一些理由。

4.2. 本征函数

前四个估计特征函数ρ^k (5)如所示图2标题中说明了每个主成分的方差比例。相应的特征值范围从0到6000以上,这表明变化很大,与图1第一个特征函数大致与产卵曲线的总体平均形状相似,在前10天急剧上升,随后相对平坦的向下倾斜。

产蛋数据的前四个特征函数ρ^kk=1,…,4(5):第一个特征函数解释了数据总方差的34.38%,第二个特征函数额外解释了16.41%,第三个特征函数另外解释了8.82%,第四个特征值额外解释了4.93%(带宽是通过曲线交叉验证单独选择的)
图2

前四个本征函数ρ^kk=1,,4 (5)产蛋数据中:第一个特征函数解释了34.38%,第二个特征函数额外解释了16.41%,第三个特征函数另外解释了8.82%,第四个特征函数又解释了4.93%的总方差(带宽是通过曲线交叉验证单独选择的)

如果我们真的这么想ρ1(t吨) =μ(t吨),并且只有一个主要成分很重要,我们得到了乘法效应模型

E类{X(X)(t吨)Z轴=z(z)}=μ(t吨)1+E类A类1Z轴=z(z)

E类{X(X)(t吨)Z轴=z(z)}=μ(t吨)1+E类A类1Z轴=z(z).

有条件地打开Z轴=z(z)因此,我们观察到准γ型行为。无条件地,

E类{X(X)(t吨)}=μ(t吨),无功功率,无功功率{X(X)(t吨)}=μ2(t吨)λ1,

而“第一变化模式”是“在平均函数的方向上”μ(·). 因此,无条件地,我们在t吨,√λ1是变异系数。

我们发现,第二个特征函数将这种形状区分为大约5天的早期峰值和随后更快的下降,为以后的进一步峰值留出了空间。第三个本征函数描述了后期的宽峰值,而高阶本征函数将其分解为一系列日益复杂的振荡。

均值函数μ^ (3)不取决于协变量,如图3平滑参数是通过最小化通过留一条曲线技术获得的平方预测误差来确定的,

体育课==1n个j个=1X(X)^()t吨j个X(X)t吨j个2/n个,
(14)

哪里X(X)^()是通过省略估计中的样本曲线。因为对于没有协变量的经典Karhunen–Loève模型X(X)^()(t吨)=μ^()(t吨)=μ^()(t吨,b条),我们可以计算最佳预测平滑参数b条^μ=参数最小值b条{体育课(b条)}; 这个产生了b条^μ=2.5天。

估计时间μ^(t)(3)的总平均产卵函数,留一条曲线交叉验证2.5天的带宽
图3

估计总平均产卵时间函数μ^(t吨) (3),留出一条曲线,交叉验证带宽为2.5天

4.3. 作为单个协变量的卵子总数

我们首先考虑单个协变量的情况。协变量Z轴是鸡蛋的总数,这是随机分量的预测值A类k.除了估计特征函数外ρk和平均函数μ我们需要得到估计α^k (8)对于随机效应链接函数αk执行这些步骤可以得到估计的功能α^k,如所示图4对于k= 1–4. 我们发现了α^1是单调递增的,几乎是线性的,而α^2在开始时略有增加,然后随着负斜率的增加而单调减少。功能α^k更高阶的(k=3,4及更高)表现平缓,表明协变量仅影响前两个特征函数的得分。实际上,我们的最终模型不会使用所有四个随机成分。

平滑的随机效应链接函数α^k(z)k=1,…,4(8),依赖于总卵子数,基于936只medflies样本(交叉验证为所有四条曲线选择的带宽为500)
图4

平滑随机效应链接功能α^k(z(z))k=1,,4 (8)基于936只水飞蓟的样本(交叉验证选择的四条曲线的带宽均为500)

给定的功能估计μ^()α^k()k=1,我们得到拟合曲面

μ^Z轴(t吨)=μ^(t吨)+k=1K(K)α^k(Z轴)ρ^k(t吨),
(15)

根据方程式(9).配合面μ^Z轴(t吨) (15),如所示图5(a)展示了一个有趣的倾斜山脊模式,横截面图5(b)研究表明,尽管随着卵子总数的增加,峰值略微向右漂移,但随着卵子总量的增加,图表大部分都向上移动。我们注意到本征函数的数目被选为K(K)=2用于图5。此选择基于对功能的目视检查α^k在里面图4以及拟合模型的预测质量。

(a) 以卵子总数为协变量的拟合曲面μ^Z(t)(15),以及(b)通过拟合曲面的横截面μ^Z
图5

(a) 配合面μ^Z轴(t吨) (15)以卵子总数为协变量,(b)通过拟合曲面的横截面μ^Z轴(t吨)对于固定在400(--)、800(●●●●)、1200(––)和1600(--)的卵子总数

留一曲线预测误差(14)是模型预测质量的有用量化。保留一条曲线的预测因素是

μ^Z轴()(t吨)=μ^()(t吨)+k=1K(K)α^k()Z轴ρ^k()(t吨).
(16)

在这里,μ^(),ρ^k()α^k()都是以与之前相同的方式获得的,但忽略了观察到的过程X(X)自身。Rice和Silverman引入了交叉验证平方和的函数形式(1991). 拟合特征函数数量依赖性的预测误差如所示表1,支持选择K(K)= 2. 对于没有协变量效应的模型,即方程式(18)英寸附录A结果为462.54。当包含协变量时,预测误差的大幅减少证明了在模型1中包含协变量效应的有效性。

表1

预测误差PE(14)拟合模型的(16)具有K(K)本征函数

K(K)1245678
体育课333.71313.22315.22316.41315.64315.62315.62315.49
K(K)1245678
体育课333.71313.22315.22316.41315.64315.62315.62315.49
表1

预测误差PE(14)拟合模型的(16)具有K(K)本征函数

K(K)1245678
体育课333.71313.22315.22316.41315.64315.62315.62315.49
K(K)1245678
体育课333.71313.22315.22316.41315.64315.62315.62315.49

4.4. 多协变量的情况

对于多协变量的情况,我们考虑产卵总数和寿命作为随机分量的预测因子A类k中描述的估算程序第3节为标准化协变量实施。分析中使用的随机成分的尺寸再次选择为K(K)= 2. 前四个估计的随机效应链接函数α^k在中图6.

随机效应α^kβk′Zk=1,…,4(7)的估计链接函数,由SPQR的QLUE实现获得,依赖于线性预测因子描述
图6

随机效应的估计链接函数α^kβkZ轴k=1,,4 (7),根据线性预测器描述,通过SPQR的QLUE实现获得

估计系数β^1=β^11,β^12β^2=β^21,β^22函数中的线性预测器α1α2显示在中表2这些估计值是通过QLUE实施的SPQR方法获得的。相关的标准误差对渐近推断很有用。对于线性预测器β1卵子总数相对来说比寿命更重要,这两个预测因子都显著且符号相反,在协变量中形成对比。对于线性预测器β2,两个预测因子的影响大致相同,形成协变量的平均值,且两个预测因素都显著。

表2

随机成分的QLUEαk(βkV(V))对应于图6

功能系数QLUE(标准误差)带宽贴合度好
α1β11(总蛋数)0.8842 (0.0031)83.78(方差)D类= 940.36
β12(终身)−0.4671 (0.0059)1.16(链接)P(P)= 916.20
α2β21(总蛋数)0.6730 (0.0520)28.22(方差)D类= 1250.75
β22(终身)0.7397 (0.0471)1.91(链接)P(P)= 936.60
功能系数QLUE(标准误差)带宽适合度
α1β11(总蛋数)0.8842 (0.0031)83.78(方差)D类= 940.36
β12(终身)−0.4671 (0.0059)1.16(链接)P(P)= 916.20
α2β21(总蛋数)0.6730 (0.0520)28.22(方差)D类= 1250.75
β22(终身)0.7397 (0.0471)1.91(链接)P(P)= 936.60
表2

随机成分的QLUEαk(βkV(V))对应于图6

功能系数QLUE(标准误差)带宽适合度
α1β11(总蛋数)0.8842 (0.0031)83.78(方差)D类= 940.36
β12(终身)−0.4671 (0.0059)1.16(链接)P(P)= 916.20
α2β21(总鸡蛋数)0.6730 (0.0520)28.22(方差)D类= 1250.75
β22(终身)0.7397 (0.0471)1.91(链接)P(P)= 936.60
功能系数QLUE(标准误差)带宽适合度
α1β11(总蛋数)0.8842 (0.0031)83.78(方差)D类= 940.36
β12(终身)−0.4671 (0.0059)1.16(链接)P(P)= 916.20
α2β21(总蛋数)0.6730 (0.0520)28.22(方差)D类= 1250.75
β22(终身)0.7397 (0.0471)1.91(链接)P(P)= 936.60

单个预测因子对响应的影响是下列系数的结果表2和中的本征函数图2对于当前的特定应用,我们在此提出了一个用于估计链路函数的必要平滑参数选择准则,该准则易于实现,并且具有良好的实用性。它基于伪似然(见Davidian和Carroll(1988)),

损益;μ(β),σ2(μ)=12=1n个日志2πσ2μ12=1n个μ2σ2μ,

哪里,μσ2分别表示响应向量、均值向量和方差向量。在这里μσ2被他们的估计所取代μ^σ^2在伪似然中,通过最大化伪似然PL来选择最佳带宽(;μ,σ2).

数据分析的结果带宽如第四列所示表2此外,在表2我们提供了SPQR的优良率统计数据。这里,我们提供了“非参数”准偏差D类和皮尔逊统计P(P),它们是渐近的χ2分布的自由度近似等于观测次数减去线性预测器中待估计参数的数量。中的值表2不要表示不适合,因为D类P(P)(比较Chiou和Müller(1998),了解更多详细信息)。

5.讨论和总结

我们开发了一种方法,将曲线数据的主成分分析扩展到存在协变量的情况。这种方法扩展了Rice和Silverman以前的工作(1991)并产生了一个灵活的函数回归模型,即函数平滑随机效应模型。该模型在均值回归和异方差结构方面具有很大的灵活性。通过模型的单索引功能,高维度保持了灵活性。通过对1000只雌性瓢虫产卵曲线的计算,表明该方法能有效地减少预测误差。提出的方法在概念上和计算上都很简单。SPQR模型的非参数组件的结构只需要一维平滑步骤。

我们的主要结果和技术是针对等间距设计得出的,其中曲线数据的可用测量是在规则网格上进行的。非等间距数据在实践中显然很重要。如果一个设计是“适度”非等距的,即没有持续的间隙第3.1节仍然可以应用;对于测量和稀疏设计中持续存在的差距,需要对提出的方法进行扩展。

通过各种方法,可以扩展到与实际相关的情况,即预测器是一个时变函数而不是一个时不变向量。其中一种方法是离散化协变量函数Z轴(t吨)变成一个(可能是高维的)矢量Z轴(1),…,Z轴(M(M))通过指定等距网格M(M)1,…,M(M)然后用这些向量作为协变量,对上述函数平滑随机效应模型进行拟合。另一种方法是将随机协变量函数投影到它们的第一个K(K)某些合适整数的特征函数K(K)'并将所得主成分得分向量用作K(K)函数平滑随机效应模型中的′维协变量向量。

在所提出的方法中,对于仅作用于随机效应的一维协变量,可以在一步内轻松一致地估计所有模型分量,并且我们提供了一个渐近一致性结果。通过将这些思想与多元和高维回归的单指标拟似然方法相结合,使用QLUE技术,实现了对多元预测因子的扩展。这允许我们同时包含未知链接和方差函数。我们的实现还提供了平滑参数的自动选择,因此只需要识别协变量,不需要对随机效应分布的性质或链接函数的形状作出任何规定。

致谢

我们感谢J.Carey教授为我们提供了medfly繁殖力数据。我们也感谢两位审稿人,一位副主编和一位联合主编所作的最有帮助的评论。本研究部分得到了国家卫生研究院拨款BS-090-PP-07、国家科学委员会拨款89-2118-M-194-001、国家科学基金会拨款DMS-99-71602和DMS-02-04869以及国家卫生研究所拨款99-SC-NIH-1028的支持。

工具书类

布伦巴克
,
B。
大米
,
J。
(
1998
)
曲线嵌套和交叉样本分析的平滑样条模型(附讨论)
.
美国统计学杂志。助理。
,
93
,
961
994
.

卡普拉
,
W.B.公司。
米勒
,
H.G.公司。
(
1997
)
响应曲线的加速时间模型
.
美国统计学杂志。助理。
,
92
,
72
83
.

凯里
,
J。
,
列多
,
第页。
,
米勒
,
H.G.公司。
,
,
J·L·。
,
J·M·。
(
1998
)
地中海果蝇雌性大队列中繁殖力年龄模式与死亡率、寿命和终生繁殖的关系
.
J.Geront。A类
,
53
,
B245型
B251型
.

凯里
,
J。
,
列多
,
第页。
,
米勒
,
H.G.公司。
,
,
J·L·。
沃佩尔
,
J。
(
1998
)
地中海果蝇雌性的双重衰老模式
.
科学类
,
281
,
996
998
.

卡斯特罗
,
体育。
,
劳顿
,
W.H.公司。
西尔维斯特
,
欧洲航空公司。
(
1986
)
具有连续样本曲线的过程的主要变化模式
.
技术计量学
,
28
,
329
337
.

,
J·M·。
米勒
,
H.G.公司。
(
1998
)
链函数和方差函数未知的拟似然回归
.
美国统计学杂志。助理。
,
93
,
1376
1387
.

,
J·M·。
米勒
,
H.G.公司。
(
1999
)
非参数拟似然
.
安。统计师。
,
27
,
36
64
.

,
J·M·。
米勒
,
H.G.公司。
(
2002
)
具有光滑连接和方差函数的多指标拟似然回归
。待发布。

达维迪安语
,
M。
卡罗尔
,
R·J。
(
1988
)
关于扩展拟似然的一个注记
.
J.R.统计。Soc.B公司
,
50
,
74
82
.

气体发生器
,
T。
膝盖
,
答:。
(
1995
)
搜索曲线样本中的结构
.
美国统计学杂志。助理。
,
90
,
1179
1188
.

气体发生器
,
T。
,
米勒
,
H.G.公司。
,
科勒
,
西。
,
莫利纳里
,
L。
普拉德
,
答:。
(
1984
)
生长曲线的非参数回归分析
.
安。统计师。
,
12
,
210
229
.

胡佛
,
D。
,
大米
,
J。
,
,
C、。
,
左旋-右旋。
(
1998
)
纵向数据时变系数模型的非参数平滑估计
.
生物特征
,
85
,
809
822
.

膝盖
,
答:。
气体发生器
,
T。
(
1992
)
用于分析代表曲线样本的数据的统计工具
.
安。统计师。
,
20
,
1266
1305
.

,
英国。
(
1991
)
降维的切片逆回归
.
美国统计学杂志。助理。
,
86
,
316
327
.

米勒
,
H.G.公司。
(
1984
)
非参数核回归的最优设计
.
统计人员。普罗巴伯。莱特。
,
2
,
285
290
.

米勒
,
H.G.公司。
,
凯里
,
J.R.公司。
,
,
D。
,
列多
,
第页。
沃佩尔
,
J·W·。
(
2001
)
繁殖潜力预测雌性地中海果蝇的寿命
.
程序。R.Soc.伦敦。B类
,
268
,
445
450
.

米勒
,
H.G.公司。
,
,
J·L·。
,
卡普拉
,
W.B.公司。
,
列多
,
第页。
凯里
,
J.R.公司。
(
1997
)
蛋白质缺乏雌性果蝇的早期死亡率激增导致地中海果蝇预期寿命的性别差异发生逆转
.
程序。阿卡德。科学。美国
,
94
,
2762
2765
.

鹧鸪
,
L。
哈维
,
P.H.公司。
(
1985
)
复制成本
.
自然
,
316
,
20
21
.

佩祖利
,
美国。
西尔弗曼
,
B.W.公司。
(
1993
)
函数数据平滑主成分分析的一些性质
.
计算。统计师。
,
8
,
1
16
.

拉姆齐
,
J.O.公司。
西尔弗曼
,
B.W.公司。
(
1997
)
功能数据分析
.
纽约
:
施普林格
.

拉姆齐
,
J.O.公司。
西尔弗曼
,
B.W.公司。
(
2002
)
应用功能数据分析
.
纽约
:
施普林格
.

大米
,
J.A.公司。
西尔弗曼
,
B.W.公司。
(
1991
)
当数据为曲线时,非参数估计均值和协方差结构
.
J.R.统计。Soc.B公司
,
53
,
233
243
.

大米
,
J.A.公司。
,
C、。
(
2001
)
非等采样噪声曲线的非参数混合效应模型
.
生物计量学
,
57
,
253
259
.

里斯
,
F、。
纳吉
,
学士学位。
(
1990
)
功能分析
.
纽约
:
多佛出版物
.

,
M。
,
韦斯
,
R.E.公司。
泰勒
,
J·M·G。
(
1996
)使用灵活的随机曲线分析获得性免疫缺陷综合征的儿科CD4计数。
申请。统计师。
,
45
,
151
163
.

斯坦尼斯瓦利斯
,
J·G·。
,
J·J。
(
1998
)
纵向数据的非参数回归分析
.
美国统计学杂志。助理。
,
93
,
1403
1418
.

附录A:随机过程L(左)2

假设观察到的曲线X(X)1,X(X)2,…,X(X)n个是随机过程的独立且同分布的实现X(X)在域上T型假设过程具有均值和协方差函数

E类{X(X)(t吨)}=μ(t吨),覆盖(cov){X(X)(),X(X)(t吨)}=γ(,t吨).

给定的函数(f),让〈(f),〉表示L(左)2(d)ν)内积,

(f),=T型(f)(t吨)(t吨)d日ν(t吨),

其中dν是上的度量T型通常被选为勒贝格测量。

假设γ到正交本征函数ρk(·)于L(左)2(d)ν)这样的话

γ(,t吨)=k=1λkρk()ρk(t吨),
(17)

具有有序特征值λ1λ2⩾ … ⩾ 0. 特征值是非负的,因为协方差核γ(,t吨)对称且非负定。正交本征函数满足〈ρ,ρ〉=1和〈ρ,ρj个>=0j个,和〈γ(,·),ρk〉 =λk ρk(),T型,k= 1,2,…. 特征值可以表示为

λk=T型T型ρk()γ(,t吨)ρk(t吨)d日ν()d日ν(t吨),

k= 1,2,…. 假设不足(17)然后,随机过程模型为随机选择的曲线提供以下Karhunen–Loève表示:

X(X)(t吨)=μ(t吨)+k=1A类kρk(t吨),
(18)

哪里ρkk特征函数。根据Mercer定理(参见Riesz和Nagy(1990)),右侧方程式(18)一致收敛于t吨T型.随机变量A类k对应于主成分,由下式给出

A类k=ρk,X(X)μ.
(19)

主要组成部分A类k是不相关的随机变量,并且满足

E类A类k=0,无功功率,无功功率A类k=λk,k=1λk<,

k第个特征值对应于k第th个主成分。发件人方程式(18)我们发现(17),X(X)(t吨)由平均函数组成μ加上加性噪声,随机部分。

主要组成部分A类k和基本函数ρk可以解释为定义随机过程关于其平均函数的变化。根据表达式(19),随机变量A类1是的投影长度X(X)μ到上面ρ1对于每个采样曲线X(X)、和A类1ρ1解释了X(X)在所有涉及单个实值随机变量的函数中。同样,函数A类2ρ2解释了无法解释的最大额外过程变化量A类1ρ1,依此类推k= 3,4,….

附录B:定理1的证明

对于方程(10),我们观察到

E类μ˜E类X(X)2=1n个2,j个,k,E类A类kA类j个ρk,ρ=1n个2,kE类A类k2ρk,ρk=O(运行)n个1.

这意味着E类(μ˜μ)=O(运行)n个1/2根据Jensen不等式,因此μ˜-μ=O第页n个-1/2。该断言源自中的属性(b)第3.5节

μ^μμ˜μ啜饮x个=1n个G公司2(x个)=1n个1G公司(x个)0.
(20)

对于方程式(11),让A类k(1)=X(X)μ^,ρ^k,A类k(2)=X(X)μ,ρk,α^k(1)(z(z))=S公司z(z),b条αk,Z轴,A类k(1)=1,,n个α^k(2)(z(z))=S公司z(z),b条αk,Z轴,A类k(2)=1,,n个.注意到E类A类k(2)Z轴=αk(Z轴),一致性假设(a)第3.5节on the smooler意味着α^k(2)(z(z))αk(z(z))=O(运行)第页τn个此外,观察线性假设(b),

α^k(2)(z(z))α^k(1)(z(z))==1n个G公司(z(z))A类k(2)A类k(1)=1n个G公司2(z(z))1/2=1n个A类k(2)A类k(1)21G公司(z(z))01/2.
(21)

请注意

啜饮A类k(2)A类k(1)啜饮X(X),ρkρ^k+μ^,ρ^kρk+ρk,μ^μ=O(运行)第页μ^μ+ρ^kρk,

自从K(K)是有限的,所以supX(X)‖<∞. Then表达式(21)意味着α^k(2)(z(z))α^k(1)(z(z))=O(运行)第页μ^μ+ρ^kρk,使用中的假设(b)第3.5节以与不等式相同的方式(20)。相同参数的第二个应用程序,使用啜饮A类k(1)A类^k=O(运行)第页1的平滑特性ρ^kμ^、和方程式(6)然后暗示结果。

附录C:半参数拟似然回归

SPQR是一种估计未知模型组件的估计程序:这些组件由两个平滑函数组成,即链接函数(·)和方差函数σ2(·)和回归参数向量β0,带║β0║=1.给出观察结果和预测因素x个,模型假设为

E类=x个T型β0,=1,,n个,无功功率,无功功率=无功功率,无功功率ε=σ2x个T型β0=σ2E类.

通过交替参数和非参数估计步骤,提出了一种三阶段迭代QLUE过程。程序可概括如下:;见邱和米勒(1998,1999)了解更多详细信息。

S公司(ν){z(z),小时;(z(z),)= 1,pt(磅)…,n个},ν⩾0,是非参数估计量的通用符号ν回归函数d的th导数νE类(Y(Y)|Z轴=z(z))/d日z(z)ν,基于散点图数据(z(z),)= 1,pt(磅)…,n个。这里是z(z)s是设计点,s是要平滑的原始数据,小时表示带宽和z(z)是要计算函数的目标级别。成为νlink函数的th导数。各种更新步骤如下。

  • 步骤1-链接函数的非参数估计步骤:对于给定β^,β^=1,链接函数及其一阶导数的估计值更新为
    ^ν(t吨;β^)=S公司(ν)t吨,b条ν;x个T型β^,=1,,n个,ν=0,1
  • 第2步-方差函数的非参数估计步骤:对于给定β^,β^=1^0(;β^),更新的非参数方差估计σ^2()通过以下方式获得
    σ^2(u个)=S公司(0)u个,b条;μ^,ε^2=1,,n个
    哪里μ^=^0x个T型β^;β^ε^2=μ^2是平方残差,作为“原始”方差估计,并基于当前模型拟合。
  • 步骤3-参数估计步骤:对于给定的^0(;β^),^1(;β^)σ^2(),β^通过求解关于β,插入链接和方差函数的当前估计值。分数的估计方程为
    U型*(β)==1n个^0η;βσ^2^0η;β^1η;βx个,
    哪里x个=(x个1,x个2,…,x个知识产权)T型,η=x个T型β和║β║=1

迭代一直运行到收敛。分层逆回归(Li,1991)估计值提供了令人满意的初始值β为了简化符号,我们表示回归参数的估计向量和迭代收敛时获得的估计链接函数

{β^,^()}=昆士兰州x个,=1,,n个.
本文根据牛津大学出版社标准期刊出版模式的条款出版和发行(https://academic.oup.com/journals/pages/open_access/funder_policies/chorus/standard_publication_model)