Mixed effects models 回顾（二）

最新推荐文章于 2024-07-10 23:17:21 发布

芥末与薄煎饼

最新推荐文章于 2024-07-10 23:17:21 发布

阅读量53

点赞数

文章标签：数学建模

本文链接：https://blog.csdn.net/weixin_42850327/article/details/131680843

版权

Mixed effects models 回顾（二）

- - Normalised resuidals
  - GLS-A protocol

A.F. Zuur et al., Mixed Effects Models and Extensions in Ecology with R, 11 Statistics

###写在Chapter 4, 5, 6前

在一般 (general) 线性模型以及一般加性模型的模型验证时，总会碰到一系列问题，诸如异质性、嵌套数据、数据时空相关等。对这些问题的处理，可以通过添加随机部分 (random part) 的形式解决：

请添加图片描述

我们不能简单的将上述数学模型统称为混合效应模型，这是不严谨的。对于异质性问题解决的模型，应被称为广义最小二乘法 (generalised least squares, GLS) 模型，其本质上是一种加权 (weights) 的线性回归。

##Chapter 4 广义最小二乘法(GLS): 处理模型异质性问题

Smith et al.(2005)报道了乌贼产卵(reproductive)及体细胞组织(somatic tissues)的季节变换模式，以此为例，构建睾丸重量(testis weight)与胴背长(dorsal mantle length)及月份关系的模型，线性回归模型为：

$Testisweight_i = intercept + DML_i + Month_i + DML_i:Month_i + ε_i$

$ε_i ~ (0, σ^2)$

在该模型中，index i是观测数(1-768)，ε_i的特征，var(ε_i) = σ²。

M1 <- lm(Testisweight ~ DML*fMONTH, data = Squid)
plot(M1, which = c(1), col = 1, add.smooth = FALSE, caption = "")
plot(Squid$fMONTH, resid(M1), xlab = "Month", ylab = "Residuals")
plot(Squid$DML, resid(M1), xlab = "DML", ylab = "Residuals")

请添加图片描述
Panel A中残差随着Fitted values的增大而增大，违反了同质性假设。另一个明显的模式是，Panel C中，残差随着DML的增大而增大，且不同月份的残差也明显不同。因此，与其假设var(ε_i) = σ²，不如允许var(ε_i) 是解释变量的函数。

用于残差方差的解释变量称为方差协变量 (variance covariate)，找到适合ε_ij的**方差协变量结构 (variance covariate structure)**是关键。常见的方差协变量结构如下表所示：

Names of the Function in R	What does it do?
VarFixed	Fixed variance
VarIdent	Different variances per stratum
VarPower	Power of the variance covariate
VarExp	Exponential of the variance covariate
VarConstPower	Constant plus power of the variance covariate
VarComb	A combination of variance functions

###Fixed Variance Structure

Fix variance structure，该结构假设：

$Testisweight_i = intercept + DML_i + Month_i + DML_i:Month_i + ε_i$

$(ε_i ∼ N(0, σ^2 \times DML_i)$

在该模型中，index i是观测数(1-768), 允许残差与DML呈正比例变化。

library(nlme)
M.lm <- gls(Testisweight ∼ DML * fMONTH, data=Squid)
vf1Fixed <- varFixed(∼DML) #允许残差与DML呈正比例变化。
M.gls1 <- gls(Testisweight ∼ DML * fMONTH, weights = vf1Fixed, data = Squid)
anova(M.lm, M.gls1)

Model df      AIC      BIC    logLik
M.lm       1 25 3752.084 3867.385 -1851.042
M.gls1     2 25 3620.898 3736.199 -1785.449

###The VarIdent Variance Structure

VarIdent Variance Structure，该结构假设：

$Testisweight_{ij} = intercept + DML_{ij} + Month_j + DML_{ij}:Month_j + residuals_{ij}$