5

I have a data set (df) in this format

index <- runif(n = 100,min = 0, max = 1)
type1 <- rep("low", 50)
type2 <- rep("high", 50)
type <- c(type1,type2)

level1 <- rep("single", 25)
level2 <- rep("multiple", 25)
level3 <- rep("single", 25)
level4 <- rep("multiple", 25)
level <- c(level1,level2,level3,level4)

block <- rep(1:5, 10)
set <- rep(1:5, 10)

df <- data.frame("index" = index,"type" = type, "level" = level, "block" = block, "set" = set)
df$block <- as.factor(df$block)
df$set <- as.factor(df$set)

I want to create a model that looks like like this

model <- lmer(index ~ type * level + (1|block) + (1|set), data = df)

However, in my original data the fit is bad because the data is bound between 0 and 1. I want to bootstrap this mixed effects model. Any idea on how to achieve boot-strapping for such a model? I want to compare this this full model with sub-models eg. without interaction, or with level or type alone. I also want with confidence intervals for the final model

Rspacer
  • 2,369
  • 1
  • 14
  • 40
  • A beta regression would be ideal, as it's specifically suited to outcomes bound between 0 and 1. – Phil Jan 23 '19 at 03:55
  • I tried `glmmTMB(index ~ type + level + (1|block), df, family=list(family="beta",link="logit"))` betareg but it gives an infinite loop warning output saying "singular fit" – Rspacer Jan 23 '19 at 04:08
  • I tried `brms::brm(index ~ type * level + (1|block) + (1|set), data = df, family = Beta())` and it seems to have done fine with it. Just a question of whether you're ok running a Bayesian model. – Phil Jan 23 '19 at 04:14
  • I have never worked with Bayesian models before. How do I read the output? How do I compare to see if type has a significant effect or not? Like in GLMMs do I have to do model reduction with and without the fixed effects, and run a anova? – Rspacer Jan 23 '19 at 13:25
  • The estimates have the same meaning in both bayesian and frequentist models, and can be taken using `summary(model)`. The intervals simply indicate the (default 95%) probability of the "true" coefficient being in that interval. So if the interval for the type coefficient includes 0, you can consider it non-significant. As for model comparison, I prefer to use information criteria to compare, such as WAIC or PSIS-LOO. The `brms` package has the `waic()`, `loo()`, and `compare_ic()` functions to that allow you to compare the models' respective WAIC or LOO scores. – Phil Jan 23 '19 at 17:02
  • @Biotechgeek you could also try `GLMMadaptive::mixed_model()` by Dimitris Rizopoulos. He also has great vignettes on GitHub. – Dion Groothof Jan 20 '22 at 14:54

1 Answers1

0

The confint() function has a method for merMod objects. The following should work:

confint(model, method = "boot", nsim = 1000)

And with multiple CPUs:

confint(model, method = "boot", nsim = 1000,
        parallel = "multicore", ncpus = 8)
Brigadeiro
  • 2,649
  • 13
  • 30
  • This produces an infinite loop saying "singular fit". Also, I have modified my question to reflect what I'm looking for – Rspacer Jan 23 '19 at 03:16