原文链接： http://tecdat.cn/?p=2596
1.混合模型是否适合您的需求？
混合模型在很多方面与线性模型相似。它估计一个或多个解释变量对因变量的影响。混合模型的输出将为解释值列表，它们的效果大小的估计值和置信区间，每种效果的p值以及至少一种模型拟合程度的度量。当您有一个变量将数据样本描述为可以收集的数据的子集时，应该使用混合模型而不是简单的线性模型。
让我们看一下正在研究的黄蜂亲属识别数据。
str(data) ## 'data.frame': 84 obs. of 6 variables: ## $ Test.ID : int 1 2 3 4 5 6 7 8 9 10 ... ## $ Observer : Factor w/ 4 levels "Charles","Michelle",..: 1 4 2 4 1 3 2 2 1 2 ... ## $ Relation : Factor w/ 2 levels "Same","Stranger": 1 1 1 1 1 1 1 1 1 1 ... ## $ Aggression: int 4 1 15 2 1 0 2 0 3 10 ... ## $ Tolerance : int 4 34 14 31 4 13 7 6 13 15 ... ## $ Season : Factor w/ 2 levels "Early","Late": 1 1 1 1 1 1 1 1 1 1 ...
我感兴趣的因变量是攻击性和宽容度。侵略性是指六十分钟内的攻击行为次数。宽容是指六十分钟内的宽容行为数量。我对关系（无论黄蜂来自相同还是不同的菌落）和季节（菌落周期的早期或晚期）对这些因变量的影响感兴趣。这些影响是“固定的”，因为无论我在何处，如何采样或采样了多少只黄蜂，我在相同变量中仍将具有相同的水平：相同的菌落与不同的菌落，以及早季与晚季。
但是，还有两个其他变量在样本之间不会保持固定。如果我在不同的年份进行采样，那么观察者的水平会有所不同。样品之间的测试ID也会有所不同，因为我总是可以重新安排哪些黄蜂参与每个实验试验。每个试验都是我当时收集的黄蜂的唯一子样本。如果我能够单独测试黄蜂，并且如果所有观察都对所有互动进行了评分，那么我将不会有任何随机效应。但是，相反，我的数据本来就是“块状”的，随机效应描述了这种块状性。
在继续之前，您还需要考虑随机效果的结构。您的随机效果是嵌套还是交叉？在我的研究中，随机效应是 _嵌套的_，因为每个观察者记录了一定数量的试，并且没有两个观察者记录了相同的试验，因此Test.ID嵌套在Observer中。但是说我收集了五个不同遗传谱系中的黄蜂。“遗传学”的随机效应与观察无关。它将与其他两个随机效应_交叉_。因此，这种随机效应将 与其他效应 _交叉_。
2.哪种概率分布最适合您的数据？
假设您已决定要运行混合模型。接下来要做的是找到最适合您数据的概率分布。有很多测试方法。请注意，负二项式和伽马分布只能处理正数，而泊松分布只能处理正整数。二项分布和泊松分布与其他分布不同，因为它们是离散的而不是连续的，这意味着它们可以量化不同的，可数的事件或这些事件的概率。现在让我们为我的Aggression变量找到一个合适的分布。
require(car) ## 正在加载所需的包: car require(MASS) # 必须为非零的分布 qqp(Aggression, "norm")
# lnorm 表示对数正态 qqp(Aggression, "lnorm")
# qqp需要估计负二项式，泊松和伽玛分布的参数。您可以使用fitdistr函数生成估算值。保存输出并提取每个参数的估计值，如下所示。 fitdistr(rAggression, "Negative Binomial")
qqp(Aggressio, "pois", estimate)
fitdistr(Aggression.t, "gamma")
查看我使用qqp生成的图。y轴表示观测值，x轴表示通过分布建模的分位数。红色实线表示理想分布拟合，红色虚线表示理想分布拟合的置信区间。您想选择最大的观测值落在虚线之间的分布。在这种情况下，这就是对数正态分布，其中只有一个观测值落在虚线之外。现在，我可以尝试拟合模型。
3.如何将混合模型拟合到您的数据
3a.如果您的数据是正态分布的
首先，请注意：如果您的数据最适合对数正态分布， 请不要对其进行_变换_。 由于变换使模型结果的解释更加困难。
如果数据呈正态分布，则可以使用线性混合模型（LMM）。该函数的第一个参数是一个公式，形式为y〜x1 + x2 ...等，其中y是因变量，而x1，x2等是解释变量。交叉随机效应的形式为（1 | r1）+（1 | r2）...，而嵌套随机效应的形式为（1 | r1 / r2）。
在这里，您可以指定混合模型将使用最大似然还是受限最大似然来估计参数。如果您的随机效应是嵌套的，或者只有一个随机效应，并且您的数据是平衡的（即，每个因子组中的样本量相似），则将REML设置为FALSE，因为您可以使用最大似然率。如果交叉了随机效果，请不要设置REML参数，因为无论如何它默认为TRUE。
为了避免这一切看起来太抽象，让我们尝试一些数据。我们将有关八哥歌曲研究的一些数据。在这项研究中，我们对雄性和雌性八哥歌曲之间的差异以及社会地位，不同的鸟类的歌唱是否不同感兴趣。我们的随机效应是社会群体。歌曲的平均音高符合正态概率分布。
str(starlings) ## 'data.frame': 28 obs. of 5 variables: ## $ Individual : Factor w/ 28 levels "B-40917","B-41205",..: 4 5 6 15 3 16 8 13 20 14 ... ## $ Sex : Factor w/ 2 levels "F","M": 2 2 2 2 2 1 1 1 1 2 ... ## $ Group : Factor w/ 5 levels "DRT1","MRC1",..: 2 5 5 4 4 4 4 4 4 4 ... ## $ Social.Rank: Factor w/ 2 levels "Breeder","Helper": 2 1 1 1 2 2 2 2 1 2 ... ## $ Mean.Pitch : num 2911 2978 3313 3268 3312 ... summary(lmm) ## Linear mixed model fit by maximum likelihood \['lmerMod'\] ## Formula: Mean.Pitch ~ Sex + Social.Rank + (1 | Group) ## Data: starlings ## ## AIC BIC logLik deviance df.resid ## 389.3 396.0 -189.7 379.3 23 ## ## Scaled residuals: ## Min 1Q Median 3Q Max ## -2.0559 -0.6272 0.0402 0.5801 2.0110 ## ## Random effects: ## Groups Name Variance Std.Dev. ## Group (Intercept) 0 0 ## Residual 44804 212 ## Number of obs: 28, groups: Group, 5 ## ## Fixed effects: ## Estimate Std. Error t value ## (Intercept) 3099.0 82.2 37.7 ## SexM 51.7 81.3 0.6 ## Social.RankHelper -45.0 82.4 -0.5 ## ## Correlation of Fixed Effects: ## (Intr) SexM ## SexM -0.630 ## Scl.RnkHlpr -0.668 0.106
让我们看看结果。首先，我们获得一些模型拟合的度量，包括AIC，BIC，对数似然度和偏差。然后我们得到由随机效应解释的方差估计。这个数字很重要，因为如果它与零没有区别，那么您的随机效应可能并不重要，您可以继续进行常规的线性模型建立。接下来，我们将对固定效应进行估算，带有标准误差。这些信息可能足以满足您的需求。一些期刊将这些模型的结果报告为带有置信区间的效应大小。当然，当我查看固定效应估算值时，我已经可以看出，性别和社会地位之间的平均音高没有差异。但是有些期刊希望您报告p值。
如果您想要一些p值，则需要使用Anova函数。
## Analysis of Deviance Table (Type II Wald chisquare tests) ## ## Response: Mean.Pitch ## Chisq Df Pr(>Chisq) ## Sex 0.4 1 0.52 ## Social.Rank 0.3 1 0.58
Anova函数进行了Wald检验，该检验告诉我们我们对性别和社会地位对音高的影响的估计p值。
拟合线性混合模型时，可能会遇到一种复杂情况。R可能会有“无法收敛”错误，通常将其表述为“没有收敛就达到了迭代限制”。这意味着您的模型有太多因素，样本量不够大，无法拟合。然后，您应该做的是从模型中删除固定效果和随机效果，然后进行比较以找出最合适的效果。一次删除固定效果和随机效果。保持固定效果不变，并一次删除一个随机效果，然后找出最合适的效果。然后保持随机效果不变，并一次删除固定效果。在这里，我只有一种随机效果，
anova(noranklmm, nosexlmm, nofixedlmm) ## Data: starlings ## Models: ## nofixedlmm: Mean.Pitch ~ 1 + (1 | Group) ## noranklmm: Mean.Pitch ~ Sex + (1 | Group) ## nosexlmm: Mean.Pitch ~ Social.Rank + (1 | Group) ## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq) ## nofixedlmm 3 386 390 -190 380 ## noranklmm 4 388 393 -190 380 0.48 1 0.49 ## nosexlmm 4 388 393 -190 380 0.00 0 1.00
请注意，该方差分析函数与我们用来评估模型中固定效果的重要性的方差分析函数不同。方差分析函数用于比较模型。p值表明模型之间没有明显的重要差异。我们还可以比较AIC值，请注意，具有最低AIC值的模型是完全没有固定影响的模型，这符合我们的理解，即性别和社会地位对歌曲的音调没有影响。无论采用哪种方法，请务必在稿件中报告用于选择最佳模型的p值或AIC值。
基于R语言混合效应模型（mixed model）案例研究-2
