glmmTMB syntax to brm() syntax
On 2024-10-28 7:43 p.m., Simon Harmel wrote:
All good now, I managed to get the 1.1.10 version installed. So, with
large enough data (ex. below), I should be able to recover the
parameters, right?
I ask because I can recover all parameters, but wonder how to
determine if my model is recovering theta=-1?
set.seed(101)
dd <- data.frame(ID = rep(1:10000, each=2),
con = rep(c("simple","complex"), 10000))
dd$VOCD <- simulate_new(
~ 0+con + (1|ID),
dispformula = ~ 0+con,
newdata = dd,
newparams = list(beta = c(0, .5), theta = -1,
betadisp = c(2,1.6)))[[1]]
glmmTMB(VOCD ~ 0+con + (1|ID),
dispformula = ~ 0+con, data = dd)
Saving that last model as g1,
log(sqrt(c(VarCorr(g1)$cond[[1]])))
or if you like tidyverse
library(tidyverse)
broom.mixed::tidy(g1, effects = "ran_pars")
|> filter(term == "sd__(Intercept)")
|> pull(estimate)
|> log()
These both give -1.155614, which seems reasonably close.
On Mon, Oct 28, 2024 at 1:33?PM Ben Bolker <bbolker at gmail.com> wrote:
You should be getting version 1.1.10 from CRAN, unless you're using
an outdated mirror, see
https://cran.r-project.org/web/packages/glmmTMB/index.html ... you will
need version 1.1.10 to avoid the error.
If you want a dispersion parameter of 1.0 then you probably want
betadisp = 0 (since the dispersion model uses a log link).
You also need to provide only a vector of length 1 for theta, since
you only have a single (scalar) random effect term. And you need to
provide two values for betadisp (since you are specifying ~con, and con
is a two-level factor), i.e.
dd$VOCD <- simulate_new(
~ 0+con + (1|ID),
dispformula = ~ con,
newdata = dd,
newparams = list(beta = c(0, 0.5), theta = -1,
betadisp = rep(0,2)))[[1]]
On 2024-10-28 1:39 p.m., Simon Harmel wrote:
Thank you, Ben. I deleted and installed glmmTMB from CRAN and it is
still version 1.1.9 on my machine with that same error showing up.
I love the 'covstruct' vignette. But I, for example, don't know what
to use in my simulate_new() call below to define the dispersion,
would that be 'sigmadisp'?
dd <- data.frame(ID = rep(1:100, each = 2),
con = rep(c("simple","complex"), 100)) %>%
group_by(ID) %>% mutate(obs=row_number())
dd$VOCD <- simulate_new(
~ 0+con + (1|ID),
dispformula = ~ con,
newdata = dd,
newparams = list(beta = c(0, 0.5), theta = rep(-1,2),
disp = 1))[[1]]
On Mon, Oct 28, 2024 at 12:10?PM Ben Bolker <bbolker at gmail.com> wrote:
That's quite old, by R standards (1.1.8 came out a year ago, latest
version is 1.1.10 -- there's a note in the NEWS for 1.1.10
<https://cran.r-project.org/web/packages/glmmTMB/news.html> about fixing
a simulate bug for the beta family
(It's simulate_new(), not new_simulate())
The docs say "(?beta?, ?betazi?, ?betadisp?, ?theta?, etc.)"; in this
case "etc." includes thetazi (random-effects parameters for ZI terms, if
any) and thetadisp (ditto, dispersion terms).
The 'covstruct' vignette is the best place to read about
parameterization of random-effects components.
On 2024-10-28 1:05 p.m., Simon Harmel wrote:
Ben, when I run your code below, I get the following error message: y values must be 0 <= y < 1. My glmmTMB version is: ?1.1.7?. Also, is there any good documentation explaining all possible names used in argument `newparams=` in glmmTMB::new_simulate()? Thank you! Simon On Thu, Oct 24, 2024 at 11:11?AM Ben Bolker <bbolker at gmail.com> wrote:
See below. The two models (glmmTMB and brms) give sufficiently
similar estimates that I'm confident that the specifications match.
set.seed(101)
library(glmmTMB)
library(brms)
library(broom.mixed)
library(tidyverse)
dd <- data.frame(ID = rep(1:100, each = 10),
TRIAL_INDEX = rep(1:10, 100),
con = rnorm(1000))
dd$pic_percent <- simulate_new(
~ con + (0+con | ID) +
(0+con | TRIAL_INDEX),
ziformula = ~1,
family = beta_family(),
newdata = dd,
newparams = list(beta = c(0, 0.5), theta = rep(-1,2),
betadisp = 1, betazi = -2))[[1]]
m1 <- glmmTMB(pic_percent ~ con + (0+con | ID) +
(0+con | TRIAL_INDEX),
data=dd,
family = beta_family(),
ziformula = ~1)
##
https://mvuorre.github.io/posts/2019-02-18-analyze-analog-scale-ratings-with-zero-one-inflated-beta-models/
m2 <- brm(
bf(pic_percent ~ con + (0+con | ID) +
(0+con | TRIAL_INDEX),
zi = ~ 1),
data=dd,
family = zero_inflated_beta()
)
(purrr::map_dfr(list(glmmTMB = m1, brms = m2), tidy, .id = "model")
|> select(model, effect, component, group, term, estimate)
|> pivot_wider(names_from = model, values_from = estimate)
)
On 10/23/24 19:13, Simon Harmel wrote:
Hello all,
I was wondering what is the closest equivalent of my glmmTMB syntax below
in brms::brm() syntax?
glmmTMBglmmTMB(pic_percent ~ con +
(0+con | ID) +
(0+con | TRIAL_INDEX),
data=DATA,
family = beta_family(),
ziformula = ~1)
Thank you,
Simon
[[alternative HTML version deleted]]
_______________________________________________ R-sig-mixed-models at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
-- Dr. Benjamin Bolker Professor, Mathematics & Statistics and Biology, McMaster University Director, School of Computational Science and Engineering * E-mail is sent at my convenience; I don't expect replies outside of working hours.
_______________________________________________ R-sig-mixed-models at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
-- Dr. Benjamin Bolker Professor, Mathematics & Statistics and Biology, McMaster University Director, School of Computational Science and Engineering
> E-mail is sent at my convenience; I don't expect replies outside of
working hours.
-- Dr. Benjamin Bolker Professor, Mathematics & Statistics and Biology, McMaster University Director, School of Computational Science and Engineering
> E-mail is sent at my convenience; I don't expect replies outside of
working hours.
Dr. Benjamin Bolker Professor, Mathematics & Statistics and Biology, McMaster University Director, School of Computational Science and Engineering > E-mail is sent at my convenience; I don't expect replies outside of working hours.