error returned by MakeADFun to fit LMM with correlated random effects and errors

[Timothée Flutre]

Hello,

I would like to adjust a linear mixed model (LMM) written in the classical way:

y = X beta + Z u + Z_e e

* y is a n-vector

* beta is a p-vector of fixed effects

u is a q-vector of random effects with u ~ N(0, G) where the G matrix is parametrized by a vector gamma containing unknown (co)variances to be estimated

* e is a n-vector of residual errors with e ~ N(0, R) where the R matrix is parametrized by a vector phi containing unknown (co)variances to be estimated.

X, Z and Z_e are design matrices.

It is similar to what glmmTMB can do, but I think glmmTMB cannot deal with correlated random effects and errors. Moreover, importantly, for future use cases, I need to be able to directly pass y as a vector, as well as X, Z and Z_e that I made myself. But for my first use case (which is an example), y equals in fact the vec operator applied to a Y matrix, and G equals the Kronecker product of a small, dense Sigma_u matrix with a known, larger A matrix. Similarly, R equals the Kronecker product of a small, dense Sigma_e with an identity matrix.

For the negative log-likelihood (nll) of my TMB objective function, I add the contrib of U using SEPARABLE, and the contrib of y|u using GMRF. I use these TMB macros to be as efficient as possible. GMRF requires the inverse of the R vcov matrix in a sparse format. I managed to write C++ code to compute it, and it compiles successfully ;)

However, it the call to MakeADFun fails.

I can share my code if anyone is interested to help on this, do not hesitate to ask.

Thanks in advance,

Tim

[Ben Bolker]

"the call to MakeADFun fails" -- it might be helpful to see *how* it
fails. Can you please post the error output? (If it's enormous you could
post just the first few and last few lines ...)

cheers
Ben Bolker

> --
> To post to this group, send email to us...@tmb-project.org. Before
> posting, please check the wiki and issuetracker at https://github.com/
> kaskr/adcomp/ <https://github.com/kaskr/adcomp/>. Please try to create a
> simple repeatable example to go with your question (e.g issues 154, 134,
> 51). Use the issuetracker to report bugs.
> ---
> You received this message because you are subscribed to the Google
> Groups "TMB Users" group.
> To unsubscribe from this group and stop receiving emails from it, send
> an email to tmb-users+...@googlegroups.com <mailto:tmb-
> users+un...@googlegroups.com>.
> To view this discussion visit https://groups.google.com/d/msgid/tmb-
> users/8955038a-9248-4bdd-9ef6-769f73719ca4n%40googlegroups.com <https://
> groups.google.com/d/msgid/tmb-
> users/8955038a-9248-4bdd-9ef6-769f73719ca4n%40googlegroups.com?
> utm_medium=email&utm_source=footer>.

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

[Timothée Flutre]

Sorry, here are the first few lines after calling gdbsource() on a script that simulates data, compiles the C++ code and run it:

Program received signal SIGABRT, Aborted.
__pthread_kill_implementation (threadid=<optimized out>, signo=signo@entry=6, no_tid=no_tid@entry=0)
    at ./nptl/pthread_kill.c:44
#0  __pthread_kill_implementation (threadid=<optimized out>, signo=signo@entry=6,
    no_tid=no_tid@entry=0) at ./nptl/pthread_kill.c:44
#1  0x00007ffff78a9f1f in __pthread_kill_internal (signo=6, threadid=<optimized out>)
    at ./nptl/pthread_kill.c:78
#2  0x00007ffff785afb2 in __GI_raise (sig=sig@entry=6) at ../sysdeps/posix/raise.c:26
#3  0x00007ffff7845472 in __GI_abort () at ./stdlib/abort.c:79
#4  0x00007fffed1e22ed in Eigen::internal::set_from_triplets<__gnu_cxx::__normal_iterator<Eigen::Triplet<double, int>*, std::vector<Eigen::Triplet<double, int>, std::allocator<Eigen::Triplet<double, int> > > >, Eigen::SparseMatrix<double, 0, int>, Eigen::internal::scalar_sum_op<double, double> > (
    begin=..., end=..., mat=..., dup_func=...)
    at /home/tflutre/lib/R/x86_64-pc-linux-gnu-library/4.2/RcppEigen/include/Eigen/src/SparseCore/SparseMatrix.h:1049
#5  0x00007fffed1b7374 in Eigen::SparseMatrix<double, 0, int>::setFromTriplets<__gnu_cxx::__normal_iterator<Eigen::Triplet<double, int>*, std::vector<Eigen::Triplet<double, int>, std::allocator<Eigen::Triplet<double, int> > > > > (this=0x7fffffff1e40, begin=..., end=...)
    at /home/tflutre/lib/R/x86_64-pc-linux-gnu-library/4.2/RcppEigen/include/Eigen/src/SparseCore/SparseMatrix.h:1110
#6  0x00007fffed1e3e67 in initSparsePatFromPackedLU<double> (L=..., Minv=..., nnz=@0x7fffffff1bfc: 0)
    at test_mvLMM.cpp:78
#7  0x00007fffed1b7b01 in makeSparseInvLogDet<double> (M=..., verbose=@0x7fffffff270c: 0, Minv=...,
    logdet=@0x7fffffff1e38: -inf) at test_mvLMM.cpp:146
#8  0x00007fffed19c40f in objective_function<double>::operator() (this=0x7fffffff29b0)
    at test_mvLMM.cpp:243
#9  0x00007fffed1822cb in getParameterOrder (data=0x555561c78978, parameters=0x555561be1168,
    report=0x55556089b118, control=0x55555556fc80)
    at /home/tflutre/lib/R/x86_64-pc-linux-gnu-library/4.2/TMB/include/tmb_core.hpp:2065

I need the sparse inverse of R to give to GMRF, a sparse matrix that I call Rinv_sp. To compute it, I first build R as a sparse matrix, R_sp, and then I compute its sparse inverse Rinv_sp using Eigen's SimplicialLDLT. The error seems to come from the initialization of the sparsity pattern of Rinv_sp.

[Timothée Flutre]

Another way is to run gdbsource("test_mvLMM.R", interactive=TRUE) in a console : it launches test_mvLMM.R that compile test_mvLMM.cpp (attached). I am using R 4.2.2.

[Kasper Kristensen]

You are almost there...

From the interactive gdb session you can now inspect the frames that look suspicious. In this case my first guess was to look at 'setFromTriplets'' (frame 5).

I then printed the 'triplets' object:

(gdb) p triplets

$2 = std::vector of length 724, capacity 724 = {{m_row = 0, m_col = 0, m_value = 1}, {m_row = 0, m_col = 0, m_value = 1}, {m_row = 0, m_col = 0,
    m_value = 1}, {m_row = 0, m_col = 207145260, m_value = 4.6344408982361827e-310}, {m_row = 1, m_col = 1, m_value = 1}, {m_row = 1, m_col = 1,
    m_value = 1}, {m_row = 1, m_col = 1, m_value = 1}, {m_row = 2, m_col = 2, m_value = 1}, {m_row = 2, m_col = 2, m_value = 1}, {m_row = 2,
    m_col = 2, m_value = 1}, {m_row = 2, m_col = 21840, m_value = 0}, {m_row = 3, m_col = 3, m_value = 1}, {m_row = 3, m_col = 3, m_value = 1}

There's definitely something wrong with this element:

 m_row = 0, m_col = 207145260, m_value = 4.6344408982361827e-310}

Hopefully that gives you a clue on what to fix.

[Timothée Flutre]

Thank you! I will follow this direction (after updating my system because for me "p triplets" does not return anything, unfortunately).

To post to this group, send email to us...@tmb-project.org. Before posting, please check the wiki and issuetracker at https://github.com/kaskr/adcomp/. Please try to create a simple repeatable example to go with your question (e.g issues 154, 134, 51). Use the issuetracker to report bugs.

---
You received this message because you are subscribed to the Google Groups "TMB Users" group.

To unsubscribe from this group and stop receiving emails from it, send an email to tmb-users+...@googlegroups.com.
To view this discussion visit https://groups.google.com/d/msgid/tmb-users/d6452fce-6395-4655-90f9-72a5f9c51a5cn%40googlegroups.com.