Попытка использовать tidy для анализа мощности и использование clmm2

Question

Я пытаюсь провести анализ мощности на анализе clmm2, который я делаю. Это код конкретной статистической модели:

test <- clmm2(risk_sensitivity ~ treat + sex + dispersal + 
sex*dispersal + treat*dispersal + treat*sex,random = id, data = datasocial, Hess=TRUE)

Теперь у меня есть следующая функция:

sim_experiment_power <- function(rep) {
  s <- sim_experiment(n_sample = 1000,
                      prop_disp = 0.10,
                      prop_fem = 0.35,
                      disp_probability = 0.75,
                      nondisp_probability = 0.90,
                      fem_probability = 0.75,
                      mal_probability = 0.90)
  broom.mixed::tidy(s) %>%
    mutate(rep = rep)
}
my_power <- map_df(1:10, sim_experiment_power)

Детали функции sim_experiment не важны, потому что они работают должным образом. Важно знать, что он дает статистический результат clmm2. Моя цель с функцией выше - провести анализ мощности. Однако я получаю следующую ошибку:

Ошибка: нет аккуратного метода для объектов класса clmm2

Я немного новичок в R, но думаю, это означает, что tidy не работает с clmm2. Кто-нибудь знает способ решения этой проблемы?

РЕДАКТИРОВАТЬ: это то, что следует за кодом, который я опубликовал выше, и в конечном итоге это то, что я пытаюсь получить.

Затем вы можете построить график распределения оценок по вашему моделированию.

ggplot(my_power, aes(estimate, color = term)) +
  geom_density() +
  facet_wrap(~term, scales = "free")

Вы также можете просто вычислить мощность как долю p-значений меньше вашего альфа.

my_power %>%
  group_by(term) %>%
  summarise(power <- mean(p.value < 0.05))

Answer 1

Для того, что вам нужно, вы можете написать функцию для возврата коэффициентов с тем же именем столбца:

library(ordinal)
library(dplyr)
library(purrr)

tidy_output_clmm = function(fit){
  results = as.data.frame(coefficients(summary(fit)))
  colnames(results) = c("estimate","std.error","statistic","p.value")
  results %>% tibble::rownames_to_column("term")
}

Затем мы применяем его, используя пример, в котором я пробую набор данных вина по порядковому номеру:

sim_experiment_power <- function(rep) {
  idx = sample(nrow(wine),replace=TRUE)
  s <- clmm2(rating ~ temp, random=judge, data=wine[idx,], nAGQ=10,Hess=TRUE)
  tidy_output_clmm(s) %>% mutate(rep=rep)
}

my_power <- map_df(1:10, sim_experiment_power)

Земельные работы:

ggplot(my_power, aes(estimate, color = term)) +
  geom_density() +
  facet_wrap(~term, scales = "free")

h ttps:https://stackru.com/images/6f3f13c21cfd5f6e5128b120862de3e07a88e1bd.png

И власть тоже:

my_power %>% group_by(term) %>% summarise(power = mean(p.value < 0.05))
    # A tibble: 5 x 2
  term     power
  <chr>    <dbl>
1 1|2        0.9
2 2|3        0.1
3 3|4        1  
4 4|5        1  
5 tempwarm   1