假设我们测量了每个被试(Subject)的身高(Height)和体重(Weight),每一年测量一次(Time)。现在我们想研究的是,身高的变化率和体重的变化率是否存在关联?这里我们先假设使用线性混合模型(LME model)的假设是满足的。我在文献上看到了不同的做法,一种是先分别对身高和体重拟合LME,提取每个被试的随机斜率(random slope,我理解就是身高或体重随时间的变化率),然后计算身高和体重随机斜率的皮尔逊相关,大概代码如下:

library(lme4)
## Fit LME for height and weight
height_mod <- lme4::lmer(Height ~ Time +  (1+Time|Subject), data = mydat)
weight_mod <- lme4::lmer(Weight ~ Time +  (1+Time|Subject), data = mydat) 
## 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)                     

我的问题是,这样的做法是合理的吗?虽然我在文献中看到的这种做法,但是文献中统计方法有误的情况也并不少,所以希望大佬能确认一下。

另外,我在文献中也看到另外一种方法,声称是直接在LME中检验的变化率的关联,但是由于描述不是很清楚,所以我不知道应该如何设置模型。

目前看起来合理,个人建议是保守一点,用spearman相关系数,因为随机斜率是否符合正态分布这个不清楚。

至于 “但是文献中统计方法有误的情况也并不少,所以希望大佬能确认一下。”的问题, 这个基本还是经验判断。 一般来说经过同行评议的期刊上的研究paper不一定方法最好最恰当,但基本逻辑还是可以保障的。

相对于你给出办法,以下补充两个思考,都是扩展。

  1. 按照一般地经验性观察,身高和体重是有关系的,所以,身高和体重可以作为一个二维的响应变量,其相关性的估计也由模型给出。多重线性混合效应模型。
  2. 按照一般地经验性观察,身高的变化率和体重的变化率在人的不同年龄段显然是不同的,比如到了 20 多岁以后,很多人不长个子了,30 岁以后体重开始涨了。线性的模型也许需要根据年龄来调整,比如某种类型的生长曲线,常见的有逻辑斯蒂 S 曲线(高中生物课好像提及这个)。多重广义或非线性混合效应模型。

可以做一下统计模拟来理解模型在不同情况下的表现。

    6 天 后

    方法上感觉没问题。不过从直觉上,更推荐用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

    nan.xiao 感谢回复,请问您的意思是,构造一些模拟数据,比如在不同真实参数的情况下,看看估计的结果如何?是这个思路吗?

    感谢楼上各位老师的热心回复!我这里延伸一个相关的问题,因为是相关联的问题,所以不再新开一贴。

    假设除了身高和体重是每年进行测量的,还有一个变量(比如,被试父母的平均身高ParentHeight)是不随时间变化的(或者只在基线测量的),想检验的问题是,被试父母的平均身高和被试身高的变化率是否存在关联(比如,父母平均身高越高,被试身高变化率越大)。当然,这是完全假想的例子,举身高的例子只是为了避免引入专业知识。

    这里似乎有两种做法:第一种是算出每个被试身高的变化率,然后和父母的平均身高进行皮尔逊相关(或者其他相关);第二种是在LME模型中检验时间和父母的平均身高的交互效应。这两种做法得到的结果的解读是很不一样的,如果是第一种,那么我们可以说父母的平均身高与子女的身高变化率是否显著相关,比如相关系数为0.6;如果是第二种,按照交互作用的常规解读,我们可以说子女的身高变化率跟父母的平均身高有关,如果父母平均身高比较高,那么子女身高变化率可能比较大。不知道各位老师如何看待这个问题?

    library(lme4)
    
    ## Method 1
    ## Fit LME for height
    height_mod <- lme4::lmer(Height ~ Time +  (1+Time|Subject), data = mydat)
    ## Extract random slope
    height_rs <- coef(height_mod)$Subject[['Time']]
    ## Correlation between random slope and ParentHeight
    cor.test(height_rs, ParentHeight) 
    
    ## Method 2
    ## Examine  the interaction between Time and ParentHeight
    height_mod <- lme4::lmer(Height ~ Time + ParentHeight + Time:ParentHeight + (1+Time|Subject), data = mydat)
    ## Using likelihood test to get the p-value
    height_mod_null <- lme4::lmer(Height ~ Time + ParentHeight + (1+Time|Subject), data = mydat)
    anova(height_mod_null, height_mod)
      12 天 后

      huyang 顶一下自己的帖子,如果不符合论坛规范,请版主删除。

        14 天 后

        您发的问题答案,显然就是两个因变量了!

        我自己读文章少,也确实没有见过。

        至于两变量的变化率比较,个人觉得就是变量的比较,
        如果真抠数学上的字眼,那就是导数比较了,这个我还真不会。。

        实际上,您把一个当自变量,另外一个当因变量,其他因素作为协变量,就可以慢慢分析了。