R语言从线性混合模型中提取选定的估计值和CI

qcuzuvrc 于 2023-05-20 发布在其他

关注(0)|答案(3)|浏览(136)

我需要帮助从我的线性混合模型中提取yr_qun变量的估计值和CI。当我运行下面的代码时，我得到
错误：
！Tibble列必须具有兼容的大小。·大小30：现有数据。·大小32：列xlmlnlx。仅回收大小为1的值。
任何帮助非常感谢

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))

来源：https://stackoverflow.com/questions/76288227/extracting-selected-estimate-and-ci-from-a-linear-mixed-model

3条答案

按热度按时间

vaqhlq811#

因为你已经加载了整个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

赞(0）回复(0）举报 2023-05-20

zdwk9cvp2#

其他答案都很好，但有一个更简单的方法来做到这一点。

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"))
)

(you可以忽略有关“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
...

赞(0）回复(0）举报 2023-05-20

6ljaweal3#

如果你看看你的组合，你应该看到的问题：

> 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

赞(0）回复(0）举报 2023-05-20

我来回答

R语言从线性混合模型中提取选定的估计值和CI

3条答案

相关问题

热门标签

最新问答

R语言 从线性混合模型中提取选定的估计值和CI

3条答案

相关问题

热门标签

最新问答

R语言从线性混合模型中提取选定的估计值和CI