R多模型多函数

Question

我编写了一个例程，从lmer模型中提取信息来计算ICC，并从lmerTest的ranova函数中获得LRT。下面的内容可以工作，但我怀疑它可以通过以下方法得到改进：(a)将两个函数组合为一个函数并返回一个列表，但我似乎无法使用purrr的map函数访问list元素，(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

# 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())