`emmeans` for a gls model does not run inside `map`

Esta pregunta se inspira en no puede utilizar emmeans dentro del mapa, y relacionados con Mapa `joint_tests` a una lista después de encajar un modelo 'gls' y `group_by` y mantener los niveles de agrupación como el nombre del marco de datos anidado

Quiero envolver varias pruebas en un flujo de trabajo.

Esto funciona para un glm modelo.

library(dplyr)
library(purrr)
library(emmeans)
library(nlme)

diamond_result <- diamonds %>%
  group_by(cut) %>%
  nest() %>% 
  ungroup %>%
  dplyr::mutate(models=map(data,~glm(price ~ x + y + z + clarity + color,data=.x)),
         jt = map(models, ~emmeans::joint_tests(.x, data = .x$data)),
         means=map(models,~emmeans::emmeans(.x,"color",data=.x$data)),
         p_cont = map(means, ~emmeans::contrast(.x, "pairwise",infer = c(T,T))),
         across(models:p_cont, stats::setNames,  .$cut)) 
> diamond_result$jt
$Ideal
 model term df1 df2 F.ratio p.value
 x            1 Inf 611.626 <.0001 
 y            1 Inf   2.914 0.0878 
 z            1 Inf 100.457 <.0001 
 clarity      7 Inf 800.852 <.0001 
 color        6 Inf 256.796 <.0001 


$Premium
 model term df1 df2  F.ratio p.value
 x            1 Inf 2074.371 <.0001

Pero la misma sintaxis no funciona para una gls modelo así que me detuve emmeans() Paso. Eventualmente, quiero joint_tests, emmeans, y contrast en el mutate paso.

diamonds_emm2 <-  diamonds %>%
  group_by(cut) %>% 
  nest() %>%
  ungroup() %>%
  dplyr::mutate(models=map(data,~gls(price ~ x + y + z + clarity,  
                              weights = varIdent(form = ~ 1|color),
                              data =.x)),
                means=map(models,~emmeans::emmeans(.x,"clarity",data=.x$data)),
                across(models:p_cont, setNames,  .$cut))
    Error: Problem with `mutate()` input `means`.
    x undefined columns selected
    ℹ Input `means` is `map(models, ~emmeans::emmeans(.x, "clarity", data = .x$data))`.
    Run `rlang::last_error()` to see where the error occurred.

Problem with `mutate()` input `means`.
x undefined columns selected
ℹ Input `means` is `map(models, ~emmeans::emmeans(.x, "clarity", data = .x$data))`.
Backtrace:
  1. `%>%`(...)
 18. base::.handleSimpleError(...)
 19. dplyr:::h(simpleError(msg, call))


Problem with `mutate()` input `means`.
x undefined columns selected
ℹ Input `means` is `map(models, ~emmeans::emmeans(.x, "clarity", data = .x$data))`.
Backtrace:
     █
  1. ├─`%>%`(...)
  2. ├─dplyr::mutate(...)
  3. ├─dplyr:::mutate.data.frame(...)
  4. │ └─dplyr:::mutate_cols(.data, ...)
  5. │   ├─base::withCallingHandlers(...)
  6. │   └─mask$eval_all_mutate(dots[[i]])
  7. ├─purrr::map(models, ~emmeans::emmeans(.x, "clarity", data = .x$data))
  8. │ └─.f(.x[[i]], ...)
  9. │   └─emmeans::emmeans(.x, "clarity", data = .x$data)
 10. │     ├─base::do.call(ref_grid, args)
 11. │     └─(function (object, at, cov.reduce = mean, cov.keep = get_emm_option("cov.keep"), ...
 12. │       ├─emmeans::recover_data(object, data = as.data.frame(data), ...)
 13. │       └─emmeans:::recover_data.gls(...)
 14. │         └─emmeans:::recover_data.call(...)
 15. │           ├─tbl[, vars, drop = FALSE]
 16. │           └─base::`[.data.frame`(tbl, , vars, drop = FALSE)
 17. │             └─base::stop("undefined columns selected")
 18. └─base::.handleSimpleError(...)
 19.   └─dplyr:::h(simpleError(msg, call))

El código funciona bien en este paso.

diamonds_emm <-  diamonds %>%
  group_by(cut) %>% nest() %>%
  mutate(models=map(data,~gls(price ~ x + y + z + clarity,  
                              weights = varIdent(form = ~ 1|color),
                              data =.x)))

¿Cómo trabajo alrededor de este problema? Gracias.

Actualización: map2 función de la respuesta de Ronak solucionó el problema en el means pasos pero no hará contrastes pares. ¿Qué me perdí?

diamonds %>%
  group_by(cut) %>% 
  nest() %>%
  mutate(models=map(data,~gls(price ~ x + y + z + clarity,  
                              weights = varIdent(form = ~ 1|color),
                              data =.x)), 
         means = map2(data, models,~emmeans::emmeans(.y,"clarity",data=.x)),
         p_cont = map2(means, ~emmeans::contrast(.y, "pairwise",infer = c(T,T)))) %>%
  ungroup %>%
  mutate(across(models:p_cont, setNames,  .$cut)) -> result
Error: Problem with `mutate()` input `p_cont`.
x object '.z' not found
ℹ Input `p_cont` is `map(means, ~emmeans::contrast(.y, "pairwise", infer = c(T, T)))`.
ℹ The error occurred in group 1: cut = "Fair".

Dar un nuevo nombre a la entrada en el p_cont paso, como ~emmeans::contrast(.z, "pairwise", infer = c(T, T))) No solucionó el problema.

Pregunta hecha hace 3 años, 4 meses, 29 días - Por codecraftsman

4 Respuestas:

  • Datos de pase así como modelo en el emmeans paso utilizando map2. Para contrasts y joint_tests puedes usar map.

    library(tidyverse)
library(emmeans)
library(nlme)

diamonds %>%
  group_by(cut) %>% 
  nest() %>%
  mutate(models=map(data,~gls(price ~ x + y + z + clarity,  
                              weights = varIdent(form = ~ 1|color),
                              data =.x))) %>%
  ungroup %>%
  mutate(means = map2(data, models,~emmeans(.y,"clarity",data=.x)),
         p_cont = map(means, contrast, "pairwise"), 
         joint_tests = map(means, joint_tests),
         across(models:joint_tests, setNames,  .$cut)) -> result

result

#   cut       data                models      means       p_cont     joint_tests           
#                                  
#1 Ideal

    Respondida el Dec 17, 2020 a las 19:28 - por scriptsorcererf493

    Votos positivos: 0 | Votos negativos: 0

  • Si no insistes en usar el purrr::map familia, sugeriría usar el nuevo (dplyr 1.0.0) rowwise estilo. Es menos confuso, ya que puede utilizar los nombres de variable/columna como es y no hay necesidad de elegir el correcto map función y averiguar la notación de la lambda. Sólo tienes que envolver la llamada de la función list(...). Sólo la última llamada across necesita ser llamado a los no agrupados data.frame.

    library(tidyverse)
library(emmeans)
library(nlme)

diamonds_emm_row <- diamonds %>%
  nest_by(cut) %>% 
  dplyr::mutate(models = list(gls(price ~ x + y + z + clarity,  
                                  weights = varIdent(form = ~ 1|color),
                                  data = data)),
                jt = list(joint_tests(models, data = data)),
                means = list(emmeans::emmeans(models, "clarity", data = data)),
                p_cont = list(emmeans::contrast(means, "pairwise", infer = c(T,T)))) %>% 
  ungroup %>% 
  mutate(across(models:p_cont, setNames, .$cut))

diamonds_emm_row
#> # A tibble: 5 x 6
#>   cut                   data models      jt                 means     p_cont    
#>                 
#> 1 Fair           [1,610 × 9]          
#> 2 Good           [4,906 × 9]          
#> 3 Very Go…      [12,082 × 9]          
#> 4 Premium       [13,791 × 9]          
#> 5 Ideal         [21,551 × 9]

    Creado el 2021-01-01 por el paquete de reprex (v0.3.0)

    Esto produce más o menos el mismo resultado que el uso purrr::map. "Más o menos", porque identical() no lo mostrará, ya que la llamada de función se guarda en los atributos (y otros lugares), y difiere dependiendo de si utiliza {dplyr}'s rowwise estilo o map notación de lambda.

    diamonds_emm_map <- diamonds %>%
  nest_by(cut) %>% 
  ungroup() %>%
  dplyr::mutate(models=map(data,~gls(price ~ x + y + z + clarity,  
                                     weights = varIdent(form = ~ 1|color),
                                     data =.x)),
                jt = map2(models, data, ~ joint_tests(.x, data = .y)),
                means = map2(data, models, ~ emmeans::emmeans(.y, "clarity", data = .x)),
                p_cont = map(means, ~emmeans::contrast(.x, "pairwise", infer = c(T,T))),
                across(models:p_cont, setNames, .$cut))

map2(diamonds_emm_row, diamonds_emm_map, all.equal, check.attributes = FALSE)
#> $cut
#> [1] TRUE
#> 
#> $data
#> [1] TRUE
#> 
#> $models
#> [1] "Component \"Fair\": Component \"call\": target, current do not match when deparsed"     
#> [2] "Component \"Good\": Component \"call\": target, current do not match when deparsed"     
#> [3] "Component \"Very Good\": Component \"call\": target, current do not match when deparsed"
#> [4] "Component \"Premium\": Component \"call\": target, current do not match when deparsed"  
#> [5] "Component \"Ideal\": Component \"call\": target, current do not match when deparsed"    
#> 
#> $jt
#> [1] TRUE
#> 
#> $means
#> [1] TRUE
#> 
#> $p_cont
#> [1] TRUE

    Creado el 2021-01-01 por el paquete de reprex (v0.3.0)

    Respondida el Dec 17, 2020 a las 19:37 - por bytebarde61c

    Votos positivos: 0 | Votos negativos: 0

  • El problema que estás experimentando con emmeans dentro de map cuando usas un modelo gls se debe a que las funciones emmeans no están tratando adecuadamente los datos debido a la estructura de los objetos generados por el modelo gls.

    Para solucionar este problema y permitir que las funciones emmeans funcionen correctamente con modelos gls dentro de un map, puedes hacer lo siguiente:

    • Para el cálculo de las medias (means), puedes usar `map2` en lugar de `map`. Esta función te permite aplicar una función a dos listas diferentes al mismo tiempo. En este caso, puedes combinar los datos y los modelos para asegurarte de que las medias se calculen correctamente:
    R
diamonds_emm2 <-  diamonds %>%
  group_by(cut) %>% 
  nest() %>%
  ungroup() %>%
  dplyr::mutate(models=map(data,~gls(price ~ x + y + z + clarity,  
                              weights = varIdent(form = ~ 1|color),
                              data =.x)),
                means=map2(data, models, ~emmeans::emmeans(.y, "clarity", data = .x$data)),
                across(models:means, setNames,  .$cut))
    • Para el cálculo de los contrastes (p_cont), también debes utilizar `map2` de manera similar para asegurarte de que los contrastes se calculen correctamente:
    R
diamonds_emm2 <-  diamonds_emm2 %>%
  mutate(p_cont = map2(data, means, ~emmeans::contrast(.y, "pairwise", infer = c(TRUE, TRUE))),
         across(models:p_cont, setNames,  .$cut))

    Al utilizar map2 en lugar de map para las funciones emmeans y contrast, puedes garantizar que se apliquen correctamente a los datos y modelos dentro de un flujo de trabajo en R. Espero que esta solución te ayude a resolver el problema que estás experimentando.

    Respondida el Dec 18, 2020 a las 19:30 - por Gpt

    Votos positivos: 0 | Votos negativos: 0

  • Problema con emmeans() para un modelo gls dentro de map()

    El problema surge porque emmeans() espera un objeto de modelo lm o glm como entrada, mientras que gls() produce un objeto de clase gls.

    Solución:

    Para emmeans(): Utiliza la función recover_data() para convertir el objeto gls en un marco de datos antes de aplicar emmeans().

    
means = map2(data, models, ~emmeans::emmeans(.y, "clarity", data = emmeans::recover_data(.x)))

    Para contrast():

    Opción 1: Utiliza emmeans() %>% contrast() para realizar contrastes directamente en los resultados de emmeans().

    
p_cont = map(means, ~.x %>% emmeans::contrast("pairwise", infer = c(T, T)))

    Opción 2: Utiliza la función lsmeans() del paquete lsmeans para realizar contrastes en objetos gls.

    
library(lsmeans)

p_cont = map(means, ~lsmeans(.x, "clarity") %>% contrast(pairs = TRUE, adjust = "tukey"))

    Ejemplo completo:

    
library(dplyr)
library(purrr)
library(emmeans)
library(nlme)
library(lsmeans)

diamond_result <- diamonds %>%
  group_by(cut) %>%
  nest() %>%
  ungroup %>%
  dplyr::mutate(
    models=map(data,~gls(price ~ x + y + z + clarity,  
                          weights = varIdent(form = ~ 1|color),
                          data =.x)),
    means = map2(data, models,~emmeans::emmeans(.y,"clarity",data=emmeans::recover_data(.x))),
    p_cont = map(means,~lsmeans(.x, "clarity") %>% contrast(pairs = TRUE, adjust = "tukey")),
    across(models:p_cont, setNames,  .$cut)
  )

    Nota: Si utilizas la opción 2 para los contrastes, necesitarás instalar y cargar el paquete lsmeans antes de ejecutar el código.

    Respondida el Dec 18, 2020 a las 19:34 - por Gemini

    Votos positivos: 0 | Votos negativos: 0