I'm using the lme
function from the nlme
package and having a hard time trying to apply it to each column of my tibble. I can successfully run it 'manually' with a single column, but fail when trying to use purrr::map
. I feel like something is right under my nose, but I'm not seeing it. Any help would be appreciated.
library(tidyverse)
library(nlme)
# Create Tibble
df <- tibble(
id = rep(1:5, length = 30),
ko = as_factor(rep(c('ctrl', 'ko1', 'ko2', 'ko3', 'ko4', 'ko5'), each = 5)),
marker1 = rnorm(30),
marker2 = rnorm(30)
)
# 'Manual' calculation
df %>% lme(marker1 ~ ko, random = ~1|id, data = .)
# Attempt at calculating lme for each column
df %>% map_if(is_double, ~ lme(.x ~ ko, random = ~ 1|id, data = .))
# Ideal would be to get the tibble of each like:
df %>% summarize(broom.mixed::tidy(lme(marker1 ~ ko, random = ~1|id, data = .)))
CodePudding user response:
This is a little bit tricky because the formulas inside lme
use their own version of non-standard evaluation. I would suggest
vars <- names(df)[map_lgl(df, is_double)]
(vars
%>% setNames(vars) ## trick for naming results
%>% map(~reformulate("ko", response = .))
%>% map(lme, random = ~1|id, data = df)
%>% map_dfr(broom.mixed::tidy, effects = "fixed", .id = "var")
)
(I'm not sure what summarize()
is doing there?)
If you were using the lme4
package you could use the refit()
function to do this a bit more efficiently ...
CodePudding user response:
Because of how lme()
interprets its formula argument, I don't think this can be done across columns. However, we could pivot the data into a longer format so that the dependent variable is in one column, with another column indicating that these should be tested separately:
library(nlme)
library(tidyverse)
tests <- df %>%
pivot_longer(marker1:marker2, names_to = 'marker') %>%
nest(data = -marker) %>%
rowwise() %>%
mutate(
lme = list(lme(value ~ ko, random = ~ 1|id, data = data)),
result = list(broom.mixed::tidy(lme))
) %>%
unnest(result)
marker data lme effect group term estimate std.error df statistic p.value
<chr> <list> <list> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
1 marker1 <tibble> <lme> fixed fixed (Inte… 2.00e-1 0.394 20 0.506 0.618
2 marker1 <tibble> <lme> fixed fixed koko1 3.46e-1 0.558 20 0.621 0.542
3 marker1 <tibble> <lme> fixed fixed koko2 -7.70e-1 0.558 20 -1.38 0.183
4 marker1 <tibble> <lme> fixed fixed koko3 1.56e-2 0.558 20 0.0281 0.978
5 marker1 <tibble> <lme> fixed fixed koko4 1.76e-1 0.558 20 0.316 0.755
6 marker1 <tibble> <lme> fixed fixed koko5 -2.27e-1 0.558 20 -0.408 0.688
7 marker1 <tibble> <lme> ran_pars id sd_(I… 2.03e-5 NA NA NA NA
8 marker1 <tibble> <lme> ran_pars Residual sd_Ob… 8.82e-1 NA NA NA NA
9 marker2 <tibble> <lme> fixed fixed (Inte… 1.75e-1 0.414 20 0.423 0.677
10 marker2 <tibble> <lme> fixed fixed koko1 -1.46e-1 0.542 20 -0.269 0.791
11 marker2 <tibble> <lme> fixed fixed koko2 -4.89e-1 0.542 20 -0.902 0.378
12 marker2 <tibble> <lme> fixed fixed koko3 -6.79e-1 0.542 20 -1.25 0.225
13 marker2 <tibble> <lme> fixed fixed koko4 -5.70e-1 0.542 20 -1.05 0.306
14 marker2 <tibble> <lme> fixed fixed koko5 -5.57e-1 0.542 20 -1.03 0.317
15 marker2 <tibble> <lme> ran_pars id sd_(I… 3.48e-1 NA NA NA NA
16 marker2 <tibble> <lme> ran_pars Residual sd_Ob… 8.58e-1 NA NA NA NA