rggplot2statisticssjplot
plot estimates with 95%CI from lmer model
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
Solution
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
- Test whether one factor is nested in another
- Renaming columns in the global environment
- How can I use a character containing the DB instead of the connectin in R?
- Identifying all columns shared in common by a set of files and then vertically stack them
- curly-curly braces - why won't this function accept an unquoted variable
- Concatenate columns from Excel file based on sheet name and column's name
- Conditioning previous values within groups in R
- Finding out which columns are same in different files
- Rcpp error, learning Rcpp and C++ on mac
- How can I subset a list, based on list element details, using sapply?
- Reactively pull data from list objects in Shiny
- read_xml not recognizing the contents inside the nodes and returning charater[0]
- Saving plots in multiple columns using ggsave() cuts 2/3 of output
- Replacing NA's within a set of similar values
- Hide initially displayed code from shiny modularized app after clicking actionButton
- edit plot attributes manually and redraw
- Reformat data to include 2 minute time intervals
- ggsurvplot not displaying number at risk for time 0
- Change width of callout-block in Quarto dashboard
- Is there a R function to pair elements from two columns of a data.frame to make a new data.frame
- Question about glm family in R
- How to flatten a list of lists?
- Is there a fast/clever way to return a logical vector if elements of a vector are in at least one interval?
- how can I combine three tables - 2 separate questionnaires used on the same cohort and one table with blood results
- Modify pop-up message in leaflet map
- Conditional bslib tooltip
- How to specify column parameters programmatically?
- Complex subsetting in R for case-control study
- Set CMU Sans Serif font on pdf quarto
- How to remove string before a certain delimiter in R?