R много моделей несколько функций - PullRequest
Новая
       39

R много моделей несколько функций

0 голосов
/ 28 января 2019

Я написал процедуру, которая извлекает информацию из моделей lmer для вычисления ICC и получения LRT из функции ranova lmerTest.То, что у меня есть ниже, работает, но я подозреваю, что это может быть улучшено путем (а) объединения двух функций в одну и возврата списка, но я не могу получить доступ к элементам списка с помощью функции карты мурлыкания, и (б) с использованием множественного преобразования/ purrr строки, чтобы получить все необходимые данные в одном месте, а не присоединяться позже.Мой код следует с использованием набора данных «Peet», предоставленного в Hox (2002) и доступного на сайте 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")])

Любая помощь очень ценится.

...