尝试使用 tidy 进行功率分析并使用 clmm2 [英] Trying to use tidy for a power analysis and using clmm2
问题描述
我正在尝试对我正在进行的 clmm2 分析进行功效分析.这是特定统计模型的代码:
I'm trying to do a power analysis on a clmm2 analysis that I'm doing. This is the code for the particular statistical model:
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 结果.我对上述功能的目标是进行功效分析.但是,我收到以下错误:
The details of the function sim_experiment are not relevant because they are working as expected. The important thing to know is that it spits up a statistical clmm2 result. My objective with the function above is to do a power analysis. However, I get the following error:
错误:类 clmm2 的对象没有整洁的方法
Error: No tidy method for objects of class clmm2
我对 R 有点陌生,但我想这意味着 tidy 不适用于 clmm2.有人知道解决此问题的方法吗?
I'm a bit new to R, but I guess it means that tidy doesn't work with clmm2. Does anyone know a work-around for this issue?
这是我在上面发布的代码之后的内容,这最终是我想要得到的.
This is what follows the code that I posted above, which is ultimately what I'm trying to get.
然后您可以绘制模拟中的估计分布.
You can then plot the distribution of estimates across your simulations.
ggplot(my_power, aes(estimate, color = term)) +
geom_density() +
facet_wrap(~term, scales = "free")
您也可以将功效计算为小于您的 alpha 的 p 值的比例.
You can also just calculate power as the proportion of p-values less than your alpha.
my_power %>%
group_by(term) %>%
summarise(power <- mean(p.value < 0.05))
推荐答案
对于你需要的,你可以写一个函数来返回具有相同列名的系数:
For what you need, you can write a function to return the coefficients with the same column name:
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")
}
然后我们使用一个例子来应用它,我以序数对葡萄酒数据集进行采样:
Then we apply it using an example where I sample the wine dataset in ordinal:
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")
权力也是如此:
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
这篇关于尝试使用 tidy 进行功率分析并使用 clmm2的文章就介绍到这了,希望我们推荐的答案对大家有所帮助,也希望大家多多支持IT屋!