基于R语言混合效应模型（mixed model）案例研究-1
https://developer.aliyun.com/article/1489317
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")
结果很好，因为两个模型之间的估算值非常相似，但是在第二个模型中，对年的置信区间明显较小，说明这个估计更好。图中可以证明第二种模型的推论，即基因型和年份是变异的主要因素。
