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

[Rolf Turner]

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

-- 
Honorary Research Fellow
Department of Statistics
University of Auckland
Stats. Dep't. (secretaries) phone:
         +64-9-373-7599 ext. 89622
Home phone: +64-9-480-4619