使用 lmLst 和间隔函数获取每位患者的截距和斜率（作为重复测量研究中的诊断）

Question

我有一个包含 24 名中风患者的重复测量数据集，我想在其中评估三种不同类型的康复 (Group) 对功能恢复评分 (Barthel_index) 的影响。每周测量每位患者的功能能力 (Time_num)，持续 8 周。

数据如下：

library(dplyr)
library(magrittr)
library(nlme)
library(lmer)

mydata <- 
structure(list(Subject = c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 4L, 
4L, 4L, 4L, 4L, 4L, 4L, 4L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 6L, 
6L, 6L, 6L, 6L, 6L, 6L, 6L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 8L, 
8L, 8L, 8L, 8L, 8L, 8L, 8L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 10L, 
10L, 10L, 10L, 10L, 10L, 10L, 10L, 11L, 11L, 11L, 11L, 11L, 11L, 
11L, 11L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 13L, 13L, 13L, 
13L, 13L, 13L, 13L, 13L, 14L, 14L, 14L, 14L, 14L, 14L, 14L, 14L, 
15L, 15L, 15L, 15L, 15L, 15L, 15L, 15L, 16L, 16L, 16L, 16L, 16L, 
16L, 16L, 16L, 17L, 17L, 17L, 17L, 17L, 17L, 17L, 17L, 18L, 18L, 
18L, 18L, 18L, 18L, 18L, 18L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 
19L, 20L, 20L, 20L, 20L, 20L, 20L, 20L, 20L, 21L, 21L, 21L, 21L, 
21L, 21L, 21L, 21L, 22L, 22L, 22L, 22L, 22L, 22L, 22L, 22L, 23L, 
23L, 23L, 23L, 23L, 23L, 23L, 23L, 24L, 24L, 24L, 24L, 24L, 24L, 
24L, 24L), Group = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L), .Label = c("A", "B", "C"), class = "factor"), 
    Time_num = c(1, 2, 3, 4, 5, 6, 7, 8, 1, 2, 3, 4, 5, 6, 7, 
    8, 1, 2, 3, 4, 5, 6, 7, 8, 1, 2, 3, 4, 5, 6, 7, 8, 1, 2, 
    3, 4, 5, 6, 7, 8, 1, 2, 3, 4, 5, 6, 7, 8, 1, 2, 3, 4, 5, 
    6, 7, 8, 1, 2, 3, 4, 5, 6, 7, 8, 1, 2, 3, 4, 5, 6, 7, 8, 
    1, 2, 3, 4, 5, 6, 7, 8, 1, 2, 3, 4, 5, 6, 7, 8, 1, 2, 3, 
    4, 5, 6, 7, 8, 1, 2, 3, 4, 5, 6, 7, 8, 1, 2, 3, 4, 5, 6, 
    7, 8, 1, 2, 3, 4, 5, 6, 7, 8, 1, 2, 3, 4, 5, 6, 7, 8, 1, 
    2, 3, 4, 5, 6, 7, 8, 1, 2, 3, 4, 5, 6, 7, 8, 1, 2, 3, 4, 
    5, 6, 7, 8, 1, 2, 3, 4, 5, 6, 7, 8, 1, 2, 3, 4, 5, 6, 7, 
    8, 1, 2, 3, 4, 5, 6, 7, 8, 1, 2, 3, 4, 5, 6, 7, 8, 1, 2, 
    3, 4, 5, 6, 7, 8), Barthel_index = c(45L, 45L, 45L, 45L, 
    80L, 80L, 80L, 90L, 20L, 25L, 25L, 25L, 30L, 35L, 30L, 50L, 
    50L, 50L, 55L, 70L, 70L, 75L, 90L, 90L, 25L, 25L, 35L, 40L, 
    60L, 60L, 70L, 80L, 100L, 100L, 100L, 100L, 100L, 100L, 100L, 
    100L, 20L, 20L, 30L, 50L, 50L, 60L, 85L, 95L, 30L, 35L, 35L, 
    40L, 50L, 60L, 75L, 85L, 30L, 35L, 45L, 50L, 55L, 65L, 65L, 
    70L, 40L, 55L, 60L, 70L, 80L, 85L, 90L, 90L, 65L, 65L, 70L, 
    70L, 80L, 80L, 80L, 80L, 30L, 30L, 40L, 45L, 65L, 85L, 85L, 
    85L, 25L, 35L, 35L, 35L, 40L, 45L, 45L, 45L, 45L, 45L, 80L, 
    80L, 80L, 80L, 80L, 80L, 15L, 15L, 10L, 10L, 10L, 20L, 20L, 
    20L, 35L, 35L, 35L, 45L, 45L, 45L, 50L, 50L, 40L, 40L, 40L, 
    55L, 55L, 55L, 60L, 65L, 20L, 20L, 30L, 30L, 30L, 30L, 30L, 
    30L, 35L, 35L, 35L, 40L, 40L, 40L, 40L, 40L, 35L, 35L, 35L, 
    40L, 40L, 40L, 45L, 45L, 45L, 65L, 65L, 65L, 80L, 85L, 95L, 
    100L, 45L, 65L, 70L, 90L, 90L, 95L, 95L, 100L, 25L, 30L, 
    30L, 35L, 40L, 40L, 40L, 40L, 25L, 25L, 30L, 30L, 30L, 30L, 
    35L, 40L, 15L, 35L, 35L, 35L, 40L, 50L, 65L, 65L)), row.names = c(NA, 
-192L), class = c("tbl_df", "tbl", "data.frame"))

head(mydata)

# A tibble: 6 x 4
  Subject Group Time_num Barthel_index
    <int> <fct>    <dbl>         <int>
1       1 A            1            45
2       1 A            2            45
3       1 A            3            45
4       1 A            4            45
5       1 A            5            80
6       1 A            6            80

为了查看每个患者的截距和斜率是否以及如何变化，我想使用 lmList 和 interval 函数绘制截距和斜率。

问题1我不明白为什么在lme4中调用lmList函数()会给我48个警告，而在nlme中调用相同的函数却不会:

lmlist <- 
  lme4::lmList(Barthel_index ~ Time_num | Subject, 
         data=mydata)
> There were 48 warnings (use warnings() to see them)

lmlist <- 
  nlme::lmList(Barthel_index ~ Time_num | Subject, 
         data=mydata)
# Works fine

问题 2 我正在尝试提取每个回归斜率的置信区间，但这给出了某些值的警告和 NaN：

lmlist <- 
  nlme::lmList(Barthel_index ~ Time_num | Subject, 
         data=mydata)
coefs <- coef(lmlist)
names(coefs) <- c("Intercepts", "Slopes")
intervals(lmlist)

> Warning message:
In summary.lm(el) : essentially perfect fit: summary may be unreliable

问题 3 现在我有了带有置信区间的新系数列表，我想绘制它们以查看不同患者之间的截距和斜率是否不同以及有多少不同。我正在尝试实现以下目标：

有什么帮助吗？谢谢。

Answer 1

对于第 3 季度，您可以使用 lattice 包中的 dotplot 函数：

require(lattice)
m0 <- lmer(Reaction ~ Days + (Days | Subject), data = sleepstudy)
dotplot(ranef(m0, condVar = TRUE))

Answer 2

Q1。警告出现在 lme4::lmList 中，因为您使用小标题作为输入：没有来自

的警告

lme4::lmList(Barthel_index ~ Time_num | Subject,
         data=as.data.frame(mydata))

（这是 lme4 中的一个无害的“不当行为”或漏洞 ...）

Q2。如果您查看系数列表，您会发现主题 5 是有问题的。该主题的数据都具有相同的响应值：因此我们无法计算线性回归拟合的置信区间也就不足为奇了...

 mydata[mydata$Subject=="5",]
# A tibble: 8 × 4
  Subject Group Time_num Barthel_index
    <int> <fct>    <dbl>         <int>
1       5 A            1           100
2       5 A            2           100
3       5 A            3           100
4       5 A            4           100
5       5 A            5           100
6       5 A            6           100
7       5 A            7           100
8       5 A            8           100

Q3plot(intervals(lmlist))

使用 lmLst 和间隔函数获取每位患者的截距和斜率（作为重复测量研究中的诊断）

Obtaining intercept and slope per patient (as diagnostics in a repeated measurements study) using the lmLst and intervals functions

r

linear-regression

lme4

mixed-models

confidence-interval