R多种型号具有多种功能

Question

我编写了一个例程，从 lmer 模型中提取信息以计算 ICC 并从 lmerTest 的 ranova 函数中获取 LRT。 我在下面拥有的内容有效，但我怀疑可以通过 （a） 将两个函数合并为一个并返回一个列表来改进它，但我似乎无法使用 purrr 的地图函数访问列表元素，以及 （b） 使用多个变异/咕噜声线在一个地方获取所有需要的数据，而不必稍后加入。 我的代码使用Hox（2002）中提供的"Peet"数据集进行操作，该数据集可在加州大学洛杉矶分校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

我认为这里的关键是根据变量type split()您的数据帧：

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