iid*d and generic0 with Cmatrix

[Gregor Gorjanc]

Hi,

I would like to fit such a model:

y_i = \mu + a1_i + a2_i + e_i

For invidual i:

[a1_i, a2_i]^T ~ N([0, 0]^T, \Sigma_i)

\Sigma_i = \Sigma * A_i

\Sigma[1,1] = \sigma_a1^2

\Sigma[1,2] = \sigma_a12

\Sigma[2,2] = \sigma_a2^2

e_i ~ N(0, \sigma_e^2)

For all i:

[a1, a2]^T ~ N(0, A \kronecker \Sigma)

e ~ N(0, I\sigma_e^2)

I can bring in the generic A matrix via generic0 or I can get Sigma via iid2k (or iidkd), but can I combine/Kronecker them?

Thanks!

Gregor

[Helpdesk (Haavard Rue)]

I would guess so... 

the model formulation is

a1_i + a2_i

and as far as I see, their dependence structure is the same, so am I missing
something?

> --
> You received this message because you are subscribed to the Google Groups "R-
> inla discussion group" group.
> To unsubscribe from this group and stop receiving emails from it, send an
> email to r-inla-discussion...@googlegroups.com.
> To view this discussion, visit
> https://groups.google.com/d/msgid/r-inla-discussion-group/2d158453-0d1d-4ad7-8054-bd7a9ba2d89en%40googlegroups.com
> .

--
Håvard Rue
he...@r-inla.org

[Gregor Gorjanc]

Hi Havard,

In the model  y = Xb + a1 + a2 + e, a1 (nx1) and a2 (nx1) are modelled with MVN(0, A \kronecker Sigma), with A being a known nxn covariance coefficient matrix (we have its inverse) and Sigma being an 2x2 unknown covariance parameter matrix (want to estimate 3 parameters, \sigma_a1^2, \sigma_a2^2, and \sigma_a1,a2).

We can use generic0 for a1 ~ MVN(0, A sigma_a1^2) and a2 ~ MVN(0, A sigma_a2^2).

We can use iid*d for a1, a2 ~ MVN(0, I \kronecker Sigma).

But I don't know how to kronecker-combine the above two into a1, a2 ~ MVN(0, A \kronecker Sigma).

I could push the A matrix into observation model (A = LL^T, a1=Lw1, a2=Lw2, y = Xb + Lw1 + Lw2 + e, and w1, w2 ~ MVN(0, I \kronecker Sigma) - I could use z-model to pass L, but then I again don't know how to kronecker-combine the z-model and iid*d model, and I would loose out on the sparsity in the inverse of A).

Thanks!

gg

[Finn Lindgren]

Maybe I’m missing something, but it doesn’t seem like the Sigma elements are identifiable, since you always use the sum of the pairs, so only

sigma_b^2 = sigma_1^2+sigma_2^2+2sigma_12

would be identifiable, so

a1+a2

Could be replaced by

b ~ N(0, A * sigma_b^2)

which is just a generic0 model.

Finn

On 8 Oct 2025, at 07:00, Gregor Gorjanc <gregor....@gmail.com> wrote:

Hi Havard,

You received this message because you are subscribed to the Google Groups "R-inla discussion group" group.

To unsubscribe from this group and stop receiving emails from it, send an email to r-inla-discussion...@googlegroups.com.

To view this discussion, visit https://groups.google.com/d/msgid/r-inla-discussion-group/239318fb-cf53-40ec-90d7-fcf0f8bd38e3n%40googlegroups.com.

[Finn Lindgren]

That said, your last option, with the Cholesky, can be implemented as an inlabru model with explicit code

solve(Li, a1 + a2)

in the predictor expression, where Li is a sparse matrix from the cholesky of A inverse.

(Not sure if that would give a dense predictor matrix though; probably…)

Finn

On 8 Oct 2025, at 07:36, Finn Lindgren <finn.l...@gmail.com> wrote:

Maybe I’m missing something, but it doesn’t seem like the Sigma elements are identifiable, since you always use the sum of the pairs, so only

[Gregor Gorjanc]

Finn, thanks for pointing out the identifiability issue, which I assume Harvard was also thinking of (saying that coefficient is the same ...). Also thanks for pointing on how to use the solve() with spare Linv in inlabru. Will think more about that angle.

Re identifiability. You are right! I wrote the initial model in haste and tried to boil it down to the simplest case based on two similar models that came up in discussions in our lab, but went too far! Apologies!

This is more appropriate version of model1 (I will just stick to one model to avoid too long e-mail), with the key point that Z1 and Z2 are different and that multiple a2 variables impact one y variable, while just one a1 variable impacts one y variable. So, handling Z1a1 is easy just through indexing, but not for Z2a2:

y = Xb + Z1a1 + Z2a2 + e, a1 (nx1) and a2 (nx1) are modelled with MVN(0, A \kronecker Sigma), with A being a known nxn covariance coefficient matrix (we have its inverse) and Sigma being an 2x2 unknown covariance parameter matrix (want to estimate 3 parameters, \sigma_a1^2, \sigma_a2^2, and \sigma_a1,a2).

We can use generic0 model for a1 ~ MVN(0, A \sigma_a1^2) and a2 ~ MVN(0, A \sigma_a2^2).

We can use iid*d model for a1, a2 ~ MVN(0, I \kronecker Sigma).

We can use z model for Za2 ~ MVN(0, A \sigma_a2^2).

But I don't know how to combine them all into one INLA/inlabru model.

I could push the A matrix into observation model (A = LL^T, a1=Lw1, a2=Lw2, y = Xb + Z1Lw1 + Z2Lw2 + e, and w1, w2 ~ MVN(0, I \kronecker Sigma) - I could use z-model to pass Z1L and Z2L (possibly via inlabru's predictor), but then I again don't know how to kronecker-combine the z-model and iid*d model, and I would loose out on the sparsity in the inverse of A).

Thanks!

gg

[Elias T. Krainski]

As Finn said, it is not identifiable. I've worked an example with the current CRAN version of the graphpcor package and INLAtools, where features to do Kronecker between 'cgeneric' models are implemented (These packages are under development, and some lines in the attached script would be simplified, as per comments therein)

I've set up an identifiable model (one can run by setting simpler = TRUE at the beginning of the script) to check it. 

To view this discussion, visit https://groups.google.com/d/msgid/r-inla-discussion-group/55994789-902F-474E-89F0-4FA617B4E7B6%40gmail.com.

[Finn Lindgren]

The graphpcor approach should work for this, yes!

And with the addition of the different Z1 and Z2 matrices, the model might be identifiable.

For inlabru definition of the predictor expression, there isn’t yet any special handling for iiid2d indexing, but one can always access the latent vector direction (label_latent) and write a general R expression. I plan to add helper features for models that are naturally written as matrix/array instead of long vectors, but for now one needs to do it manually.

Finn

On 8 Oct 2025, at 08:49, Elias T. Krainski <eliask...@gmail.com> wrote:



To view this discussion, visit https://groups.google.com/d/msgid/r-inla-discussion-group/CAA0dk0FMfp-BWehau1qb8G5upOY%3D-LBxYQRBe24UySir_tvx5A%40mail.gmail.com.

<generic0kronecker2d.R>

[Finn Lindgren]

Actually, there _is_ a kind of support for the matrix-vector-product-sum in inlabru, by using a combination of bm_sum() and bm_matrix() for the component mapper.

Probably easier to discuss this offline…

Finn

On 8 Oct 2025, at 08:54, Finn Lindgren <finn.l...@gmail.com> wrote:

The graphpcor approach should work for this, yes!