Trying to use tidy for a power analysis and using clmm2

Question

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)

Now, I have the following function:

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)

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:

Error: No tidy method for objects of class 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?

EDIT: 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")

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

you can pull out the results with coefficients(summary()), convert to a dataframe with the rownames as one of the columns? — StupidWolf, Feb 13 '20 at 11:53
Probably? I'm following someone else code - I edited my post to show what would come next. As I said, I'm a bit new to R, sorry if my reply sounds dumb. — GWasp, Feb 13 '20 at 12:07

score 4 · Accepted Answer · answered Feb 13 '20 at 14:22

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)

Plotting works:

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

And so does power:

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

Hi @BenBolker, you mean the example or code? Sure thing for both.. — StupidWolf, Feb 13 '20 at 15:12
huh, that's weird. The ordinal (`clmm`) tidiers are in broom, not broom.mixed ... — Ben Bolker, Feb 13 '20 at 15:19

Trying to use tidy for a power analysis and using clmm2

1 Answers1

Linked