方法上感觉没问题。不过从直觉上,更推荐用lavaan
包,基于sem的逻辑做潜增长模型(能用模型一次估计完成的,就不估计第二次)。
这是生成模拟数据测试的结果。设置的增长率相关性为0.2,分别使用了基于lme4和lavaan的两种方法。
n = 500 # 500个人
timepoints = 4 # 4个观测时点
Time = rep(0:3, times = n) # 时点标记
Subject = rep(1:n, each=4) # 个体标记
Height_intercept_b0 <- Weight_intercept_b0 <- .5 # 截距
Height_slope_b0 <- Weight_slope_b0 <- .25 # 变化率
# 截距和变化率的随机效应
# 这里假设只有Height_slope和Weight_slope存在相关性,cor=0.2
randomEffectsCorr = matrix(c(
1, 0, 0, 0,
0, 1, 0, .2,
0, 0, 1, 0,
0, .2, 0, 1
), ncol=4)
randomEffects = MASS::mvrnorm(n, mu=c(0,0,0,0), Sigma = randomEffectsCorr, empirical=T) |>
data.frame()
colnames(randomEffects) = c('Height_int', 'Height_slope', 'Weight_int', 'Weight_slope')
randomEffects <- randomEffects |>
as_tibble() |>
mutate(Subject = row_number())
# 将截距和变化率合成dataframe
library(tidyverse)
d <- tibble(
Subject = Subject,
Time = Time,
) |>
left_join(randomEffects, by = "Subject") |>
# 生成Height和Weight的模拟数据
mutate(Height = Height_intercept_b0 + Height_int + (Height_slope_b0 + Height_slope) * Time + rnorm(n*timepoints, 0, .5),
Weight = Weight_intercept_b0 + Weight_int + (Weight_slope_b0 + Weight_slope) * Time + rnorm(n*timepoints, 0, .5))
#-------------
# lme4的方法
#-------------
library(lme4)
## Fit LME for height and weight
height_mod <- lme4::lmer(Height ~ Time + (1+Time|Subject), data = d)
weight_mod <- lme4::lmer(Weight ~ Time + (1+Time|Subject), data = d)
## Extract random slopes
height_rs <- coef(height_mod)$Subject[['Time']]
weight_rs <- coef(weight_mod)$Subject[['Time']]
## Correlation between random slopes
cor.test(height_rs, weight_rs) # 0.2034596
#-------------
# lavaan的方法
#-------------
library(lavaan)
d_wide <- d |>
select(Subject, Time, Height, Weight) |>
pivot_wider(names_from = Time, values_from = c(Height, Weight))
model = "
# intercept and slope with fixed coefficients
Height_i =~ 1*Height_0 + 1*Height_1 + 1*Height_2 + 1*Height_3
Height_s =~ 0*Height_0 + 1*Height_1 + 2*Height_2 + 3*Height_3
Weight_i =~ 1*Weight_0 + 1*Weight_1 + 1*Weight_2 + 1*Weight_3
Weight_s =~ 0*Weight_0 + 1*Weight_1 + 2*Weight_2 + 3*Weight_3
"
growthCurveModel = growth(model, data=d_wide)
summary(growthCurveModel) # 0.218