基于R语言混合效应模型（mixed model）案例研究

原文链接： 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值。

3b.如果您的数据不是正态分布的

您会看到，用于估计模型中影响大小的REML和最大似然法做出了不适用于数据的正态假设，因此您必须使用其他方法进行参数估计。问题在于，存在许多替代的估算方法，每种估算方法都使用不同的R包运行，并且很难确定哪种方法合适。

首先，我们需要测试是否可以使用惩罚拟似然（PQL）。PQL是一种灵活的技术，可以处理非正常数据，不平衡设计和交叉随机效应。但是，如果您的因变量符合离散计数分布（例如泊松或二项式）且均值小于5，或者您的因变量为二元变量，则会产生偏差估计。

Aggression变量适合对数正态分布，该分布不是离散分布。这意味着我们可以继续使用PQL方法。但是在继续之前，让我们回到转变为正态的问题。

将分布设置为对数正态，我们将族设置为高斯，并将链接设置为log。


##     lmList
summary(PQL)
## Linear mixed-effects model fit by maximum likelihood
##  Data: recog 
##   AIC BIC logLik
##    NA  NA     NA
## 
## Random effects:
##  Formula: ~1 | Observer
##         (Intercept)
## StdDev:      0.3312
## 
##  Formula: ~1 | Test.ID %in% Observer
##         (Intercept) Residual
## StdDev:      0.5295    7.128
## 
## Variance function:
##  Structure: fixed weights
##  Formula: ~invwt 
## Fixed effects: Aggression.t ~ Relation + Season 
##                   Value Std.Error DF t-value p-value
## (Intercept)       1.033    0.5233 55   1.974  0.0535
## RelationStranger  1.210    0.4674 55   2.589  0.0123
## SeasonLate       -1.333    0.5983 23  -2.228  0.0359
##  Correlation: 
##                  (Intr) RltnSt
## RelationStranger -0.855       
## SeasonLate       -0.123  0.000
## 
## Standardized Within-Group Residuals:
##      Min       Q1      Med       Q3      Max 
## -4.86916 -0.29958 -0.08012  0.14280  5.93336 
## 
## Number of Observations: 84
## Number of Groups: 
##              Observer Test.ID %in% Observer 
##                     4                    28

该模型表明季节对侵略性有影响，也就是说，在菌落周期后期收集的黄蜂比早期收集的黄蜂侵略性小。这也表明黄蜂之间的关系有影响。他们相对于巢友更可能对陌生人有侵略性。我将这些统计数据与估计值，标准误，t值和p值一起报告。

那么，如果您的因变量的平均值小于5，或者您有一个二元因变量，而您不能使用PQL，该怎么办？在这里，您可以使用两种选择：拉普拉斯（Laplace）近似和马尔可夫链蒙特卡罗算法（MCMC）。拉普拉斯近似值最多可以处理3个随机效果。除此之外，您还必须使用MCMC。

让我们从一个可以使用拉普拉斯逼近的例子开始。我们将使用学生在学校的学习情况的数据。出于本示例的目的，我将仅将数据子集化为几个感兴趣的变量，并将“ repeatgr”变量简化为二元因变量。


str(bdf)
## 'data.frame':    2287 obs. of  4 variables:
##  $ schoolNR: Factor w/ 131 levels "1","2","10","12",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ Minority: Factor w/ 2 levels "N","Y": 1 2 1 1 1 2 2 1 2 2 ...
##  $ ses     : num  23 10 15 23 10 10 23 10 13 15 ...
##  $ repeatgr: Factor w/ 3 levels "0","1","2": 1 1 1 1 1 1 1 1 2 1 ...

假设我们要找出是否属于少数民族和社会经济地位会影响学生复读成绩的可能性。我们的因变量是“ repeatgr”，指示学生是否重复了成绩。少数族身份是二元Y / N类别，社会经济地位由“ ses”表示，其数字范围为10至50，其中50是最富有的。我们的随机因素是“ schoolNR”，它代表从中采样学生的学校。因为因变量是二元的，所以我们需要具有二项式分布的广义线性混合模型，并且由于我们的随机效应少于五个，因此可以使用Laplace近似。

严格来说，Laplace近似是一种称为Gauss-Hermite正交（GHQ）的参数估算方法的特例，需要进行一次迭代。GHQ由于重复迭代而比Laplace更为准确，但在第一次迭代后变得不那么灵活，因此您只能将它用于一个随机效果。我们唯一的随机效应是“ schoolNR”。


summary(GHQ)
## Generalized linear mixed model fit by maximum likelihood (Adaptive
##   Gauss-Hermite Quadrature, nAGQ = 25) [glmerMod]
##  Family: binomial ( logit )
## Formula: repeatgr ~ Minority + ses + ses * Minority + (1 | schoolNR)
##    Data: bdf
## 
##      AIC      BIC   logLik deviance df.resid 
##   1672.8   1701.4   -831.4   1662.8     2282 
## 
## Scaled residuals: 
##    Min     1Q Median     3Q    Max 
## -0.964 -0.407 -0.314 -0.221  5.962 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  schoolNR (Intercept) 0.264    0.514   
## Number of obs: 2287, groups: schoolNR, 131
## 
## Fixed effects:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   -0.45194    0.20763   -2.18     0.03 *  
## MinorityY      0.47957    0.47345    1.01     0.31    
## ses           -0.06205    0.00798   -7.78  7.5e-15 ***
## MinorityY:ses  0.01196    0.02289    0.52     0.60    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) MnrtyY ses   
## MinorityY   -0.400              
## ses         -0.905  0.369       
## MinortyY:ss  0.299 -0.866 -0.321
Laplace <- glmer(repeatgr ~ Minority + ses + ses * Minority + (1 | schoolNR),
    data = bdf, family = binomial(link = "logit"))  # Contrast to the Laplace approximation
## Warning: Model failed to converge with max|grad| = 0.0458484 (tol = 0.001)
summary(Laplace)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial ( logit )
## Formula: repeatgr ~ Minority + ses + ses * Minority + (1 | schoolNR)
##    Data: bdf
## 
##      AIC      BIC   logLik deviance df.resid 
##   1672.8   1701.5   -831.4   1662.8     2282 
## 
## Scaled residuals: 
##    Min     1Q Median     3Q    Max 
## -0.960 -0.407 -0.316 -0.222  5.950 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  schoolNR (Intercept) 0.258    0.508   
## Number of obs: 2287, groups: schoolNR, 131
## 
## Fixed effects:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   -0.45456    0.20611   -2.21    0.027 *  
## MinorityY      0.48005    0.47121    1.02    0.308    
## ses           -0.06191    0.00791   -7.83  4.9e-15 ***
## MinorityY:ses  0.01194    0.02274    0.53    0.600    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) MnrtyY ses   
## MinorityY   -0.400              
## ses         -0.906  0.369       
## MinortyY:ss  0.299 -0.866 -0.321

您可以看到在这种情况下，Laplace和GHQ之间没有重要区别。两者都表明，社会经济课对学生复读成绩的可能性有非常显着的影响，尽管即使采用logit变换，我们也可以看到影响量很小。但是，使用此方法时，还需要考虑其他一些因素。当用于过度分散的数据时，即合并的残差远大于残差自由度时，它变得不准确。使用此方法时，应检查模型以确保数据不会过度分散。


    ## n×n方差-协方差矩阵中方差参数的数量


    vpars <- function(m) {
        nrow(m) * (nrow(m) + 1)/2
    }
    # 接下来计算剩余自由度

    rdf <- nrow(model.frame(model)) - model.df
    # 提取皮尔逊残差


    prat <- Pearson.chisq/rdf
    # 生成一个p值。如果小于0.05，则数据过于分散。

这些学校数据并不分散。如果是，可以将过度分散建模为随机效应，每个观察值具有一个随机效应水平。在这种情况下，我将使用学生ID作为随机效果变量。

我们继续讨论泊松数据的均值太小或因变量是分类的情况，并且我们有五个或五个以上随机效应。考虑有关大麦农民的这些数据。想象一下，要使大麦收成产生利润，大麦收成的收入必须大于140。


farmers <- farmers[c(1:5, 8)]
str(farmers)
## 'data.frame':    102 obs. of  6 variables:
##  $ year    : int  1901 1901 1901 1901 1902 1902 1902 1902 1902 1902 ...
##  $ farmer  : Factor w/ 18 levels "Allardyce","Dooley",..: 10 5 3 18 10 5 18 17 4 12 ...
##  $ place   : Factor w/ 16 levels "Arnestown","Bagnalstown",..: 3 16 14 11 3 16 11 4 8 6 ...
##  $ district: Factor w/ 4 levels "CentralPlain",..: 2 2 1 1 2 2 1 1 4 4 ...
##  $ gen     : Factor w/ 2 levels "Archer","Goldthorpe": 1 1 1 1 1 1 1 1 1 1 ...
##  $ profit  : num  1 1 1 1 1 1 1 1 1 1 ...

在这种情况下，假设我们根本没有解释变量。我们对哪些随机效应会影响利润感兴趣。毕竟，随机效应是改变因变量方差的因素。为此，我们将处理随机效应，而且还为随机效应提供置信区间。

所有模型都对数据中方差的分布进行假设，但是在贝叶斯方法中，这些假设是明确的，因此我们需要指定这些假设的分布。在贝叶斯统计中，我们称这些先验。


##  Iterations = 3001:12991
##  Thinning interval  = 10
##  Sample size  = 1000 
## 
##  DIC: 76.61 
## 
##  G-structure:  ~year
## 
##      post.mean l-95% CI u-95% CI eff.samp
## year      18.2    0.266     68.6     18.1
## 
##                ~farmer
## 
##        post.mean l-95% CI u-95% CI eff.samp
## farmer      1.93   0.0789     6.68     99.7
## 
##                ~place
## 
##       post.mean l-95% CI u-95% CI eff.samp
## place      1.54   0.0711     4.84      173
## 
##                ~gen
## 
##     post.mean l-95% CI u-95% CI eff.samp
## gen      12.1   0.0874     28.6     1000
## 
##                ~district
## 
##          post.mean l-95% CI u-95% CI eff.samp
## district      1.78   0.0639     5.63     1000
## 
##  R-structure:  ~units
## 
##       post.mean l-95% CI u-95% CI eff.samp
## units         1        1        1        0
## 
##  Location effects: profit ~ 1 
## 
##             post.mean l-95% CI u-95% CI eff.samp pMCMC
## (Intercept)      3.47    -1.16    10.75      220  0.14

如果我们看一下均值和置信区间，我们可以看到，真正影响利润可变性的唯一随机效应是年份和基因型。也就是说，大麦收成获利的可能性每年之间以及基因型之间变化很大，但是在农民或地方之间变化不大。

结束：了解您的数据

除非您熟悉这些分析，否则您将无法真正知道哪种分析对您的数据适用，而熟悉它们的最佳方法是绘制它们。通常，我的第一步是绘制我感兴趣的变量的密度图。


ggplot(recog, aes(x = Aggression)+ geom_density() +

在这里，我按季节和关系（两个固定效应）划分了数据。我们可以立即看到数据集包含一个极端正的异常值；大多数观测值都介于0到20之间。我们还可以看到，后期观测值的很大一部分等于零。

绘图对于评估模型拟合也很重要。您可以通过各种方式绘制拟合值来判断适合的模型最能描述数据。一个简单的应用是绘制模型的残差。如果您将模型想象为通过数据散点图的最佳拟合线，则残差为散点图中各点与最佳拟合线之间的距离。如果模型适合，则将残差与拟合值作图，则应该看到随机散布。如果散布不是随机的，则意味着还有其他随机或固定的影响。

让我们尝试绘制拟合八哥的歌曲音高的混合模型的残差。此图所做的是创建一条表示零的水平虚线：与最佳拟合线的零偏差平均值。它还创建了一条实线，代表与最佳拟合线的残差。我希望实线覆盖虚线。

结果很好：与最佳拟合线的偏差趋于零。

让我们尝试一种以图形方式比较MCMC模型的技术。我们将回到大麦农户的示例。

Trace(MCMC, log = TRUE)

这些随机效果看起来不像白噪声。因此，让我们尝试使用更多迭代来重新拟合模型。这需要更多的计算量，但会产生更准确的结果。

Trace(MCMC2, log = TRUE)

现在，看起来都更接近直线周围的白噪声，这表明模型更好。现在，让我们比较两个模型之间随机效应方差的置信区间。

# 将两个模型的估计值和置信区间放在一起


rbind (covariances, Gcovariances)
# 创建一个数据框架，其中包含模型和随机效应的因素

data.frame(coint, model = factor(el1", "model2"), eac),
    random.effect = dimnames(i)
## Warning: some row.names duplicated: 6,7,8,9,10 --> row.names NOT used
# 将估计和置信区间绘制为按模型分组的箱形图


ggplot(conf.int+ geom_crossbar(aes(y.95..CI,
    y.95..CI= model= "dodge")