[R-sig-ME] glmmTMB syntax to brm() syntax

[Simon Harmel]

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)

On Mon, Oct 28, 2024 at 1:33 PM Ben Bolker <bbolker using 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 using 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 using 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 using 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 using 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.
>