[R-sig-ME] Fixing non-symetrix sigma matrix

[Simon Harmel]

Thank you so much all for your responses.
Ben, I wonder what would be the equivalent of my mediator model in brms?

mediator <- glmmTMB(pic_percent ~ con +
                      (0+con | ID) +
                      (0+con | TRIAL_INDEX),
                    data=df,
        family = beta_family(),
        ziformula = ~1)

On Mon, Oct 21, 2024, 9:28 PM Ben Bolker <bbolker using gmail.com> wrote:

>
>    The mediate() function simply isn't set up to work with glmmTMB
> objects.  Rolf's fix takes care of the fact that vcov() applied to a
> glmmTMB object returns a list of covariance matrices (for the fixed
> effect parameters of the conditional distribution, zero-inflation model,
> dispersion model) rather than a single matrix.
>
>    Beyond that, though, if you look at the code for mediate() you'll see
> a whole bunch of logic that checks whether each of the two models
> inherits from classes `glm`, `gam`, `survreg`, `merMod`, etc. ... by
> structuring their code this way (rather than using S3 methods), the
> authors have made it hard to adapt their code to new model types without
> digging all the way into the code and writing new components for the new
> model type ...
>
>    You could ask the package maintainers to add these capabilities, or
> try to find someone else who can do it for you ... (it might not be too
> hard, but the mediate() function is 1400+ lines long ... the way to go
> at it would be to look for everywhere that `isMer.m` or `isMer.y` (main
> or mediation model is from lme4) and write parallel code that works with
> glmmTMB objects instead ...)
>
>    Ben Bolker
>
>
>
> On 10/21/24 18:57, Rolf Turner wrote:
> >
> > On Mon, 21 Oct 2024 13:34:56 -0500
> > Simon Harmel <sim.harmel using gmail.com> wrote:
> >
> >> Hello Mixed-Models experts,
> >>
> >> I'm trying to fit the following reproducible mediation model called
> >> final. But I get an error saying: ...sigma must be a symmetric
> >> matrix...
> >>
> >> Could you please advise how I can possibly fix this error?
> >>
> >> Many thanks,
> >> Simon
> >>
> >> Reproducible R code:
> >> df <-
> >> read.csv("
> https://raw.githubusercontent.com/fpqq/w/refs/heads/main/t.csv")
> >>
> >> library(glmmTMB)
> >> library(mediation)
> >>
> >> mediator <- glmmTMB(pic_percent ~ con +
> >>                        (0+con | ID) +
> >>                        (0+con | TRIAL_INDEX),
> >>                      data=df,
> >>          family = beta_family(),
> >>          ziformula = ~1)
> >>
> >> outcome <- glmmTMB(acc ~ con + pic_percent +
> >>                        (0+con+pic_percent | Q),
> >>                       data = df,
> >>                     family = binomial())
> >>
> >> final <- mediate(mediator, outcome, sims=50,
> >>                      treat="con", mediator="pic_percent")
> >>
> >
> > Thank you for the clearly structured reproducible example.
> >
> > I must preface my remarks by saying that this is definitely a case of
> > the blind leading the partially-sighted; I am no "mixed model expert",
> > and I really haven't a clue about "mediation".
> >
> > Having said that, there seems to me (in my ignorance) to be several
> > things wrong in the code of mediate().  If we do
> >
> >      sink("mediate.R")
> >      mediate
> >      sink()
> >
> > we find at line 469 of mediate.R the code:
> >
> >      MModel <- rmvnorm(sims, mean = MModel.coef, sigma = MModel.var.cov)
> >
> > This seems to be incorrect.  MModel.var.coef should be a 2 x 2 matrix,
> > but it is a list.  The first entry of this list seems to be the
> > required covariance matrix.  Similarly MModel.coef is a list; the
> > second two entries seem to be empty lists.  The first entry is a list of
> > length 2, each entry being a matrix.  It is not at all clear to me what
> > role these matrices play.  If one modifies the call to rmvnorm() to,
> > e.g.
> >
> >      MModel <- rmvnorm(sims, mean = MModel.coef[[1]][[1]],
> >                        sigma = MModel.var.cov[[1]])
> >
> > then the code runs a bit further, but rmvnorm() hangs up at the end
> > with an error.  The call to sweep()
> >
> >      retval <- sweep(retval, 2, mean, "+")
> >
> > complains about a dimension miss-match, but then throws an error
> > "non-numeric argument to binary operator" which mystifies me completely.
> >
> > At this point I give up, and hope that someone who knows what they are
> > doing will take over.
> >
> > cheers,
> >
> > Rolf Turner
> >
>
> --
> 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
>

	[[alternative HTML version deleted]]