我需要帮助从我的线性混合模型中提取
yr_qun
变量的估计值和 CI。当我运行以下代码时,我得到
错误:
! Tibble 列必须具有兼容的大小。 • 规模 30:现有数据。 • 32 号:专栏
。 ℹ 只有大小一的值被回收。ci_lower
非常感谢任何帮助
library(lme4) # linear mixed-effects models
library(lmerTest) # test for linear mixed-effects models
library(gtsummary)
library(sjPlot)
library(ggplot2)
library(tidyverse)
library(tibble)
m1 <- lmer(distance ~ yr_qun+age+sex + (1 | id), data = trajectories)
summary(m1)
ci <- as.data.frame(confint(m1))
m1_estimate <- summary(m1)$coefficients[2:length(unique(trajectories$yr_qun)),1]
target <- setdiff(intersect(rownames(ci),rownames(summary(m1)$coefficients)),"(Intercept)")
m1_ci <- ci[target,]
summary(m1)
result <- tibble(name = names(m1_estimate),
est = as.numeric(m1_estimate),
ci_lower = as.numeric(m1_ci[,1]),
ci_upper = as.numeric(m1_ci[,2]))
这是我的虚拟数据
structure(list(id = c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 5L, 5L,
5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 6L, 6L, 6L, 6L, 6L, 6L, 7L, 7L,
7L, 7L, 7L, 8L, 8L, 8L, 8L, 8L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 10L,
10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 11L, 11L, 11L, 11L,
11L, 11L, 11L, 11L, 11L, 11L, 12L, 12L, 12L, 12L, 12L, 12L, 12L,
12L, 12L, 12L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L,
14L, 14L, 14L, 14L, 14L, 14L, 14L, 14L, 14L, 14L, 15L, 15L, 15L,
15L, 15L, 15L, 15L, 15L, 15L, 15L, 16L, 16L, 16L, 16L, 16L, 16L,
16L, 16L, 16L, 16L, 17L, 17L, 17L, 17L, 17L, 17L, 17L, 17L, 17L,
17L, 18L, 18L, 18L, 18L, 18L, 18L, 18L, 18L, 18L, 18L, 19L, 19L,
19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 20L, 20L, 20L, 20L, 20L,
20L, 20L, 20L, 20L, 20L), year = c(2001L, 2002L, 2003L, 2004L,
2005L, 2006L, 2007L, 2008L, 2009L, 2010L, 2001L, 2002L, 2003L,
2004L, 2005L, 2006L, 2007L, 2008L, 2009L, 2010L, 2001L, 2002L,
2003L, 2004L, 2005L, 2006L, 2007L, 2008L, 2009L, 2010L, 2001L,
2002L, 2003L, 2004L, 2005L, 2006L, 2007L, 2008L, 2009L, 2010L,
2001L, 2002L, 2003L, 2004L, 2005L, 2006L, 2007L, 2008L, 2009L,
2010L, 2001L, 2002L, 2003L, 2004L, 2005L, 2006L, 2001L, 2002L,
2003L, 2004L, 2005L, 2006L, 2007L, 2008L, 2009L, 2010L, 2004L,
2005L, 2006L, 2007L, 2008L, 2009L, 2010L, 2001L, 2002L, 2003L,
2004L, 2005L, 2006L, 2007L, 2008L, 2009L, 2010L, 2001L, 2002L,
2003L, 2004L, 2005L, 2006L, 2007L, 2008L, 2009L, 2010L, 2001L,
2002L, 2003L, 2004L, 2005L, 2006L, 2007L, 2008L, 2009L, 2010L,
2001L, 2002L, 2003L, 2004L, 2005L, 2006L, 2007L, 2008L, 2009L,
2010L, 2001L, 2002L, 2003L, 2004L, 2005L, 2006L, 2007L, 2008L,
2009L, 2010L, 2001L, 2002L, 2003L, 2004L, 2005L, 2006L, 2007L,
2008L, 2009L, 2010L, 2001L, 2002L, 2003L, 2004L, 2005L, 2006L,
2007L, 2008L, 2009L, 2010L, 2001L, 2002L, 2003L, 2004L, 2005L,
2006L, 2007L, 2008L, 2009L, 2010L, 2001L, 2002L, 2003L, 2004L,
2005L, 2006L, 2007L, 2008L, 2009L, 2010L, 2001L, 2002L, 2003L,
2004L, 2005L, 2006L, 2007L, 2008L, 2009L, 2010L, 2001L, 2002L,
2003L, 2004L, 2005L, 2006L, 2007L, 2008L, 2009L, 2010L), distance = c(15,
20, 21.5, 23, 21, 21.5, 24, 25.5, 20.5, 24, 21, 20, 21.5, 23,
21, 21.5, 24, 25.5, 20.5, 24, 21, 20, 21.5, 23, 21, 21.5, 24,
25.5, 20.5, 24, 21, 20, 21.5, 23, 21, 21.5, 24, 25.5, 20.5, 24,
21, 20, 21.5, 23, 21, 21.5, 24, 25.5, 20.5, 24, 21, 20, 21.5,
23, 21, 21.5, 21, 20, 21.5, 23, 21, 21.5, 24, 25.5, 20.5, 24,
23, 21, 21.5, 24, 25.5, 20.5, 24, 15, 20, 21.5, 23, 21, 21.5,
24, 25.5, 20.5, 24, 15, 20, 21.5, 23, 21, 21.5, 24, 25.5, 20.5,
24, 15, 20, 21.5, 23, 21, 21.5, 24, 25.5, 20.5, 24, 15, 20, 21.5,
23, 21, 21.5, 24, 25.5, 20.5, 24, 15, 20, 21.5, 23, 21, 21.5,
24, 25.5, 20.5, 24, 15, 20, 21.5, 23, 21, 21.5, 24, 25.5, 20.5,
24, 15, 20, 21.5, 23, 21, 21.5, 24, 25.5, 20.5, 24, 15, 20, 21.5,
23, 21, 21.5, 24, 25.5, 20.5, 24, 15, 20, 21.5, 23, 21, 21.5,
24, 25.5, 20.5, 24, 15, 20, 21.5, 23, 21, 21.5, 24, 25.5, 20.5,
24, 15, 20, 21.5, 23, 21, 21.5, 24, 25.5, 20.5, 24), age = c(8L,
9L, 10L, 11L, 12L, 13L, 14L, 15L, 16L, 17L, 11L, 12L, 13L, 14L,
15L, 16L, 17L, 18L, 19L, 20L, 15L, 16L, 17L, 18L, 19L, 20L, 21L,
22L, 23L, 24L, 24L, 25L, 26L, 27L, 28L, 29L, 30L, 31L, 32L, 33L,
6L, 7L, 8L, 9L, 10L, 11L, 12L, 13L, 14L, 15L, 9L, 10L, 11L, 12L,
13L, 14L, 6L, 7L, 8L, 9L, 10L, 11L, 12L, 13L, 14L, 15L, 18L,
19L, 20L, 21L, 22L, 23L, 24L, 28L, 40L, 41L, 42L, 43L, 44L, 45L,
46L, 47L, 48L, 28L, 29L, 30L, 31L, 32L, 33L, 34L, 35L, 36L, 37L,
28L, 29L, 30L, 31L, 32L, 33L, 34L, 35L, 36L, 37L, 28L, 29L, 30L,
31L, 32L, 33L, 34L, 35L, 36L, 37L, 28L, 29L, 30L, 31L, 32L, 33L,
34L, 35L, 36L, 37L, 28L, 29L, 30L, 31L, 32L, 33L, 34L, 35L, 36L,
37L, 28L, 29L, 30L, 31L, 32L, 33L, 34L, 35L, 36L, 37L, 28L, 29L,
30L, 31L, 32L, 33L, 34L, 35L, 36L, 37L, 28L, 29L, 30L, 31L, 32L,
33L, 34L, 35L, 36L, 37L, 28L, 29L, 30L, 31L, 32L, 33L, 34L, 35L,
36L, 37L, 28L, 29L, 30L, 31L, 32L, 33L, 34L, 35L, 36L, 37L),
Quintile = structure(c(5L, 2L, 3L, 3L, 2L, 2L, 4L, 2L, 5L,
5L, 1L, 4L, 2L, 5L, 4L, 3L, 3L, 4L, 3L, 3L, 1L, 3L, 1L, 2L,
1L, 5L, 2L, 4L, 1L, 4L, 3L, 2L, 5L, 3L, 4L, 4L, 3L, 1L, 4L,
3L, 4L, 1L, 4L, 4L, 5L, 1L, 5L, 2L, 2L, 2L, 3L, 5L, 3L, 3L,
4L, 1L, 3L, 1L, 1L, 5L, 2L, 4L, 1L, 3L, 2L, 4L, 1L, 3L, 4L,
5L, 3L, 3L, 1L, 2L, 2L, 3L, 1L, 2L, 3L, 5L, 5L, 2L, 5L, 2L,
2L, 3L, 1L, 2L, 3L, 5L, 5L, 2L, 5L, 2L, 2L, 3L, 1L, 2L, 3L,
5L, 5L, 2L, 5L, 2L, 2L, 3L, 1L, 2L, 3L, 5L, 5L, 2L, 5L, 2L,
2L, 3L, 1L, 2L, 3L, 5L, 5L, 2L, 5L, 2L, 2L, 3L, 1L, 2L, 3L,
5L, 5L, 2L, 5L, 2L, 2L, 3L, 1L, 2L, 3L, 5L, 5L, 2L, 5L, 2L,
2L, 3L, 1L, 2L, 3L, 5L, 5L, 2L, 5L, 2L, 2L, 3L, 1L, 2L, 3L,
5L, 5L, 2L, 5L, 2L, 2L, 3L, 1L, 2L, 3L, 5L, 5L, 2L, 5L, 2L,
2L, 3L, 1L, 2L, 3L, 5L, 5L, 2L, 5L), levels = c("Q1", "Q2",
"Q3", "Q4", "Q5"), class = "factor"), sex = structure(c(2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 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, 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, 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), levels = c("F", "M"), class = "factor"), yr_qun = structure(c(3L,
4L, 8L, 11L, 13L, 16L, 21L, 22L, 27L, 30L, 1L, 6L, 7L, 12L,
15L, 17L, 20L, 24L, 26L, 29L, 1L, 5L, 31L, 31L, 31L, 31L,
19L, 24L, 25L, 30L, 2L, 4L, 9L, 11L, 15L, 18L, 20L, 22L,
27L, 29L, 3L, 4L, 9L, 12L, 15L, 16L, 21L, 22L, 25L, 28L,
2L, 6L, 8L, 11L, 15L, 16L, 2L, 4L, 7L, 12L, 13L, 18L, 19L,
23L, 25L, 30L, 10L, 14L, 18L, 21L, 23L, 26L, 28L, 1L, 4L,
8L, 10L, 13L, 17L, 21L, 24L, 25L, 30L, 1L, 4L, 8L, 10L, 13L,
17L, 21L, 24L, 25L, 30L, 1L, 4L, 8L, 10L, 31L, 31L, 31L,
24L, 25L, 31L, 31L, 31L, 8L, 10L, 13L, 17L, 21L, 24L, 25L,
30L, 31L, 31L, 8L, 10L, 13L, 17L, 21L, 24L, 25L, 30L, 31L,
31L, 8L, 10L, 13L, 17L, 21L, 24L, 25L, 30L, 31L, 31L, 8L,
10L, 13L, 17L, 21L, 24L, 25L, 30L, 31L, 31L, 8L, 10L, 13L,
17L, 21L, 24L, 25L, 30L, 31L, 31L, 8L, 10L, 13L, 17L, 21L,
24L, 25L, 30L, 31L, 31L, 8L, 10L, 13L, 17L, 21L, 24L, 25L,
30L, 31L, 31L, 8L, 10L, 13L, 17L, 21L, 24L, 25L, 30L), levels = c("2001_low",
"2001_medium", "2001_high", "2002_low", "2002_medium", "2002_high",
"2003_low", "2003_medium", "2003_high", "2004_low", "2004_medium",
"2004_high", "2005_low", "2005_medium", "2005_high", "2006_low",
"2006_medium", "2006_high", "2007_low", "2007_medium", "2007_high",
"2008_low", "2008_medium", "2008_high", "2009_low", "2009_medium",
"2009_high", "2010_low", "2010_medium", "2010_high", ""), class = "factor")), class = "data.frame", row.names = c(NA,
-183L))
当你已经加载了整个
tidyverse
,也许只是加入系数和 CI:
library(lme4) # linear mixed-effects models
library(dplyr)
library(stringr)
m1 <- lmer(distance ~ yr_qun+age+sex + (1 | id), data = trajectories)
left_join(
as_tibble(summary(m1)$coefficients, rownames = "term"),
as_tibble(confint(m1), rownames = "term")
) %>% filter(str_starts(term,"yr_qun"))
#> Joining with `by = join_by(term)`
#> # A tibble: 30 × 6
#> term Estimate `Std. Error` `t value` `2.5 %` `97.5 %`
#> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 yr_qun2001_medium 3.81 1.01 3.76 1.95 5.61
#> 2 yr_qun2001_high 0.946 1.15 0.823 -1.17 3.00
#> 3 yr_qun2002_low 2.83 0.785 3.61 1.40 4.23
#> 4 yr_qun2002_medium 2.02 1.47 1.37 -0.604 4.74
#> 5 yr_qun2002_high 2.72 1.14 2.39 0.658 4.75
#> 6 yr_qun2003_low 4.13 1.16 3.57 2.03 6.18
#> 7 yr_qun2003_medium 4.56 0.714 6.39 3.24 5.84
#> 8 yr_qun2003_high 4.24 1.14 3.71 2.18 6.28
#> 9 yr_qun2004_low 6.04 0.728 8.29 4.70 7.34
#> 10 yr_qun2004_medium 5.94 0.998 5.95 4.11 7.72
#> # ℹ 20 more rows
创建于 2023-05-19 与 reprex v2.0.2
其他答案都不错,但有更简单的方法。
library(broom.mixed)
result <- (
tidy(m1, effects = "fixed", conf.int = TRUE, conf.method = "profile")
|> select(name = term, est = estimate, ci_lower = conf.low, ci_upper = conf.high)
|> filter(str_starts(name,"yr_qun"))
)
(你可以忽略关于“ddf.method ignored ...”的警告消息:这是应该在包中修复的东西...)
<chr> <dbl> <dbl> <dbl>
1 yr_qun2001_medium 3.81 1.81 5.81
2 yr_qun2001_high 0.946 -1.32 3.22
3 yr_qun2002_low 2.83 1.28 4.38
4 yr_qun2002_medium 2.02 -0.887 4.93
5 yr_qun2002_high 2.72 0.471 4.97
6 yr_qun2003_low 4.13 1.84 6.41
7 yr_qun2003_medium 4.56 3.15 5.97
8 yr_qun2003_high 4.24 1.98 6.50
9 yr_qun2004_low 6.04 4.60 7.48
10 yr_qun2004_medium 5.94 3.97 7.91
...
如果你看看你的组合,你应该看到问题:
> rownames(ci)
[1] ".sig01" ".sigma" "(Intercept)"
[4] "yr_qun2001_medium" "yr_qun2001_high" "yr_qun2002_low"
[7] "yr_qun2002_medium" "yr_qun2002_high" "yr_qun2003_low"
[10] "yr_qun2003_medium" "yr_qun2003_high" "yr_qun2004_low"
[13] "yr_qun2004_medium" "yr_qun2004_high" "yr_qun2005_low"
[16] "yr_qun2005_medium" "yr_qun2005_high" "yr_qun2006_low"
[19] "yr_qun2006_medium" "yr_qun2006_high" "yr_qun2007_low"
[22] "yr_qun2007_medium" "yr_qun2007_high" "yr_qun2008_low"
[25] "yr_qun2008_medium" "yr_qun2008_high" "yr_qun2009_low"
[28] "yr_qun2009_medium" "yr_qun2009_high" "yr_qun2010_low"
[31] "yr_qun2010_medium" "yr_qun2010_high" "yr_qun"
[34] "age" "sexM"
> names(m1_estimate)
[1] "yr_qun2001_medium" "yr_qun2001_high" "yr_qun2002_low"
[4] "yr_qun2002_medium" "yr_qun2002_high" "yr_qun2003_low"
[7] "yr_qun2003_medium" "yr_qun2003_high" "yr_qun2004_low"
[10] "yr_qun2004_medium" "yr_qun2004_high" "yr_qun2005_low"
[13] "yr_qun2005_medium" "yr_qun2005_high" "yr_qun2006_low"
[16] "yr_qun2006_medium" "yr_qun2006_high" "yr_qun2007_low"
[19] "yr_qun2007_medium" "yr_qun2007_high" "yr_qun2008_low"
[22] "yr_qun2008_medium" "yr_qun2008_high" "yr_qun2009_low"
[25] "yr_qun2009_medium" "yr_qun2009_high" "yr_qun2010_low"
[28] "yr_qun2010_medium" "yr_qun2010_high" "yr_qun"
> target
[1] "yr_qun2001_medium" "yr_qun2001_high" "yr_qun2002_low"
[4] "yr_qun2002_medium" "yr_qun2002_high" "yr_qun2003_low"
[7] "yr_qun2003_medium" "yr_qun2003_high" "yr_qun2004_low"
[10] "yr_qun2004_medium" "yr_qun2004_high" "yr_qun2005_low"
[13] "yr_qun2005_medium" "yr_qun2005_high" "yr_qun2006_low"
[16] "yr_qun2006_medium" "yr_qun2006_high" "yr_qun2007_low"
[19] "yr_qun2007_medium" "yr_qun2007_high" "yr_qun2008_low"
[22] "yr_qun2008_medium" "yr_qun2008_high" "yr_qun2009_low"
[25] "yr_qun2009_medium" "yr_qun2009_high" "yr_qun2010_low"
[28] "yr_qun2010_medium" "yr_qun2010_high" "yr_qun"
[31] "age" "sexM"
注意,
target
有元素age
和sexM
,它们不在m1
中。
您可以重新创建
target
,这样它只包含带有子字符串qun
的变量:
intersect(rownames(ci),rownames(summary(m1)$coefficients))[grepl( '*qun*', intersect(rownames(ci),rownames(summary(m1)$coefficients)))] -> target
现在:
m1_ci <- ci[target,]
result <- tibble(name = names(m1_estimate),
est = as.numeric(m1_estimate),
ci_lower = as.numeric(m1_ci[,1]),
ci_upper = as.numeric(m1_ci[,2]))
> result
# A tibble: 30 × 4
name est ci_lower ci_upper
<chr> <dbl> <dbl> <dbl>
1 yr_qun2001_medium 3.81 1.95 5.61
2 yr_qun2001_high 0.946 -1.17 3.00
3 yr_qun2002_low 2.83 1.40 4.23
4 yr_qun2002_medium 2.02 -0.604 4.74
5 yr_qun2002_high 2.72 0.658 4.75
6 yr_qun2003_low 4.13 2.03 6.18
7 yr_qun2003_medium 4.56 3.24 5.84
8 yr_qun2003_high 4.24 2.18 6.28
9 yr_qun2004_low 6.04 4.70 7.34
10 yr_qun2004_medium 5.94 4.11 7.72
# … with 20 more rows
# ℹ Use `print(n = ...)` to see more rows