全文链接：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 中。然而说我收集了五个不同遗传谱系中的黄蜂。“遗传学”的随机效应与察看无关。它将与其余两个随机效应_穿插_。因而，这种随机效应将 与其余效应 _穿插_。
视频线性混合效应模型 LMM,Linear Mixed 和 R 语言实现
**，时长 12:13
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 轴示意通过散布建模的分位数。红色实线示意现实散布拟合，红色虚线示意现实散布拟合的置信区间。您想抉择最大的观测值落在虚线之间的散布。在这种状况下，这就是对数正态分布，其中只有一个观测值落在虚线之外。当初，我能够尝试拟合模型。
点击题目查阅往期内容
R 语言用 Rshiny 摸索 lme4 狭义线性混合模型（GLMM）和线性混合模型（LMM）
左右滑动查看更多
01
02
03
04
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.58Anova 函数进行了 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")后果很好，因为两个模型之间的估算值十分类似，然而在第二个模型中，对年的置信区间显著较小，阐明这个预计更好。图中能够证实第二种模型的推论，即基因型和年份是变异的次要因素。
本文摘选 《 基于 R 语言混合效应模型（mixed model）案例钻研 》，点击“ 浏览原文”获取全文残缺材料。
点击题目查阅往期内容
R 语言建设和可视化混合效应模型 mixed effect model
【视频】线性混合效应模型 (LMM,Linear Mixed Models) 和 R 语言实现案例
R 语言 线性混合效应模型实战案例
R 语言用潜类别混合效应模型 (Latent Class Mixed Model ,LCMM) 剖析老年痴呆年龄数据
R 语言贝叶斯狭义线性混合（多层次 / 程度 / 嵌套）模型 GLMM、逻辑回归剖析教育留级影响因素数据 R 语言预计多元标记的潜过程混合效应模型（lcmm）剖析心理测试的认知过程
R 语言因子实验设计 nlme 拟合非线性混合模型剖析有机农业施氮程度
R 语言非线性混合效应 NLME 模型 (固定效应 & 随机效应) 反抗哮喘药物茶碱动力学钻研
R 语言用线性混合效应（多程度 / 档次 / 嵌套）模型剖析腔调高下与礼貌态度的关系
R 语言 LME4 混合效应模型钻研老师的受欢迎水平 R 语言 nlme、nlmer、lme4 用（非）线性混合模型 non-linear mixed model 剖析藻类数据实例
R 语言混合线性模型、多层次模型、回归模型剖析学生均匀问题 GPA 和可视化
R 语言线性混合效应模型（固定效应 & 随机效应）和交互可视化 3 案例
R 语言用 lme4 多层次（混合效应）狭义线性模型（GLM），逻辑回归剖析教育留级考察数据 R 语言 线性混合效应模型实战案例
R 语言混合效应逻辑回归（mixed effects logistic）模型剖析肺癌数据
R 语言如何用潜类别混合效应模型（LCMM）剖析抑郁症状
R 语言基于 copula 的贝叶斯分层混合模型的诊断准确性钻研
R 语言建设和可视化混合效应模型 mixed effect model
R 语言 LME4 混合效应模型钻研老师的受欢迎水平
R 语言 线性混合效应模型实战案例
R 语言用 Rshiny 摸索 lme4 狭义线性混合模型（GLMM）和线性混合模型（LMM）
R 语言基于 copula 的贝叶斯分层混合模型的诊断准确性钻研
R 语言如何解决线性混合模型中畸形拟合 (Singular fit) 的问题
基于 R 语言的 lmer 混合线性回归模型
R 语言用 WinBUGS 软件对学术能力测验建设档次（分层）贝叶斯模型
R 语言分层线性模型案例
R 语言用 WinBUGS 软件对学术能力测验（SAT）建设分层模型
应用 SAS，Stata，HLM，R，SPSS 和 Mplus 的分层线性模型 HLM
R 语言用 WinBUGS 软件对学术能力测验建设档次（分层）贝叶斯模型
SPSS 中的多层（等级）线性模型 Multilevel linear models 钻研整容手术数据
用 SPSS 预计 HLM 多层（档次）线性模型模型