Linear mixed effect model with group

I'm trying to fit a lmer model using dplyr::group_by to not fit the model for each of my species separately.

I found this line of code that seems to work but, I don't know how to visualize the results.

library("lme4")
data(Orthodont,package="nlme")

ort_test <- Orthodont %>% group_by(Sex) %>%
  do(model = lmer(.,formula=distance~age (1|Subject)))

and this is what I get out of this

# A tibble: 2 × 2
# Rowwise: 
  Sex    model     
  <fct>  <list>    
1 Male   <lmrMdLmT>
2 Female <lmrMdLmT>

Can you help me to get the info from the ort_test$model column?

Thanks!!!

CodePudding user response：

We could use tidy from broom.mixed package

library(tidyr)
library(dplyr)
ort_test %>%
   mutate(out = list(broom.mixed::tidy(model))) %>%
   ungroup %>% 
   select(Sex, out) %>% 
   unnest(out)

-output

# A tibble: 8 × 7
  Sex    effect   group    term            estimate std.error statistic
  <fct>  <chr>    <chr>    <chr>              <dbl>     <dbl>     <dbl>
1 Male   fixed    <NA>     (Intercept)       16.3      1.13       14.5 
2 Male   fixed    <NA>     age                0.784    0.0938      8.36
3 Male   ran_pars Subject  sd__(Intercept)    1.63    NA          NA   
4 Male   ran_pars Residual sd__Observation    1.68    NA          NA   
5 Female fixed    <NA>     (Intercept)       17.4      0.859      20.2 
6 Female fixed    <NA>     age                0.480    0.0526      9.12
7 Female ran_pars Subject  sd__(Intercept)    2.07    NA          NA   
8 Female ran_pars Residual sd__Observation    0.780   NA          NA

CodePudding user response：

The new reframe function from dplyr v1.1 is a perfect fit for doing this kind of work. Read more about it here. Assuming that you're only interested in the model coefficients, you could do the following:

Orthodont |> 
  group_by(Sex) |> 
  reframe(
    lmer(distance ~ age   (1 | Subject)) |> 
      summary() |> 
      (`$`)("coefficients") |> 
      as_tibble(rownames = "term", .name_repair = janitor::make_clean_names)
  )
#> # A tibble: 4 × 5
#>   sex    term        estimate std_error t_value
#>   <fct>  <chr>          <dbl>     <dbl>   <dbl>
#> 1 Male   (Intercept)   16.3      1.13     14.5 
#> 2 Male   age            0.784    0.0938    8.36
#> 3 Female (Intercept)   17.4      0.859    20.2 
#> 4 Female age            0.480    0.0526    9.12