Mixed-Effects模型

混合效应模型介绍

在统计中,一个效果是在预测变量的特定设置下影响响应变量的值的任何东西。效果被转换为模型参数。在线性模型中,效应成为系数,表示模型项的比例贡献。在非线性模型中,效应通常有特定的物理解释,并出现在更一般的非线性组合中。

固定的影响表示总体参数,假设每次收集的数据都是相同的。估计固定效应是回归建模的传统领域。随机效应,则为样本相关的随机变量。在建模中,随机效应就像附加误差项,必须指定它们的分布和协方差。

例如,考虑一个从血液中去除药物的模型。模型使用时间t作为预测药物浓度的指标C作为响应。非线性模型项C0e- - - - - -rt结合参数C0r,分别表示初始浓度和消除速率。如果数据是跨多个个体收集的,那么可以合理地假设消除率是一个随机变量r取决于个人,在总体均值附近变化 r ¯ 。这个词C0e- - - - - -rt就变成了

C 0 e ( r ¯ + ( r r ¯ ) ] t = C 0 e ( β + b ) t ,

在哪里β= r ¯ 是固定效应和b= r r ¯ 是随机效应。

当数据落入自然组时,随机效应很有用。在药物消除模型中,这些群体只是被研究的个体。更复杂的模型可能会根据个体的年龄、体重、饮食等对数据进行分组。虽然这些群体并不是研究的重点,但是在模型中加入随机效应可以将推论的可靠性扩展到个体的特定样本之外。

Mixed-effects模型同时考虑固定和随机效应。与所有回归模型一样,它们的目的是将响应变量描述为预测变量的函数。然而,混合效应模型可以识别样本子组之间的相关性。通过这种方式,它们在完全忽略数据组和用单独的模型对每个组进行拟合之间提供了一种折衷。

Mixed-Effects模型层次结构

假设一个非线性回归模型的数据落入不同的组= 1,…,。要说明模型中的组,请编写responsej在集团为:

y j = f ( φ , x j ) + ε j

yj的反应,xj是一个预测向量,φ为模型参数向量,εj是测量或加工误差。该指数j范围从1到n,在那里n是观察组的个数吗。这个函数f指定模型的形式。通常,xj仅仅是观察时间吗tj。误差通常被认为是独立的、恒等的、正态分布的、方差恒定的。

参数的估计值φ描述总体,假设这些估计对所有群体都是一样的。但是,如果估计值因组而异,则模型变为

y j = f ( φ , x j ) + ε j

在混合效应模型中,φ可能是固定效应和随机效应的组合:

φ = β + b

随机效应b通常被描述为多元正态分布,具有均值为零和协方差Ψ。估计固定效应β随机效应的协方差Ψ提供不假定参数的总体的描述φ组间是相同的。估计随机效应b还提供了数据中特定组的描述。

模型参数不需要用单独的效果来识别。一般来说,设计矩阵一个B是用来识别参数与固定和随机效应的线性组合:

φ = 一个 β + B b

如果组间的设计矩阵不同,则模型变为

φ = 一个 β + B b

如果设计矩阵在不同的观测值之间也不同,模型就变成

φ j = 一个 j β + B j b y j = f ( φ j , x j ) + ε j

一些特定于组的预测器xj可以不随观察而改变吗j。调用这些v,模型变成

y j = f ( φ j , x j , v ) + ε j

指定Mixed-Effects模型

假设一个非线性回归模型的数据落入不同的组= 1,…,。(特别地,假设组没有嵌套。)为该数据指定一个通用的非线性混合效应模型:

  1. 定义特定于组的模型参数φ作为固定效果的线性组合β和随机效应b

  2. 定义响应值y作为一个非线性函数f参数和组特定的预测变量X

该模型是:

φ = 一个 β + B b y = f ( φ , X ) + ε b N ( 0 , Ψ ) ε N ( 0 , σ 2 )

非线性混合效应模型的这个公式使用如下符号:

φ 组特定模型参数的向量
β 一个固定效果的向量,建模的人口参数
b 多元正态分布群特异性随机效应的向量
一个 用于组合固定效果的特定组的设计矩阵
B 一种用于组合随机效果的组专用设计矩阵
X 一组特定预测值的数据矩阵
y 特定于组的响应值的数据向量
f 的一般实值函数φX
ε 组特有误差的向量,被认为是独立的、恒等的、正态分布的和独立的b
Ψ 随机效应的协方差矩阵
σ2 误差方差,假设各观测值之间为常数

例如,考虑一个从血液中去除药物的模型。模型包含两个重叠的阶段:

  • 一个初始阶段p在此期间,药物浓度与周围组织达到平衡

  • 第二阶段在此期间,药物从血液中消失

获取多个个体的数据,模型为

y j = C p e r p t j + C e r t j + ε j ,

在哪里yj观察到的浓度在个体内吗在时间tj。该模型允许对不同的个体进行不同的采样时间和不同数量的观察。

消除率rpr必须是积极的身体有意义。通过引入日志速率来实现这一点Rp=日志(rp),R=日志(r),并重新参数化模型:

y j = C p e 经验值 ( R p ) t j + C e 经验值 ( R ) t j + ε j

在建立混合效应模型时,选择随机效应模型的参数是一个重要的考虑因素。一种技术是向所有参数添加随机效应,并使用对其方差的估计来确定它们在模型中的重要性。另一种方法是将模型分别适合于每一组,不存在随机效应,并观察参数估计的变化。如果估计值在不同的组之间有很大的差异,或者如果每个组的置信区间有最小的重叠,那么该参数就是随机效应的一个很好的候选者。

引入固定效果β和随机效应b对于所有模型参数,将模型重新表述如下:

y j = ( C ¯ p + ( C p C ¯ p ) ] e 经验值 ( R ¯ p + ( R p R ¯ p ) ] t j + ( C ¯ + ( C C ¯ ) ] e 经验值 ( R ¯ + ( R R ¯ ) ] t j + ε j = ( β 1 + b 1 ) e 经验值 ( β 2 + b 2 ) t j + ( β 3. + b 3. ) e 经验值 ( β 4 + b 4 ) t j + ε j

用一般模型表示:

β = ( β 1 β 4 ) , b = ( b 1 b 4 ) , y = ( y 1 y n ) , X = ( t 1 t n ) ,

在哪里n是个体的观察次数吗。在这种情况下,设计矩阵一个B是,至少最初是,4×4的单位矩阵。可以根据需要改变设计矩阵,以引入个别效果或时间依赖性的权重。

拟合模型并估计协方差矩阵Ψ往往会导致进一步的改进。一个相对较小的随机效应方差估计表明,它可以从模型中删除。同样地,对某些随机效应的协方差的相对较小的估计表明,完全的协方差矩阵是不必要的。由于随机效应是不可观测的,Ψ必须间接估计。指定对角或块对角的协方差模式Ψ可以提高拟合算法的收敛性和效率。

统计和机器学习工具箱功能nlmefitnlmefitsa对数据拟合一般非线性混合效应模型,估计固定效应和随机效应。函数也估计协方差矩阵Ψ随机效应。额外的诊断输出允许您评估模型参数的数量和适合度之间的权衡。

指定协变量模型

如果模型指定Mixed-Effects模型假设一个依赖群体的协变量,如权重(w)模型变成:

( φ 1 φ 2 φ 3. ) = ( 1 0 0 0 1 0 0 0 1 w 0 0 ) ( β 1 β 2 β 3. β 4 ) + ( 1 0 0 0 1 0 0 0 1 ) ( b 1 b 2 b 3. )

因此,参数φ中的任何个人th组是:

( φ 1 φ 2 φ 3. ) = ( β 1 + β 4 * w β 2 β 3. ) + ( b 1 b 2 b 3. )

若要指定协变量模型,请使用“FEGroupDesign”选择。

“FEGroupDesign”是一个p-by-q-by-m指定不同的p-by-q为每个固定效果设计矩阵组。使用前面的例子,数组类似于以下:

  1. 创建数组。

    模型参数个数(Phi) num_params = 3;协变量数% num_cov = 1;假设数据集中的组数为7,num_groups = 7;数组的协变量值协变量= [75;52个;66;55;70;58;62); A = repmat(eye(num_params, num_params+num_cov),... [1,1,num_groups]); A(1,num_params+1,1:num_groups) = covariates(:,1)
  2. 创建一个结构体使用指定的设计矩阵。

    选项。FEGroupDesign =一个;
  3. 指定nlmefit(或nlmefitsa),如使用nlmefitsa和nlmefitsa的混合效果模型

选择nlmefitnlmefitsa

统计和机器学习工具箱提供了两个功能,nlmefitnlmefitsa用于拟合非线性混合效应模型。每个函数提供不同的功能,这些功能可以帮助您决定使用哪个功能。

近似的方法

nlmefit为拟合非线性混合效应模型提供了以下四种近似方法:

  • LME的-使用线性混合效应模型在当前条件估计的可能性βB。这是默认值。

  • “RELME”-在当前的条件估计下,使用线性混合效应模型的限制似然βB

  • ‘佛’-无随机效应的一阶拉普拉斯近似。

  • “FOCE”的条件估计的一阶拉普拉斯近似B

nlmefitsa提供一种附加的逼近方法,随机逼近期望最大化(SAEM)[24]分三步:

  1. 模拟:生成随机效应的模拟值b从后验密度p(b|Σ)鉴于目前的参数估计。

  2. 随机逼近:更新日志似然函数的期望值,从前面的步骤中获取其值,并向模拟随机效应计算的日志似然的平均值移动一部分。

  3. 最大化步骤:在给定随机效应的模拟值的情况下,选择新的参数估计来最大化对数似然函数。

这两个nlmefitnlmefitsa尝试寻找参数估计以使似然函数最大化,这是很难计算的。nlmefit通过各种方法逼近似然函数,并使近似函数最大化来处理该问题。它使用传统的优化技术,这些技术依赖于收敛标准和迭代限制。

nlmefitsa另一方面,以这样一种方式模拟参数的随机值:从长远来看,它们收敛于使精确似然函数最大化的值。结果是随机的,传统的收敛测试不适用。因此nlmefitsa提供选项以在模拟过程中绘制结果,并多次重新启动模拟。您可以使用这些特性来判断结果是否已收敛到您希望的精度。

参数特定于nlmefitsa

以下参数是特定于nlmefitsa。大多数控制随机算法。

  • Cov0-协方差矩阵的初始值ψ。必须是一个r——- - - - - -r正定矩阵。如果为空,则默认值取决于的值BETA0

  • ComputeStdErrors- - - - - -真正的计算系数估计值的标准误差并将其存储在输出中统计数据结构,或者(默认)省略此计算。

  • LogLikMethod-指定逼近日志似然值的方法。

  • NBurnIn-不重新计算参数估计的初始老化迭代次数。默认是5。

  • NIterations-控制算法的三个阶段每个阶段执行多少次迭代。

  • NMCMCIterations-马尔科夫链蒙特卡罗(MCMC)迭代次数。

模型和数据需求

的能力有一些不同nlmefitnlmefitsa。因此,有些数据和模型对于任何一个函数都是可用的,但是有些可能要求您只选择其中的一个函数。

  • 误差模型- - - - - -nlmefitsa万博1manbetx支持各种错误模型。例如,响应的标准差可以是常数,与函数值成比例,或者两者的组合。nlmefit在响应的标准差为常数的假设下拟合模型。其中一个误差模型,“指数”,指定响应的日志具有恒定的标准偏差。您可以使用以下工具来适应这些模型nlmefit通过提供日志响应作为输入,并通过重写模型函数来生成非线性函数值的日志。

  • 随机效应-两个函数都将数据与带有参数的非线性函数相匹配,参数可以是简单的标量值,也可以是协变量的线性函数。nlmefit允许线性函数的任意系数同时具有固定和随机效应。nlmefitsa万博1manbetx仅对线性函数的常数(截距)系数支持随机效应,而对斜率系数不支持。在in的例子中指定协变量模型,nlmefitsa只能将前三个beta值视为随机效应。

  • 模型形式- - - - - -nlmefit万博1manbetx支持非常通用的模型规范,对将固定系数和随机效应与模型参数关联起来的设计矩阵几乎没有限制。nlmefitsa更严格的:

    • 固定效应设计在每个群体(对于每个个体)中必须是恒定的,因此不支持依赖于观测的设计。万博1manbetx

    • 随机效应设计对于整个数据集必须是常数,因此既不支持观测相关设计,也不支持群体相关设计。万博1manbetx

    • 正如前面提到的下随机效应,在随机效应设计中,坡度系数不能指定随机效应。这意味着设计必须由0和1组成。

    • 随机效应设计不能对多个系数使用相同的随机效应,也不能对单个系数使用多个随机效应。

    • 固定效应设计对于多个参数不能使用相同的系数。这意味着它在每个列中最多只能有一个非零值。

    如果你想用nlmefitsa对于协变量效应是随机的数据,将协变量直接包含在非线性模型表达式中。不包括在固定或随机效应设计矩阵的协变量。

  • 收敛-如所述模型形式,nlmefitnlmefitsa有不同的方法来测量收敛性。nlmefit采用传统的优化措施,并且nlmefitsa提供诊断以帮助您判断随机模拟的收敛性。

在实践中,nlmefitsa倾向于更健壮,在困难的问题上更不容易失败。然而,nlmefit可能会更快地收敛于它所收敛的问题。一些问题可能从组合策略中获益,例如通过运行nlmefitsa一段时间内获得合理的参数估计值,并使用这些值作为额外迭代的起点nlmefit

使用混合效果模型的输出函数

Outputfcn场的选项结构指定求解器在每次迭代后调用的一个或多个函数。通常,可以使用输出函数在每次迭代中绘制点,或者显示算法的优化量。设置一个输出函数:

  1. 把输出函数写成MATLAB®文件函数或本地函数。

  2. 使用statset设置的值Outputfcn要成为函数句柄,即函数的名称前面加上@的迹象。例如,如果输出函数是outfun.m,命令

    options = statset('OutputFcn', @outfun);

    指定OutputFcn成为…的把柄outfun。要指定多个输出函数,请使用以下语法:

    options = statset('OutputFcn',{@outfun, @outfun2});
  3. 调用优化函数选项作为输入参数。

有关输出函数的示例,请参见样例输出函数

输出函数的结构

输出函数的函数定义行如下:

停止= outfun (β,状态,状态)
在哪里

  • β是当前固定效果。

  • 状态是一个包含当前迭代数据的结构。领域的地位详细描述结构。

  • 状态是算法的当前状态。算法状态列出可能的值。

  • 停止国旗是吗真正的取决于优化例程应该退出还是继续。看到停止标志为更多的信息。

求解程序将输入参数的值传递给outfun在每一个迭代。

领域的地位

下表列出了的字段状态结构:

描述
过程
  • “ALT”-交替优化算法的线性混合效果或限制线性混合效果近似

  • “一圈”-优化一阶或一阶条件估计的拉普拉斯近似

迭代 从0开始的整数。
内心的

结构中描述内部迭代状态的结构ALT腿上程序和字段:

  • 过程——当过程“ALT”:

    • “PNLS”(惩罚非线性最小二乘法)

    • LME的(线性mixed-effects估计)

    • “没有”

    过程“一圈”,

    • “PNLS”(惩罚非线性最小二乘法)

    • “PLM”(异形可能性最大化)

    • “没有”

  • 状态-以下其中一项:

    • “init”

    • “通路”

    • “完成”

    • “没有”

  • 迭代-从0开始的整数,或。为nlmefitsa对于老化迭代,输出函数在每个迭代之后调用,且为负值STATUS.iteration

fval 当前对数似然
ψ 当前随机效应协方差矩阵
θ 的当前参数化ψ
均方误差 当前误差方差

算法状态

下表列出了可能的值状态:

状态 描述

“init”

算法在第一次迭代前处于初始状态。

“通路”

算法在迭代的末尾。

“完成”

在最后一次迭代后,算法处于最终状态。

下面的代码演示了输出函数如何使用的值状态决定在当前迭代中执行哪些任务:

切换状态大小写“iter”%根据需要更新plot或gui,大小写“init”%设置为plot,大小写“done”%清除plot、gui或final plot,否则结束

停止标志

输出参数停止国旗是吗真正的。该标志告诉解决程序应该退出还是继续。下面的示例展示了使用停止国旗。

停止基于中间结果的优化输出函数可以根据传入参数的值在任何迭代中停止评估。例如,下面的代码基于存储在“fval”状态结构域:

停止= outfun(beta,status,state)停止= false;%检查loglikelihood是否大于132。如果状态。fval> -132 stop = true; end

停止基于GUI输入的迭代如果您设计一个GUI来执行nlmefit迭代时,可以在用户单击a时停止输出函数停止GUI上的按钮。例如,下面的代码实现了一个对话框来取消计算:

函数retval = stop_outfcn(beta,str,status)如果isequal(str.inner.state,'none')开关(status)大小写'init' % Initialize dialog stop = false;h = msgbox('按停止键取消计算',…'NLMEFIT: Iteration 0 ');按钮= findobj (h,“类型”,“uicontrol”);set(button,'String','STOP','Callback',@stopper) pos = get(h,'Position');pos(3) = 1.1 * pos(3);在对话框标题集中显示迭代次数(h,'Name',sprintf('NLMEFIT: iteration %d',…str.iteration drawnow));大小写‘完成’%删除对话框删除(h); end end if stop % Stop if the dialog button has been pressed delete(h) end retval = stop; function stopper(varargin) % Set flag to stop when button is pressed stop = true; disp('Calculation stopped.') end end

样例输出函数

nmlefitoutputfcn是样本统计和机器学习工具箱的输出函数吗nlmefitnlmefitsa。它使用固定效果(β)和随机效应的方差(诊断接头(STATUS.Psi))。为nlmefit,图中还包括对数似然(STATUS.fval)。

nlmefitoutputfcn是默认的输出函数吗nlmefitsa。用在nlmefit,在选项结构中指定函数句柄:

= statset('OutputFcn', @nlmefitoutputfcn,…)beta = nlmefit(…,'Options', opt,…)
为了防止nlmefitsa在使用这个函数时,为输出函数指定一个空值:
选择= statset (“OutputFcn”,[],…)β= nlmefitsa(…,“选项”,选择,…)
nlmefitoutputfcn停止nlmefitnlmefitsa如果你关闭它产生的数字。