美文网首页遗传改良（育种）
Isik2011 GLMM

Isik2011 GLMM

作者: 董八七 | 来源:发表于2019-04-02 17:43 被阅读0次

    Isik F (2011) Generalized Linear Mixed Models An Introduction for Tree Breeders and Pathologists.

    前言

    广义线性混合模型（GLMM）在过去几年引起了相当多的关注。 “广义”一词指的是响应变量的非正态分布，而“混合”一词指的是除了回归分析的通常固定效应之外的随机效应。随着SAS，R和ASReml等现代统计软件包的发展，可以为更多的受众提供大量的统计分析。然而，除了能够处理更复杂的模型之外，用户还有责任了解这些高级工具的工作原理。

    本次研讨会的目的是通过首先讨论统计线性模型的一些假设和不足，然后给出自然科学中常见情况的用例，介绍广义线性混合模型。

    • 第一部分回顾了简单和多变量的线性模型和回归分析。使用SAS REG软件解决了两个数值示例。

    • 第二部分通过将随机效应添加到线性模型来呈现线性混合模型。使用SAS MIXED过程提供了一个简单的数值示例。

    • 第三（最后）部分介绍了广义线性模型。使用SAS GLIMMIX过程和ASReml软件呈现二进制和计数数据的两个说明性示例。

    线性模型

    线性模型（回归）通常用于建模单个变量y（称为响应或因变量 response or dependent）与一个或多个预测变量，独立或解释变量（predictor, independent or explanatory）X1,...,Xp之间的关系。当p = 1时，它被称为简单回归，但当p> 1时，它被称为多元回归。

    回归分析可用于评估响应变量上的解释变量之间的关系。它也是预测未来观测或仅描述数据结构的有用工具。

    首先从一个简单的例子开始，假设y是树的重量，预测变量是高度（X1）和树的年龄（X2）。通常，数据将以如下阵列的形式提供
    \begin{array}{ccc}{y_{1}} & {x_{11}} & {x_{12}} \\ {y_{2}} & {x_{21}} & {x_{22}} \\ {\vdots} & {\vdots} & {\vdots} \\ {y_{n}} & {x_{n 1}} & {x_{n 2}}\end{array}

    其中yi是第i个树的观察值，n是观察数。
    有无数种方法来模拟响应和解释变量之间的关系。 但是，为了简单起见，可以通过参数中的线性函数对关系进行建模，如下所示
    y=\beta_{0}+\beta_{1} X_{1}+\beta_{2} X_{2}+\varepsilon
    其中\beta_{i}是未知参数，\varepsilon是残差项。因此，问题就简化为三个参数的估计。
    请注意，在线性模型中，参数线性输入，但预测变量不一定必须是线性的。 例如，考虑以下两个函数：
    y=\beta_{0}+\beta_{1} \exp \left(X_{1}\right)+\beta_{2} \log \left(X_{2}\right)+\varepsilon
    y=\beta_{0}+X_{1}^{\beta_{1}}+X_{2} \exp \left(\beta_{2}\right)+\varepsilon
    第一个是线性，但第二个不是。
    使用矩阵表示，上述示例的回归方程可写为：
    \mathbf{y}=\mathbf{X} \boldsymbol{\beta}+\boldsymbol{\varepsilon}
    \mathbf{y}=\left(y_{1}, \dots, y_{n}\right)^{T}，\boldsymbol{\beta}=\left(\beta_{0}, \beta_{1}, \beta_{2}\right)^{T}，\boldsymbol{\varepsilon}=\left(\varepsilon_{0}, \dots, \varepsilon_{n}\right)^{T}，设计矩阵X是
    \mathbf{X}=\left[ \begin{array}{ccc}{1} & {x_{11}} & {x_{12}} \\ {1} & {x_{21}} & {x_{22}} \\ {\vdots} & {\vdots} & {\vdots} \\ {1} & {x_{n 1}} & {x_{n 2}}\end{array}\right]
    可以使用最小二乘法来进行\beta的估计。 也就是说，我们将\hat{\beta}定义为\beta的最佳估计值，即最小化误差平方和的意义。
    \sum_{i=1}^{n} \varepsilon_{i}^{2}=\boldsymbol{\varepsilon}^{\mathrm{T}} \boldsymbol{\varepsilon}=(\mathbf{y}-\mathbf{X} \boldsymbol{\beta})^{T}(\mathbf{y}-\mathbf{X} \boldsymbol{\beta})
    对\beta微分，并令其等于零，可以表明\hat{\beta}满足正态方程：
    \mathbf{X}^{\mathrm{T}} \mathbf{X} \widehat{\boldsymbol{\beta}}=\mathbf{X}^{\mathbf{T}} \mathbf{y}
    因此，令\mathbf{X}^{\mathrm{T}} \mathbf{X}可逆，
    \widehat{\boldsymbol{\beta}}=\left(\mathbf{X}^{\mathrm{T}} \mathbf{X}\right)^{-1} \mathbf{X}^{\mathbf{T}} \mathbf{y}
    到目前为止，我们还没有假设任何残差\varepsilon的分布形式。 通常的假设残差是正态分布，实际上这通常但并非完全是一个合理的假设。
    如果我们假设残差是独立的并且相同的正态分布，均值为0，方差为\sigma^2，也就是说\boldsymbol{\varepsilon} \sim N\left(0, \sigma^{2} \mathbf{I}\right)，然后观测值的期望是
    \mathbf{y} \sim N\left(\mathbf{X} \boldsymbol{\beta}, \sigma^{2} \mathbf{I}\right)
    和参数的期望是
    \widehat{\boldsymbol{\beta}} \sim N\left(\boldsymbol{\beta},\left(\mathbf{X}^{\mathrm{T}} \mathbf{X}\right)^{-1} \sigma^{2}\right)

    线性回归示例（略）

    我们想探索身高和体重之间的线性依赖关系。 线性模型可以使用SAS GLM或REG程序运行。

    多元线性回归（略）

    现在考虑我们还想将树龄作为另一个预测因子。 也就是说，我们想要评估年龄和体重之间的关系，同时保持身高的效果不变。 在简单回归中，直线适合数据，而在多个回归中，要拟合第p维平面，其中p是解释变量的数量。 为了可视化数据和拟合值，需要3D图。

    线性混合模型

    线性混合模型是包含固定效应和随机效应的统计模型。 这些模型广泛用于生物和社会科学。 在矩阵表示法中，线性混合模型可以表示为
    \mathbf{y}=\mathbf{X} \boldsymbol{\beta}+\mathbf{Z} \mathbf{\gamma}+\boldsymbol{\varepsilon}

    y是观测的nx1向量
    \beta是px1固定效应向量
    \gamma是qx1随机效应向量，
    \varepsilon是nx1随机误差项向量
    X是nxp固定效应设计矩阵，将观测值观测y关联到\beta Z是nxq随机效应设计矩阵，将观测值观测y关联到\gamma

    我们将假设\gamma和\varepsilon是不相关的随机变量，分别具有零均值和协方差矩阵G和R。
    \begin{array}{ll}{E[\gamma]=0,} & {\operatorname{Var}[\gamma]=G} \\ {E[\varepsilon]=0,} & {\operatorname{Var}[\varepsilon]=R}\end{array}
    \operatorname{cov}(\varepsilon, Y)=0
    因此，观测向量y的期望和方差（V）由下式给出：
    \begin{aligned} E[\mathbf{y}] &=\mathbf{X} \boldsymbol{\beta} \\ \operatorname{Var}[\mathbf{y}]=\mathbf{V} &=\mathbf{Z G Z}^{\mathrm{T}}+\mathbf{R} \end{aligned}
    理解V矩阵是使用混合模型的一个非常重要的组成部分，因为它包含两个随机变量源，并定义这些模型与普通最小二乘（OLS）计算的区别。
    如果只有随机效应模型（例如随机区组设计），则G矩阵是主要焦点。 另一方面，对于重复测量或用于或空间分析，R矩阵是相关的。

    如果我们也假设随机项是正态分布：
    \boldsymbol{\gamma} \sim N(\mathbf{0}, \mathbf{G}), \quad \boldsymbol{\varepsilon} \sim N(\mathbf{0}, \mathbf{R})
    然后，观测向量将为正态分布\mathbf{y} \sim N(\mathbf{X} \boldsymbol{\beta}, \mathbf{V})。
    对于上述一般线性混合模型，Henderson的混合模型方程（MME）可用于找到\hat{\beta}和\hat{\gamma}，分别为\beta的最佳线性无偏估计（BLUE）和\gamma的最佳线性无偏预测（BLUP）。
    \left[ \begin{array}{c} X'R^{-1}X & X'R^{-1}Z\\ Z'R^{-1}X & Z'R^{-1}Z+G^{-1} \end{array} \right] \left[ \begin{array}{c} \hat{\beta}\\ \hat{\gamma} \end{array} \right] = \left[ \begin{array}{c} X'R^{-1}y\\ Z'R^{-1}y \end{array} \right]
    模型的解也可以写成：
    \begin{array}{c}{\widehat{\boldsymbol{\beta}}=\left(\mathbf{X}^{\mathbf{T}} \mathbf{V}^{-1} \mathbf{X}\right)^{-1} \mathbf{X}^{\mathbf{T}} \mathbf{V}^{-1} \mathbf{y}} \\ {\hat{\mathbf{\gamma}}=\mathbf{G} \mathbf{Z}^{\mathrm{T}} \mathbf{V}^{-1}(\mathbf{y}-\mathbf{X} \widehat{\boldsymbol{\beta}})}\end{array}
    如果G和R矩阵是已知的，则广义最小二乘可以估计固定效应\beta的任何线性组合。 然而，通常情况下这些矩阵是未知的，因此必须使用用于拟合线性模型的复杂迭代算法来估计它们。

    请考虑以下示例。 假设我们收集了在两个不同地点测量的不同树木生长的数据。 我们可以假设树木来自大群体，这是一个合理的假设，因此，我们将它们视为随机的。

    Tree (t) Location (l) Height (y)
    1 1 87
    2 2 84
    3 2 75
    4 1 90
    5 2 79

    线性混合模型可以用矩阵表示法表示如下
    \left[ \begin{array}{c}{87} \\ {84} \\ {75} \\ {90} \\ {79}\end{array}\right]=\left[ \begin{array}{cc}{1} & {0} \\ {0} & {1} \\ {0} & {1} \\ {1} & {0} \\ {0} & {1}\end{array}\right] \left[ \begin{array}{l}{l_{1}} \\ {l_{2}}\end{array}\right]+\left[ \begin{array}{ccccc}{1} & {0} & {0} & {0} & {0} \\ {0} & {1} & {0} & {0} & {0} \\ {0} & {0} & {1} & {0} & {0} \\ {0} & {0} & {0} & {1} & {0} \\ {0} & {0} & {0} & {0} & {1}\end{array}\right] \left[ \begin{array}{c}{t_{1}} \\ {t_{2}} \\ {t_{3}} \\ {t_{4}} \\ {t_{5}}\end{array}\right]+\left[ \begin{array}{c}{\varepsilon_{1}} \\ {\varepsilon_{2}} \\ {\varepsilon_{3}} \\ {\varepsilon_{4}} \\ {\varepsilon_{5}}\end{array}\right]
    具有矩阵形式\mathbf{y}=\mathbf{X} \boldsymbol{\beta}+\mathbf{Z} \mathbf{\gamma}+\boldsymbol{\varepsilon}。请注意，我们将地点视为固定效应。 矩阵X和Z将表型观测值与地点和随机树效应联系起来。

    我们需要对模型的方差分量做一些假设。 在这种情况下，我们假设树彼此独立，并且残差是独立的。 因此，方差矩阵的结构可表示为：
    \mathbf{R}=\mathbf{I}_{\mathrm{n}} \sigma_{\varepsilon}^{2} \quad \mathbf{G}=\mathbf{I}_{\mathrm{n}} \sigma_{\gamma}^{2} \quad \mathbf{V}=\mathbf{Z G Z}^{\mathrm{T}}+\mathbf{R}
    \mathbf{I}_{\mathrm{n}}是n x n的单位矩阵。此外，我们假设：
    \boldsymbol{\gamma} \sim N(\mathbf{0}, \mathbf{G}) \quad \boldsymbol{\varepsilon} \sim N(\mathbf{0}, \mathbf{R}) \quad \mathbf{Y} \sim N(\mathbf{X} \boldsymbol{\beta}, \mathbf{V})
    我们可以使用R软件或使用矩阵代数的任何其他软件来计算模型解。

    广义线性模型

    回想一下，线性模型，正如其名称所述，有如下假设：

    • 因变量和固定效应之间的关系可以通过线性函数建模，
    • 方差不是均值的函数，
    • 随机项遵循正态分布
      在某些情况下，可能违反任何或甚至所有这些假设。
      广义线性模型（GLM）是普通最小二乘回归的扩展。 GLM一般化线性回归，通过1.允许通过链接函数把线性模型与响应变量关联，并且2.允许每个测量的方差的大小是其预测值的函数。
      在GLM中，假设因变量y的每个结果是来自从指数族中的特定分布。该族最常见的分布是二项式，泊松和正态分布。
      分布的均值依赖于自变量X，通过反向链接函数（g^{-1}）。
      E(\mathbf{Y})=\boldsymbol{\mu}=g^{-1}(\mathbf{X} \boldsymbol{\beta})=g^{-1}(\mathbf{\eta})
      E(\mathbf{Y})是y的期望，\eta=\mathbf{X} \boldsymbol{\beta}是线性预测子（linear predictor），是未知参数\beta的线性组合。
      在此框架中，方差V通常是均值的函数：
      \operatorname{Var}(\mathbf{Y})=V(\boldsymbol{\mu})=V\left(g^{-1}(\mathbf{X} \boldsymbol{\beta})\right)=V\left(g^{-1}(\mathbf{\eta})\right)
      作为广义线性模型的一个例子，考虑一个回归分析，其中响应变量遵循伯努利分布（结果分别为1或0，概率为p和1-p），我们想要估计某个事件的概率p。
      让我们关注线性模型假设，这对于二元响应变量可能是不合理的。
    • 首先明确违反了响应变量正态分布的假设。
    • 其次，不满足均值和方差之间的恒定方差和独立性的假设。例如，预测概率p给出方差p（1-p），其明显取决于平均值，因此是非常数的。
    • 平均响应不一定是效应值的线性函数。在许多实际情况中，线性假设并不成立。
    • 最后，也许是最重要的，在线性模型中，预测可以取任何值，而对于某些响应变量，我们希望在给定边界之间进行预测（例如，0和1用于疾病）。
      对于所有上述事实，线性模型不是用于分析二元特征的合适工具。广义线性模型解决了线性模型的不足。
      考虑将单个种子放在不同盆中的潮湿土壤上的情况。在该实验中，盆在不同温度T下保持数天D。在任意天数之后，检查种子并且如果其已经发芽则记录结果y = 1，否则y = 0。发芽概率p可以通过以下形式的线性函数建模：
      \eta=\beta_{0}+\beta_{1} D+\beta_{2} T
      其中\eta是线性预测子，\beta_{i0}是要估计的参数。 请注意，线性预测子没有任何值保持在0和1之间。
      但是，在GLM中，概率通过反向链接函数限制在区间（0,1）：
      p=\frac{1}{1+e^{-\eta}}=g^{-1}(\eta)
      值得注意的是，p是二项分布的y的期望值。 还要注意p和线性预测子间的非线性关系由反向链接函数建模。 在这种特殊情况下，链接函数是逻辑logistic链接函数或logit：
      \eta=\log \left(\frac{p}{1-p}\right)=g(p)
      反向链接函数是p=[1+\exp (-\eta)]^{-1}，然后，线性预测子采用以下形式：\eta=-8+0.19 T+0.37 D。使用该模型，我们可以估计在给定温度下发病概率超过80％所需的天数。 在一些简单的代数之后，可以证明在温度10下需要至少20天才能达到0.80发芽的概率。
      以上示例没有随机效应，因此它是广义线性模型（GLM），而不是广义混合模型（GLMM）。

    计数数据示例 - 感染的树数

    我们之前已经考虑了结果变量是数值和二项的示例。 结果变量通常是数字但是以计数的形式。 有时它是一系列罕见事件，例如在一段时间内在人口中发生的火炬松梭形锈病的新病例数。
    泊松概率函数适用于没有上限范围的计数数据，即y可以是任何非负整数。泊松pmf由函数描述：
    \operatorname{Pr}(y=k)=\frac{\lambda^{k} e^{-\lambda}}{k !}, \quad k=0,1,2, \dots
    对于这样的分布，lambda是平均值，方差也是lambda。 这个lambda通常被称为泊松“强度”，因为它是平均事件计数。
    作为一个假设的例子，我们可以使用与三个月内新感染的树木数量相关的数据，其中包括年龄A和身高H。
    我们的兴趣在于将lambda建模为给定家庭和位置的年龄（A）和高度（H）的函数。 使用线性模型可能导致负强度，lambda<0，这是没有意义的。 Poisson的自然链接函数g(\lambda)是对数，因此模型可以如下：
    \eta=\beta_{0}+\beta_{1} A+\beta_{2} H
    链接函数及其反函数由下式给出：
    \begin{array}{l}{\eta=\ln (\lambda)=g(\lambda)} \\ {\lambda=e^{\eta}=g^{-1}(\eta)}\end{array}

    广义线性混合模型

    广义线性混合模型（GLMM）是广义线性模型（GLM）的扩展，其中线性预测子除了固定效应之外还包含随机效应。
    GLMM的期望是：
    E[\mathbf{y} | \mathbf{u}]=g^{-1}(\mathbf{X} \boldsymbol{\beta}+\mathbf{Z u})=g^{-1}(\mathbf{\eta})
    将固定和随机效应组合以形成线性预测子：
    \eta=X \beta+Z u
    通过添加残差向量获得观测向量y的模型，如下：
    \mathbf{Y}=\mathbf{\eta}+\boldsymbol{\varepsilon}=\mathbf{X} \boldsymbol{\beta}+\mathbf{Z u}+\boldsymbol{\varepsilon}
    线性预测子和观测向量之间的关系被建模为
    \mathbf{y}\left|\mathbf{u} \sim\left(g^{-1}(\mathbf{\eta}), \mathbf{R}\right)\right.
    上述表示法表示，给定u，y的条件分布具有均值g^{-1}(\mathbf{\eta})和方差R。条件分布\mathbf{y}|\mathbf{u}通常称为误差分布。 请注意，，是给条件响应\mathbf{y}|\mathbf{u}指定分布，而不是像GLM中一样指定y的分布。该公式也称为条件模型规范。
    最后，观测值的方差矩阵由下式给出：
    \mathrm{V}(\mathbf{y})=\mathrm{E}[\mathrm{V}(\mathbf{y} | \mathbf{u})]+\mathrm{V}[\mathrm{E}(\mathbf{y} | \mathbf{u})]=\mathbf{A}^{1 / 2} \mathbf{R} \mathbf{A}^{1 / 2}+\mathbf{Z G Z}^{\mathrm{T}}
    其中矩阵A是包含模型的方差函数的对角矩阵。
    广义线性混合模型包含几种重要类型的统计模型。 例如，

    • 线性模型：无随机效应，一致identity链接函数和正态分布
    • 广义线性模型：不存在随机效应
    • 线性混合模型：随机效应，一致链接函数和正态分布
      已经开发了广义线性混合模型以解决线性混合模型的不足。
      在许多情况下，隐含的假设不合适。
      例如，线性混合模型假设因变量y的平均值与固定和随机效应之间的关系可以通过线性函数建模。例如，在模拟疾病发病率方面，这种假设是值得怀疑的。
      线性混合模型的另一个假设是方差不是均值的函数，随机效应遵循正态分布。当分析0/1性状时，例如患病（1）或不患病（0），违反了恒定方差的假设。在这种情况下，响应变量是二项式。因此，对于预测的疾病发病率，方差是平均值的函数。
      正态性假设对二元性状无效。一个随机变量，只能取两个值，零或1。相反，正态分布是钟形曲线，可以取任何实数。
      最后，线性混合模型的预测可以取任何值。而二元变量的预测是有界的（0,1）或对于计数变量，它不能取负值。

    概率，Odds和Odds Ratio

    为了理解输出并正确解释结果，我们需要熟悉概率，Odds和Odds Ratio。示例：

    No Yes Sum
    Sample1 30 70 100
    Sample2 20 180 200
    Sum 50 250 300

    样本1：Yes的概率=70/100=0.7，No的概率=30/100=0.3。
    odds=p/(1-p)=0.7/0.3=2.3。我们预计在例1中出现次数是非出现次数的2.3倍。
    样本2：Yes的概率=180/200=0.9，No的概率=20/200=0.1。
    odds=p/(1-p)=0.9/0.1=9。我们预计在例2中出现次数是非出现次数的9倍。
    样本2对样本1的odds比率是OR=9/2.3=3.86。
    在样本2中获得结果的odds（是）几乎是样本1中结果的4倍。

    二项和泊松回归模型中的过度分散

    当数据比某些参考模型下预期的更分散时，会产生过度分散的结果。由于二项或泊松回归模型分析的计数数据的方差是平均值的函数，因此可能会出现这种情况。即，
    \operatorname{Var}[\mathrm{Y}]=\mathrm{f}(\mathrm{E}[\mathrm{Y}]) * \varphi
    对于这两种分布，尺度参数phi的值都是1。
    为了理解过度分散意味着什么，首先回顾线性回归模型（用普通最小二乘法计算）。在正态分布下，由于均值和方差不相关，数据不会过度分散。线性模型（\mathbf{y}=\mathbf{X} \boldsymbol{\beta}^\prime+\boldsymbol{\varepsilon}）的期望值为
    \mathbf{e} \sim \mathrm{NID}\left(0, \sigma^{2}_\mathrm{e}\right)
    残差方差（\sigma^{2}_{e}）被假定为所有协变量线性组合的常数。这个方差是根据数据估计的，可以假设任何大于零的值，不管平均值是什么。因此，假设响应值具有恒定的方差：
    Var(y)=\sigma^{2}_{e}*1
    正态误差和同一identity连接函数（线性回归）的方差函数var(\mu)=1。这个方差对所有的y_i都是常数。
    对于广义线性模型：
    g(\mu)=X^{\prime} \beta
    \mu=E(y)，g是连接函数。y的方差是：
    Var(\mathrm{y})=\varphi * \mathrm{V}(\mu)
    即，观测值的方差等于某个常数（尺度参数）乘以y平均值的函数。

    • 对于二元变量y，方差是其均值的乘法函数：Var(y)=- \mu(1-\mu)
    • 在泊松分布下，方差是均值本身Var(y)=\mathrm{E}[y]=\mu
      在这两种情况下，观察到的计数都有方差，这些方差是依赖于平均值的函数。也就是说，y的方差取决于y的期望值，这是根据数据估计的。
      当任一模型在假设数据是由二项式分布或泊松过程生成的情况下拟合时，尺度参数phi自动设置为1。这就是为什么当我们拟合泊松分布或二项式分布时，asreml输出或sas genmod过程输出中的误差方差为1.00的原因。值1不是误差方差，但它是一个尺度参数，不应作为方差分量来计算二项分布的遗传力。
      对于二项式和泊松回归模型，在所选模型适当的前提下，估计协方差矩阵（以及参数估计的标准误差）。数据的变化可能比分布假设预期的要多。这被称为过度分散（也称为异质性），通常发生在观测结果相关或从“群clusters”中收集的时候。
      为了确定给定模型数据中可能存在的过度分散，将偏差除以其自由度：这称为分散参数。如果偏差合理地“接近”自由度（即尺度参数=1），则缺乏过度分散的证据。
      分散参数（或尺度偏差）=偏差/df
      大于1的尺度参数并不一定意味着存在过度分散。这也可能指示其他问题，例如不正确指定的模型（省略的变量、交互作用或非线性项）、不正确指定的函数形式（可适当使用加法而不是乘法模型）以及影响或异常的观察结果。
      如果您认为您已经正确地指定了模型，并且比例估计值大于1，那么就可以得出您的数据过于分散的结论。您应该能够确定数据过度分散的可能原因。如果不纠正过度分散，则标准误差的估计值太小，从而导致有偏差的推断（即，您将观察到比您应该观察到的更小的p值，从而产生更多的I型误差）。因此，置信区间也不正确。
      当你有“正确”的模型时，异常值不是问题，而且尺度偏差很大，SAS程序（genmod、glimmix、nlmixed）或ASREML有多种选择来纠正过度分散。

    例1：随机分组中的二项式计数

    （修改自SAS Glimmix帮助系统）
    在广义线性模型空间预测的背景下，Gotway和Stroup（1997）分析了一项农艺学田间试验的数据。研究人员研究了16个小麦品种（条目）对麻蝇侵染的抗性。他们在8x8网格上以随机的完整块设计排列这些品种。该布局的每个4x4象限构成一个块。感兴趣的结果是受损植物的数量（Y_{ij}）占该单元上生长的植物总数（n_{ij}）。换句话说，确定条目（植物品种）之间是否存在显著差异。两个下标标识块（i=1，…，4）和条目（j=1，…，16）。
    以下SAS语句创建数据集。变量lat和lng表示8x8网格上实验单元的坐标。

    data HessianFly;
label Y = 'No. of damaged plants'
n = 'No. of plants';
input block entry lat lng n Y @@;
datalines;
1 14 1 1 8 2 1 16 1 2 9 1
1 7 1 3 13 9 1 6 1 4 9 9
1 13 2 1 9 2 1 15 2 2 14 7
1 8 2 3 8 6 1 5 2 4 11 8
1 11 3 1 12 7 1 12 3 2 11 8
1 2 3 3 10 8 1 3 3 4 12 5
1 10 4 1 9 7 1 9 4 2 15 8
1 4 4 3 19 6 1 1 4 4 8 7
2 15 5 1 15 6 2 3 5 2 11 9
2 10 5 3 12 5 2 2 5 4 9 9
2 11 6 1 20 10 2 7 6 2 10 8
2 14 6 3 12 4 2 6 6 4 10 7
2 5 7 1 8 8 2 13 7 2 6 0
2 12 7 3 9 2 2 16 7 4 9 0
2 9 8 1 14 9 2 1 8 2 13 12
2 8 8 3 12 3 2 4 8 4 14 7
3 7 1 5 7 7 3 13 1 6 7 0
3 8 1 7 13 3 3 14 1 8 9 0
3 4 2 5 15 11 3 10 2 6 9 7
3 3 2 7 15 11 3 9 2 8 13 5
3 6 3 5 16 9 3 1 3 6 8 8
3 15 3 7 7 0 3 12 3 8 12 8
3 11 4 5 8 1 3 16 4 6 15 1
3 5 4 7 12 7 3 2 4 8 16 12
4 9 5 5 15 8 4 4 5 6 10 6
4 12 5 7 13 5 4 1 5 8 15 9
4 15 6 5 17 6 4 6 6 6 8 2
4 14 6 7 12 5 4 7 6 8 15 8
4 13 7 5 13 2 4 8 7 6 13 9
4 3 7 7 9 9 4 10 7 8 6 6
4 2 8 5 12 8 4 11 8 6 9 7
4 5 8 7 11 10 4 16 8 8 15 7
;

    如果侵扰是独立于实验单元的，并且一个单元内的所有植物都具有相同的侵扰倾向，则为二项随机变量。

    GLM分析

    让我们先考虑一个独立二项式计数的标准广义线性模型。SAS声明如下：

    proc glimmix data=HessianFly;
class block entry;
model y/n = block entry / solution;
run;

    值得注意的是，glimmix过程支持响应变量的两种语法。此示例使用事件/测试语法。变量y表示n个伯努利试验中的成功（事件）数。
    当使用事件/试验语法时，glimmix过程会自动选择二项分布作为响应分布。一旦确定了分布，该过程将为模型选择链接函数。二项数据的默认链接是logit链接。上述声明等同于以下声明。

    proc glimmix data=HessianFly ;
class block entry;
model y/n = block entry/solution dist=binomial link=logit ;
run;

    “模型信息”表描述了用于拟合统计模型的模型和方法。

    glmm-随机block效应

    “拟合统计”表（皮尔逊比率=2.37）中指出的过度分散有几个可能的原因。数据可能不遵循二项式分布，一个或多个重要影响可能没有在模型中说明，或数据是正相关的。
    如果省略了重要的固定效果，那么您可能需要考虑将它们添加到模型中。由于这是一个设计的实验，所以除了代表处理和误差控制设计结构的block效应和entry效应外，不期望进一步的效应是合理的。过度分散的原因必须在其他地方。如果过度分散源于观测值之间的相关性，则应适当调整模型。相关性可以有多个来源。首先，实验单元内的植物可能没有独立反应。
    如果某一特定植物的侵扰概率被同一单元内相邻植物的侵扰所改变，则侵扰计数不是二项式的，应使用不同的概率模型。第二个可能的关联源是缺乏独立的实验单元。即使处理是随机分配给各个单位，它们也可能没有独立的反应。例如，共享空间土壤效应可能是潜在因素。以下分析考虑了这些空间效应。
    首先，假设环境影响以block尺度运行。通过使block效应随机，边缘响应将相互关联，因为block内的观测具有相同的随机效应。从不同block的观察结果将保持不相关，本着block间独立随机的精神。
    下一组SAS语句适合具有随机块效应的广义线性混合模型（GLMM）：

    proc glimmix data=HessianFly;
class block entry;
model y/n = entry / solution;
random block;
run;

    在存在随机效应和条件二项式分布的情况下，proc glimmix不使用最大概率进行估计。相反，glimmix程序采用一种限制性（残差）伪似然算法。
    “维度”表与以前的模型不同。“维度”表显示有一个G边参数，即随机区组效应的方差。
    请注意，虽然块效应有四个级别，但只估计一个方差分量。但是，Z矩阵有四列，对应于块效果的四个级别。因为在RANDOM语句中没有使用subject=选项，所以glimmix过程将这些数据视为来自一个具有64个观察值的单个subject。
    “优化信息”表表明，采用拟牛顿法求解优化问题。这是GLMM模型的默认方法。
    广义卡方统计量测量最终模型中的残差平方和，其自由度比率是对平均模型观测的可变性的度量。过分散参数（2.25）仍大于1。
    下表中的随机块效应方差很小。随机块模型不能提供适当的弥散调整。
    由于块方差分量很小，与标准GLM相比，III型检验对中的品种效应影响很小。
    在实验设计中，研究人员对每一个数据都有行-列（经度和纬度）值，以解释空间尺度的微观位置变化。

    具有平滑空间趋势的分析

    如果环境效应在小于块大小的空间尺度上运行，则随机块模型不能提供适当的调整。从实验区域的粗略布局来看，单凭随机块效应并不足以解释数据的过度分散，这并不奇怪。我们需要在proc glimmix中添加一个可乘的分散成分（random _residual_;）。
    这样的过度分散分量不影响参数估计，只影响它们的标准误差。另一方面，真正的随机效应会影响参数估计及其标准误差。

    Ods graphics on/noborder ;
proc glimmix data=HessianFly PLOTS=All;
class entry ;
model y/n = entry / solution ddfm=contain;
random _residual_ / subject=intercept type=sp(exp)(lng lat);
lsmeans entry / plots=mean(sliceby=entry join);
run;

    RANDOM语句中的关键字_residual_指示glimmix过程对R矩阵进行建模。这里，R被建模为一个指数协方差结构矩阵。subject=intercept选项意味着所有观测都被认为是相关的。块效果可以保留在模型中。但是相关的R结构吸收了所有的变化，并且没有留下任何（零方差）来解释块效应。我们从最终模型中去掉了块效果。plots=all选项生成模型诊断图。
    空间过程的门槛，即潜在剩余效应的方差，估计为2.5315。空间过程实际范围的三分之一是0.9052。
    与之前的分析相比，entry效应的F值（3.6）急剧下降。 平滑的空间变化解释了品种之间的一些变化。

    结论

    在这个例子中，考虑了三个模型来分析具有二项式结果的随机区组设计。如果数据是相关的，则标准的广义线性模型通常将指示相对于二项分布的过度分散。在该示例中考虑了两个方法以解决这种过度分散。首先，包含G侧随机效应间接模拟相关性;它是通过在同一块的响应中共享随机效应而引发的。
    其次，R侧空间协方差结构直接模拟协变。
    在广义线性（混合）模型中，这两种建模方法可能导致不同的推论，因为模型具有不同的解释。随机块效应在链接（logit）尺度上建模，空间效果在平均尺度上建模。
    只有在线性混合模型中，两个尺度相同。

    带有ASReml的GLMM

    ASReml广泛用于遗传数据分析设计实验。我们使用以下ASReml代码分析了Hessian Fly数据。块效应再次随机拟合。

    Title: Hessianfly.
#Damaged,Plants,block,Entry,lat,lng
#2,8,1,14,1,1
#1,9,1,16,1,2
y N
block *
entry *
lat *
lng *
yRatio !=y !/N
Hessianfly.csv !SKIP 1
y !bin !TOTAL=N ~ mu entry !r block

    下面给出了主输出文件（HessianFly.ASR）文件的部分输出结果。

    Distribution and link: Binomial; Logit Mu=P=1/(1+exp(-XB))
V=Mu(1-Mu)/N
Warning: The LogL value is unsuitable for comparing GLM models
Notice: 1 singularities detected in design matrix.
1 LogL=-46.8426 S2= 1.0000 48 df Dev/DF= 2.671
2 LogL=-46.8446 S2= 1.0000 48 df Dev/DF= 2.671
3 LogL=-46.8446 S2= 1.0000 48 df Dev/DF= 2.671
4 LogL=-46.8446 S2= 1.0000 48 df Dev/DF= 2.671
5 LogL=-46.8446 S2= 1.0000 48 df Dev/DF= 2.671
Final parameter values 1.0000
Deviance from GLM fit 48 128.23
Variance heterogeneity factor [Deviance/DF] 2.67

    异质性因子[Deviance/DF]给出了一些指示，即离散分布与数据的拟合程度。 大于1的值表明数据过度分散，即数据值比所选分布下预期的可变性更大。

    - - - Results from analysis of y - - -
Source Model terms Gamma Component Comp/SE % C
Variance 64 48 1.00000 1.00000 0.00 0 F
Wald F statistics
Source of Variation NumDF DenDF F-inc P-inc
7 mu 1 48.0 4.64 0.036
4 entry 15 48.0 6.87 <.001
Finished: 26 Jul 2011 13:38:06.968 LogL Converged

    entry的F值（植物品种）很大且显著。
    我们可以通过在ASReml中使用!DISPERSION限定符来调整异质性（过度分散）。 如果分析人员不提供，则从残差估计分散参数。 以下是考虑过度分散的模型。

    y !bin !TOTAL=n !dispersion ~ mu entry !r block

    在调整方差的异质性后，我们看到一个小得多的entryF检验。它仍然很显著。 植物品种的预测不会改变，但标准误差会发生变化。

    ASReml的空间R结构

    ASReml功能强大，可以模拟GLMM模型的R侧，特别适合空间和power模型。 在下面的示例中，我们使用R的二维自回归阶1（AR1 x Ar1）相关结构来说明数据中的异质性作为演示。

    !PART 2
! Generalized Linear Mixed Model with spatial R
y !bin !TOTAL=N ~ mu entry mv
1 2 0
8 row AR1 0.1
8 column AR1 0.1

    主要结果如下

     Deviance from GLM fit              48      141.87
 Variance heterogeneity factor [Deviance/DF]     2.96
          - - - Results from analysis of y - - -
 Model_Term                             Sigma         Sigma   Sigma/SE   % C
 ar1(lat).ar1(lng)               64 effects
 lat                      AR_R    1  0.136060      0.136060       1.11   4 P
 lng                      AR_R    1  0.914424E-01  0.914424E-01   0.80  -2 P
                                    Wald F statistics
     Source of Variation           NumDF     DenDF    F-inc            P-inc
   8 mu                                1       5.9     4.60            0.076
   4 entry                            15      48.0    13.50            <.001

    残差的AR1结构不成功，如相关很小（0.0709）所示。 方差异质性因子仍为2.71。 下面给出了基于AR1×AR1的残差样本变异函数。 行或列方向没有明显的趋势，但曲折表面表明数据的异质性。

    相关文章

      网友评论

        本文标题：Isik2011 GLMM

        本文链接：https://www.haomeiwen.com/subject/iphysqtx.html

        延伸阅读

        深度阅读

        • 您也可以注册成为美文阅读网的作者，发表您的原创作品、分享您的心情！

        栏目导航

        热点阅读

        遗传改良（育种）

        关于我们|服务条款|联系我们|Isik2011 GLMM|投稿指南|网站地图|RSS订阅|排版工具|手机版

        提供经典美文摘抄,优美散文欣赏,现代诗歌精选,短篇小说,心情随笔,表白情书范文,故事会在线阅读欣赏

        Copyright © 2014-2023 Haomeiwen.com All Rights Reserved. 好美文阅读网 版权所有

        备案信息：桂公网安备 45052102000051号 · 桂ICP备13007215号-3

        本站所收录作品、热点评论等信息部分来源互联网，目的只是为了系统归纳学习和传递资讯

        所有作品版权归原创作者所有，与本站立场无关，如不慎侵犯了你的权益，请联系我们告知，我们将做删除处理！