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

[Ben Bolker]

   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.