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

[Ben Bolker]

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

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