plot estimates with 95%CI from lmer model

Question

I need help in plotting the estimates and 95% CI from Lmer model, i am using plot_model from sjplot function.

library(lme4)       # linear mixed-effects models
library(lmerTest)   # test for linear mixed-effects models
library(gtsummary)
library(sjPlot)
library(ggplot2)

names(trajectories)
library(tidyverse)


str(trajectories$yr_qun)

trajectories <- trajectories %>% 
  mutate(yr_qun = yr_qun %>%
           fct_relevel("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"))

m1 <- lmer(distance ~ yr_qun + (1 | id), data = trajectories)
summary(m1)

p <- plot_model(m1, order.terms = rev(1:29)) + coord_cartesian() 
p

i would like to have a plot as shown in the following picture

here is dummy data

dput(trajectories)
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))

any help much appreciated

score 1 · Accepted Answer · edited May 07 '23 at 05:27

1

First, You could extract estimate and 95% CI, then create a data.frame for a plot.
Next, year and class was extracted from yr_qun~, using str_extract().

library(lmerTest)   # test for linear mixed-effects models
library(gtsummary)

library(ggplot2)
library(dplyr)
library(tidyr)
library(forcats)
library(stringr)
library(ggplot2)

m1 <- lmer(distance ~ yr_qun + (1 | id), data = trajectories)
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,]

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 <- result |> 
  filter(name != "yr_qun") |> 
  mutate(year = str_extract(name,"[:digit:]+"),
         class = str_extract(name, "[:alpha:]+$"))

result |> 
  ggplot(aes(x=year, y=est, color = class, group=class))+
  geom_line()+
  geom_point()+
  geom_errorbar(aes(ymin = ci_lower, ymax = ci_upper),
                alpha=.5) +
  theme_minimal()

^{Created on 2023-05-04 with reprex v2.0.2}

edited May 07 '23 at 05:27

skpak

79
6

answered May 04 '23 at 00:38

YH Jang

1,306
5
15

My pleasure. If you like the answer, please accept it :) – YH Jang May 04 '23 at 08:59
sure I'll accept it shortly. just running on real data. – skpak May 04 '23 at 09:32
@ Y H Jang, i want to clarify the dimension in the coefficients in m1_estimate <- summary(m1)$coefficients[2:31,1], in my real data i have 30 levels in variable so i changed to this m1_estimate <- summary(m1)$coefficients[2:30,1], but i get an error of incorrect number of dimensions. any help would be highly appreciated. thank you – skpak May 04 '23 at 09:50
I see. How about `summary(m1)$coefficients[2:length(unique(trajectories$yr_qun) ),1]`? Then you don't have to manually set the index – YH Jang May 04 '23 at 10:08
thanks that worked. thank you. moving ahead with code – skpak May 04 '23 at 10:31
got some connection issues so unable to run the rest of the code. so I'll update you once my connection to real data is restored. thanks for all your help. – skpak May 04 '23 at 10:44
hi got stuck again with code 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])) when i run this code, i get an error of Error in tibble() tibble columns must have compatible sizes. size 145: existing data. size29: column 'ci_lower'. only values of size one are recycled. Backtrace: 1.tibble::tibble(...) – skpak May 04 '23 at 11:23
I edited the answer. Check this out. – YH Jang May 04 '23 at 11:57
many thanks for the prompt help and edit. i updated my code with new ci, target and m1_ci but still get the same error i get an error of Error in tibble() tibble columns must have compatible sizes. size 145: existing data. size30: column 'ci_lower'. only values of size one are recycled. Backtrace: 1.tibble::tibble(...) – skpak May 04 '23 at 12:36
Seems like the length of data problem. Unfortunately, I couldn't help you without your real data. – YH Jang May 05 '23 at 00:38
thank you much for all your help. Due to sensitivety of the data even I don't have full access to data. I ll try to figure out. – skpak May 05 '23 at 06:47
@ YH Jang, I am trying to figure out where 145: existing data comes from. I checked the size for both m1_estimate and ci_lower used in Tibble, both have the same size of 29. how should i figure it out. – skpak May 05 '23 at 08:50
Have you check the `m1_estimate`? I guess the error is to do with the estimate. What you only need in `summary(m1)` is the `Estimate` column. 145 means maybe you are binding all `summary(m1)` components (i.e., not only `Estimate` but `std.error`, `df`, `t value`, and `p`) – YH Jang May 05 '23 at 23:06
@ YH Jang, you are right. the 145 figure was binding all components. I modified the code as below. result <- tibble(name = names(m1_estimate[,1]), est = as.numeric(m1_estimate[,1]), ci_lower = as.numeric(m1_ci[,1]), ci_upper = as.numeric(m1_ci[,2])) – skpak May 06 '23 at 11:01
thanks it worked. I realize regression output changed the variable label in the summary output. thought it is a factor variable and the same categories as in the variable yr_qun above (see above screenshot). what should I do so my regression output doesn't change the labels? – skpak May 06 '23 at 11:40
i need your help to extract some info from the names variable in the above codes for extracting year and class from yr_qun~, using str_extract(). in the real data, the the text in names column is "intrcn.yr.highexposurecat2_revised2013 low exposure "Some examples of text are attached in the picture above. I have to extract the last two (low exposure or high exposure) or one word ( middle) and second I want to extract the year from this mutate(year = str_extract(name,"[:digit:]+"), class = str_extract(name, "[:alpha:]+$")). any help much appreciated – skpak May 08 '23 at 23:32
Why don't you ask as another question? It's quite hard to answer by your comments. Plus asking multiple questions is not allowed in stackoverflow. – YH Jang May 09 '23 at 00:03
thanks YH Jang, i asked separate question https://stackoverflow.com/questions/76206568/extracting-last-word-and-year-using-str-extract-in-r – skpak May 09 '23 at 06:52
I just answered to your new question. And I would appreciate if you accept the answer of this question. – YH Jang May 09 '23 at 07:05
thank you so much for all your help. have a wonderful day ahead – skpak May 09 '23 at 07:59
how can i add asteric to estimates in the plot? – skpak May 18 '23 at 12:03

plot estimates with 95%CI from lmer model

1 Answers1