diff --git a/README.md b/README.md
index bf55c2f434f0d31541b83aa14755b089cef51416..c199afbc3dc2a115b6155283b2cb24b260753f4d 100644
--- a/README.md
+++ b/README.md
@@ -20,11 +20,12 @@ It might also be used to distribute development versions of `spaMM`. However, us
 
 The `spaMM` package was developed first to fit mixed-effect models with spatial correlations, which commonly occur in ecology, but it has since been developed into a more general package for inferences under models with or without spatially-correlated random effects, including multivariate-response models. To make it competitive to fit large data sets, `spaMM` has distinct algorithms for three cases: sparse precision, sparse correlation, and dense correlation matrices, and is efficient to fit geostatistical, autoregressive, and other mixed models on large data sets. Notable features include: 
 
-- Fitting spatial and non-spatial correlation models: **geostatistical** models with random-effect terms following the `Matern` as well as the much less known `Cauchy` correlation models, **autoregressive** models described by an `adjacency` matrix, AR(_p_) and ARMA(_p_,_q_) time-series models (`ARp` and `ARMA`), or an **arbitrary given** precision or correlation matrix (`corrMatrix`). Conditional spatial effects can be fitted,  as in (say) `Matern(female|...) + Matern(male|...)` to fit distinct random effects for females and males (e.g., [Tonnabel et al., 2021](https://doi.org/10.1111/mec.15833)). Brave users can even define their own parametric correlation models, to be fitted as any other random effect (the `corrFamily` feature).
+- Fitting spatial and non-spatial correlation models: **geostatistical** models with random-effect terms following the `Matern` as well as the much less known `Cauchy` correlation models, **autoregressive** models described by an `adjacency` matrix, AR(_p_) and ARMA(_p_,_q_) time-series models (`ARp` and `ARMA`), or an **arbitrary given** precision or correlation matrix (`corrMatrix`). Conditional spatial effects can be fitted,  as in (say) `Matern(female|...) + Matern(male|...)` to fit distinct random effects for females and males (e.g., [Tonnabel et al., 2021](https://doi.org/10.1111/mec.15833)). "Composite" random effects that combine features of such autocorrelated random effects and of random-coefficient models (say, `Matern(age|...)`) can be fitted. Brave users can even define their own parametric correlation models, to be fitted as any other random effect (the `corrFamily` feature).
 - A further class of spatial correlation models, "Interpolated Markov Random Fields" (`IMRF`) covers widely publicized approximations of Matérn models ([Lindgren et al. 2011](http://doi.org/10.1111/j.1467-9868.2011.00777.x)) and the multiresolution model of [Nychka et al. 2015](https://doi.org/10.1080/10618600.2014.914946). 
 - Symmetric and antisymmetric **dyadic interaction** effects (such as considered in so-called Bradley-Terry models or in diallel experiments) can be fitted as fixed or as random effects (see e.g. `X.antisym`, `diallel` or `antisym` documentations)  
 - Allowed response families include beta response, beta-binomial, the Conway-Maxwell-Poisson (`COMPoisson`), and two negative binomial families. Zero-truncated variants of the `poisson` and negative-binomial families are handled;
-- All the above features combined in multivariate-response models;
+- All the above features combined in multivariate-response models, including random effects correlated over different response variables;
+- ML and REML fits (see below for comments on likelihood approximations);
 - A replacement function for `glm`, useful when the latter (or even `glm2`) fails to fit a model;
 - A syntax close to that of `glm` or [`g`]`lmer`. 
 - Many extractor methods similar to those in `stats` or `nlme`/`lmer`, and functions for inference beyond the fits, such as `confint()` for confidence intervals of fixed-effect parameters, `predict()` and related functions for point prediction and prediction variances, and compatibility with functions from other packages such as `multcomp::glht()` and `lmerTest` procedures providing F tests using Satterthwaite method (see `` `post-fit` `` and `anova` documentation items); 
@@ -57,9 +58,9 @@ The performance of likelihood ratio tests based on `spaMM` fits, and the impact
     Rousset F., Ferdy J.-B. (2014) [Testing environmental and genetic effects in the presence of spatial autocorrelation](http://onlinelibrary.wiley.com/doi/10.1111/ecog.00566/abstract). Ecography, 37: 781-790.
 Also available here is the [Supplementary Appendix G](http://kimura.univ-montp2.fr/%7Erousset/spaMM/RoussetF14AppendixG.pdf) from that paper, including comparisons with a trick that has been uncritically used to constrain the functions `lmer` and `glmmPQL` to analyse spatial models.
 
-For some substantial use of various features of `spaMM`, see e.g. the [IsoriX project](https://github.com/courtiol/IsoriX), or a story about [social dominance in hyaenas](https://doi.org/10.1038/s41559-018-0718-9), or [yet another depressing story about climate change](https://doi.org/10.1038/s41467-019-10924-4), or [the life-history of mothers of twins](https://doi.org/10.1038/s41467-022-30366-9), or a comparison of prediction by LMMs and by random-forest methods (in supplementary material of [a paper on protected area personnel](https://doi.org/10.1038/s41893-022-00970-0)), or [analyses of dyadic interactions in mandrills](https://doi.org/10.7554/eLife.79417). 
+For some substantial use of various features of `spaMM`, see e.g. the [IsoriX project](https://github.com/courtiol/IsoriX), or a story about [social dominance in hyaenas](https://doi.org/10.1038/s41559-018-0718-9), or [yet another depressing story about climate change](https://doi.org/10.1038/s41467-019-10924-4), or [the life-history of mothers of twins](https://doi.org/10.1038/s41467-022-30366-9), or a comparison of prediction by LMMs and by random-forest methods (in supplementary material of [a paper on protected area personnel](https://doi.org/10.1038/s41893-022-00970-0)), or [analyses of dyadic interactions in mandrills](https://doi.org/10.7554/eLife.79417), or [comparative analyses by heteroscedastic models with phylogenetic autocorrelation](https://doi.org/10.1038/s42003-024-06165-x). 
 
-Initial development drew inspiration from work by Lee and Nelder on _h_-likelihood (e.g. [Lee, Nelder & Pawitan](https://doi.org/10.1201/9781420011340), 2006; [Lee & Lee](http://dx.doi.org/10.1007/s11222-011-9265-9) 2012; see also [Molas and Lesaffre](http://dx.doi.org/10.1002/sim.3852), 2010), and `spaMM` retains from that work several distinctive features, such as specific methods to fit models with non-gaussian random effects, structured dispersion models with random effects, and implementation of several variants of Laplace and PQL approximations. However, later versions have increasingly relied on additional insights. Notably, the default likelihood approximation now goes beyond those discussed in these works, and is the same Laplace approximation as in `TMB` ([Kristensen et al., 2016](https://doi.org/10.18637/jss.v070.i05)) and packages based on `TMB`, in particular where this departs from what is discussed in the _h_-likelihood literature (i.e., for GLM families with non-canonical link, or response families not of the GLM class). 
+Initial development drew inspiration from work by Lee and Nelder on _h_-likelihood and more elaborate approximations of likelihood (e.g. [Lee, Nelder & Pawitan](https://doi.org/10.1201/9781420011340), 2006; [Lee & Lee](http://dx.doi.org/10.1007/s11222-011-9265-9) 2012; see also [Molas and Lesaffre](http://dx.doi.org/10.1002/sim.3852), 2010). The latter two references, and `spaMM` itself, may actually be more widely understood as applications of Laplace approximations of likelihood than as applications of a distinctive _h_-likelihood concept. `spaMM` retains from Lee & Nelder's work several distinctive features, such as a concept of restricted likelihood applicable beyond LMMs, specific methods to fit models with non-gaussian random effects, structured dispersion models with random effects, and implementation of several variants of Laplace and PQL approximations. However, it has departed from their work in various ways. Notably, the default likelihood and restricted likelihood approximations now go beyond those discussed in these works. For ML fits, it is the same Laplace approximation as in `TMB` ([Kristensen et al., 2016](https://doi.org/10.18637/jss.v070.i05)) and packages based on `TMB`, because `TMB` and `spaMM` (with default arguments) use the observed Hessian matrix of ``joint likelihood'' where the _h_-likelihood literature only considers the expected Hessian matrix. This makes a difference for GLM families with non-canonical link, or for response families not of the GLM class. 
 
 ## Credits
 Initial development was supported by a PEPS grant from the CNRS and University of Montpellier.