R数据分析：潜增长模型LGM的做法和解释，及其与混合模型对比

今天收到了北京大学老师打来的电话，问我如果没有被数据科学方向的导师录取，愿不愿意去读生物统计的博士。

我婉拒了，些许遗憾，但不后悔，原因全是个人选择，读博挺好的，但是我决定换一种环境，去工作了。

从去年11月开始申请，到一系列的纠结，到现在做下决定，确实是释然了很多。读书很自由，但读书并不是适合每一个人的选择，或者说不是适合一个人特定时期的选择，也许工作不顺意又想去读也说不定。

真的是越长大越体会到人生重要决定时的艰难迷茫。怎么选都会遗憾，只有勇敢的走下去吧。

感慨一下哈，今天继续给大家写潜增长模型。

随机效应

要理解随机效应，首先还是要理解嵌套，比如同一个学校的学生，或者同一个人的纵向测量数据，这些数据我们不能想当然地认为它们之间是独立的，但是线性回归的前提假设便是观测的独立性，所以对于嵌套数据一般不能用普通的回归分析法。

对于嵌套数据我们要用混合效应模型，混合效应模型之前的文章写了很多，大家可以自己翻翻哦，为啥叫混合效应模型呢？

They are mixed, because there is a mixture of fixed effects and random effects.

因为这个模型可以同时估计固定效应和随机效应，所以就叫它混合效应模型，而其中的固定效应就和我们普通回归分析的回归系数一样，是我们主要关心的东西。

加上随机效应就是为了让我们能够把嵌套数据的变异分解的更加的清晰：我们可以让每个学校有每个学校自己的特征（随机截距或随即斜率），也可以让纵向测量中的每个人有自己的特征（随机截距或随即斜率），之所以要加随机效应，是我们想让嵌套数据的变异分解地更好从而使得我们的固定效应估计更准确。

混合效应模型的构成

给大家写写混合效应模型构成的实例，比如现在我要研究y和x的关系，但是数据是从学生中收集的，学生来自不同的地区，那么这就是一个嵌套，学生嵌套在不同的地区。我们规定学生用i表示，地区用c表示，拟合的模型如下：

上面这个模型中有一个随机截距b0c，就是说不同的地区学生的起始x允许不一样：

怎么不一样呢？

就是这个随机截距b0c为整体截距b0和不同地区截距uc的和（地区水平的变异可大可小，但我们要考虑它，这就是随机效应的作用），如下式：

我们会假定所有的uc都是服从正态分布的。

此时我们把随机截距模型改写一下就成了下面这个样子：

或者：

看第二个式子你就会明白随机截距只是在常规回归模型上多了一个截距项uc，所以叫做随机截距。

把随机截距放在结构方程模型中理解

刚刚写了随机截距是服从正态分布的，回想我们在结构方程中潜变量也是服从正态分布的，有没有什么联想呢？

是不是可以把随机截距当作潜变量处理呢？

我们来试试，首先模拟数据，我们模拟一个纵向数据，共500个观测，4个时间点：

set.seed(1234)
n = 500
timepoints = 4
time = rep(0:3, times=n)
subject = rep(1:n, each=4)

我们会有固定效应，还有随机效应，对于本例来讲，测量是嵌套在每个观测上的（每个观测测了4次），所以我们的随机效应便是个体水平的截距或斜率，我们让固定效应的截距和斜率分别为 0.5 和 0.25，随机截距和随机斜率的相关为0.2

intercept = .5
slope = .25
randomEffectsCorr = matrix(c(1,.2,.2,1), ncol=2)
randomEffects = MASS::mvrnorm(n, mu=c(0,0), Sigma = randomEffectsCorr, empirical=T) %>%
data.frame()
colnames(randomEffects) = c('Int', 'Slope')

于是此时我们的数据长这样：

我们接着模拟自变量和应变量，我们的思路就是把随机效应加在固定效应上，就是把随机截距和固定截距相加，把随机斜率和固定斜率相加，自变量我们就用时间time就行，应变量的模拟代码如下：

sigma = .5
y1 = (intercept + randomEffects$Int[subject]) + # 随机截距
(slope + randomEffects$Slope[subject])*time + #随机斜率
rnorm(n*timepoints, mean=0, sd=sigma)

d = data.frame(subject, time, y1)

好了，现在这个d便是我们要用的数据集，它长这样：

因为这个数据集是我们亲自模拟出来的，所以我们知道它既有随机斜率又有随机截距，而且我们还知道固定效应和随机效应大小分别是多少，我们现在就来逆操作，看看我们拟合一个混合模型能不能得到我们预想的系数：

library(lme4)
mixedModel = lmer(y1 ~ time + (1 + time|subject), data=d) # 带有随机截距和随机斜率的混合模型
## summary(mixedModel)

运行上面的代码，就可以见证奇迹了：

我们可以看到，我们设定的固定斜率为0.25，图中结果输出为0.261，我们设定的固定截距为0.5，图中结果输出为0.487，因为我们本身还加了随机误差，除过这个随机误差之后可以说是完全吻合。

跑潜增长模型

上面依然是给大家回忆和混合效应模型，希望大家对混合模型已经没问题了，我们接着看潜增长模型：

在潜增长模型中我们会把随机截距和随机斜率当作潜变量处理，相应地我们会固定住因子载荷，看下面这个示意图，我们有随机截距和随机斜率，我们还把载荷都给固定了，注意，我们把不同时间点的截距的载荷都固定为1，斜率载荷的固定一定要反映出我们数据的实际间隔，通常是从0开始。

我们做潜增长模型的时候每个时间点的测量都是以变量纳入的，所以我们的数据得是宽型数据，转换方法如下：

dWide = d %>%
spread(time, y1) %>%
rename_at(vars(-subject), function(x) paste0('y', x))
head(dWide)

我们的数据就成了每个观测分别在不同的时间点时的y值，此时时间点变成了变量而非变量的水平（看不懂这句话的话去看看tidydata）。

数据处理好之后我们就可以来跑增长模型了，用到的是lavaan包中的growth方法：

model = "
# 将随机效应作为潜变量
i =~ 1*y0 + 1*y1 + 1*y2 + 1*y3
s =~ 0*y0 + 1*y1 + 2*y2 + 3*y3
"
growthCurveModel = growth(model, data=dWide)
summary(growthCurveModel)

运行上面的代码便可以得到增长模型的结果，如下图：

大家注意，结果很多地方都是空的，我们需要看的结果就是随机效应的截距项：

结果中有我们的潜变量i和s的截距，这个东西啥意思呢？

回忆一下，我们现在是在跑随机效应，我们是将随机截距和随机斜率当作了潜变量，那么这个i和s的截距就是整体的均值，大家可能又糊涂了，就是说，我们的随机效应是服从正态分布的，它的均值就是固定效应，随机效应会在固定效应周围，而这个潜变量i和s的截距就是均值，所以你就可以把i和s的截距认为是固定效应，希望我写明白了，还不懂的话私信我吧。

还有点怀疑？

哈哈，没事我们来验证下嘛。记住上图的截距分别是0.487和0.261。

我们跑潜增长模型的数据和混合效应的数据完全是一样的，我模拟数据的时候是规定截距和斜率分别是0.5和0.25的，我们i和s的截距是0.487和0.261，考虑我加了随机误差，所以这些值是完全一致的，我们还可以用以下的代码调出来固定效应：

fixef(mixedModel)

看到没，依然是0.487和0.261，和我们的增长模型中的i和s的截距是一模一样的。所以说，大家记住，增长模型是跑的随机效应，所以我们要关心截距输出，截距输出就是整体均值，也就是固定效应。

增长模型的结果输出中还有一部分结果是方差的估计：

这个该怎么看呢？这个不重要，但是大家可以了解一下

潜增长模型给每个时间点都会估计一个方差，但是在混合模型中认为不同时点的方差是一样的，所以你会看到在混合模型中std.Dev为0.48，而在增长模型中每个时点的variance都不一样，平均起来也是0.48.这个就是混合模型和增长模型的区别之一。

为了大家更好地理解上面的叙述，再给大家写个例子：

我们用同样的数据我们做增长模型的时候把方差固定住：

#将增长模型的变量方差固定就和混合模型结果一样了
model = "
# intercept and slope with fixed coefficients
i =~ 1*y0 + 1*y1 + 1*y2 + 1*y3
s =~ 0*y0 + 1*y1 + 2*y2 + 3*y3

y0 ~~ resvar*y0
y1 ~~ resvar*y1
y2 ~~ resvar*y2
y3 ~~ resvar*y3
"

growthCurveModel = growth(model, data=dWide)
summary(growthCurveModel)

结果如下：

还是同样的数据，我们再对比其与混合模型的结果：

可以看出来，增长模型固定方差后就和混合模型的结果是一样的。也就是说大家记住，增长模型的和混合模型的区别就是增长模型认为不同时间变量的变异是不一样的。

小结