R 多款多功能

Question

我编写了一个例程，从 lmer 模型中提取信息来计算 ICC 并从 lmerTest 的 Ranova 函数中获取 LRT。我在下面的工作但我怀疑它可以通过(a)将两个函数合并为一个并返回一个列表来改进，但我似乎无法使用 purrr 的 map 函数访问列表元素，以及(b)使用多个 mutate/purrr 行将所有需要的数据集中在一个地方，而不必稍后加入。我的代码使用 Hox (2002) 中提供的“Peet”数据集，可在 UCLA IDRE 站点获得:

library(foreign)
library(lme4)
library(tidyverse)
library(purrr)


#Peet family data described and used in Hox
peet.dat<-read.dta("https://stats.idre.ucla.edu/stat/stata/examples/mlm_ma_hox/peetmis.dta")

names(peet.dat)

#convert to long format
peet.long.dat <- peet.dat %>%
  tidyr::gather(type, score, -family,-sex,-person) %>%
  arrange(type)

names(peet.long.dat)

#need two functions, one for the MLM estimates and the other for
#ranova p-test for variance--merge later by type

aov_model <- function(df) {
  lmr.model <- lmerTest::lmer(score~ 1 + (1|family), data=df)
}

aov_test <- function(df) {
  lmr.model <- lmerTest::lmer(score~ 1 + (1|family), data=df)
  ll.test <- lmerTest::ranova(lmr.model)
}

#get the model estimates
models <- peet.long.dat %>%
  nest(-type) %>%
  mutate(aov_obj = map(data, aov_model),
         summaries = map(aov_obj, broom.mixed::tidy)) %>%
  unnest(summaries, .drop = T) %>%
  select(type, effect, estimate, term) %>%
  filter(effect != "fixed") %>%
  mutate(variance = estimate^2) %>%
  select(-estimate, -effect) %>%
  spread(term, variance) %>%
  rename(group.var = `sd__(Intercept)`, residual = `sd__Observation`) %>%
  mutate(ICC = group.var/(group.var+residual))


models

#get the ranova LRTs
tests <- peet.long.dat %>%
  nest(-type) %>%
  mutate(test_obj = map(data, aov_test),
         test_summaries = map(test_obj, broom.mixed::tidy)) %>%
  unnest(test_summaries, .drop = T) %>%
  filter(!is.na(LRT))

#join estimates with LRT p values
models %>% left_join(tests[c("type","p.value")])

非常感谢任何帮助。

Answer 1

我认为这里的关键是split()您的 data.frame 基于变量 type :

# convert to list by type
peet.ls <- peet.dat %>%
  tidyr::gather(type, score, -family,-sex,-person) %>%
  split(.$type)

# map to fit models on subsets and return summaries
peet.ls %>%
  map(function(df.x) {
    # fit the model
    lmr_model <- lmerTest::lmer(score~ 1 + (1|family), data = df.x)
    
    #get the model estimates
    mlm_est <- lmr_model %>% 
      broom.mixed::tidy() %>%
      select(effect, estimate, term) %>%
      filter(effect != "fixed") %>%
      mutate(variance = estimate^2) %>%
      select(-estimate, -effect) %>%
      spread(term, variance) %>%
      rename(group.var = `sd__(Intercept)`, 
             residual = `sd__Observation`) %>%
      mutate(ICC = group.var/(group.var+residual))
    
    # get the ranova LRTs & add to other estimates
    mlm_est$p.value <- lmr_model %>%
      lmerTest::ranova() %>%
      broom.mixed::tidy() %>%
      filter(!is.na(LRT)) %>%
      pull(p.value)
  
    # return summaries
    mlm_est
  }) %>%
  # combine data.frames and add the variable 'type'
  bind_rows(.id = "type") %>%
  select(type, everything())