gam with MRF penalty crashes R session

[mluerig]

I'm running a series of gam models using mgcv, and one of them just silently crashes the current R session. Just to give you an idea of the structure (three different subsets run fine):

mod3p_ld1_f_v = gam(prob ~ 
                      s(LD1, k=5) +
                      s(species, bs="mrf", xt=list(penalty=phylo_penalty)), 
                    family = betar(link = "logit"),
                    method="REML",
                    data=data_anc_fit[subset=="ld1_f_v"]
)

I'm not sure how I would debug this, since there is no error message. I'm happy to share the data though for you to reproduce what happened.

> sessionInfo()
R version 4.4.2 (2024-10-31 ucrt)
Platform: x86_64-w64-mingw32/x64
Running under: Windows 11 x64 (build 26100)

Matrix products: default


locale:
[1] LC_COLLATE=English_Germany.utf8  LC_CTYPE=English_Germany.utf8    LC_MONETARY=English_Germany.utf8 LC_NUMERIC=C                     LC_TIME=English_Germany.utf8    

time zone: America/New_York
tzcode source: internal

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] mvMORPH_1.2.1       subplex_1.9         corpcor_1.6.10      ggtree_3.14.0       phytools_2.4-4      maps_3.4.2.1        ape_5.8-1           lmerTest_3.1-3      lme4_1.1-37         Matrix_1.7-1        sjPlot_2.8.17      
[12] MRFtools_0.0-2      caret_7.0-1         lattice_0.22-6      AICcmodavg_2.3-4    factoextra_1.0.7    FactoMineR_2.11     gratia_0.10.0       mgcv_1.9-1          nlme_3.1-166        scales_1.3.0        corrplot_0.95      
[23] egg_0.4.5           gridExtra_2.3       ggimage_0.3.3       plotly_4.10.4       pals_1.10           ggnewscale_0.5.1    ggridges_0.5.6      cowplot_1.1.3       googlesheets4_1.1.1 stringi_1.8.7       moments_0.14.1     
[34] reshape_0.8.9       lubridate_1.9.4     forcats_1.0.0       stringr_1.5.1       dplyr_1.1.4         purrr_1.0.4         readr_2.1.5         tidyr_1.3.1         tibble_3.2.1        ggplot2_3.5.2       tidyverse_2.0.0    
[45] MASS_7.3-61         terra_1.8-42        data.table_1.17.0   pacman_0.5.1       

loaded via a namespace (and not attached):
  [1] splines_4.4.2           ggplotify_0.1.2         cellranger_1.1.0        datawizard_1.0.2        hardhat_1.4.1           pROC_1.18.5             rpart_4.1.23            lifecycle_1.0.4         Rdpack_2.6.4           
 [10] sf_1.0-20               pbmcapply_1.5.1         doParallel_1.0.17       globals_0.16.3          insight_1.1.0           flashClust_1.01-2       magrittr_2.0.3          remotes_2.5.0           spam_2.11-1            
 [19] mapproj_1.2.11          DBI_1.2.3               minqa_1.2.8             multcomp_1.4-28         quadprog_1.5-8          expm_1.0-0              yulab.utils_0.2.0       nnet_7.3-19             TH.data_1.1-3          
 [28] sandwich_3.1-1          ipred_0.9-15            lava_1.8.1              ggrepel_0.9.6           listenv_0.9.1           tidytree_0.4.6          units_0.8-7             performance_0.13.0      parallelly_1.43.0      
 [37] codetools_0.2-20        DT_0.33                 tidyselect_1.2.1        aplot_0.2.5             ggeffects_2.2.1         farver_2.1.2            stats4_4.4.2            googledrive_2.1.1       jsonlite_2.0.0         
 [46] glassoFast_1.0.1        e1071_1.7-16            survival_3.7-0          iterators_1.0.14        emmeans_1.11.0          foreach_1.5.2           tools_4.4.2             treeio_1.30.0           Rcpp_1.0.14            
 [55] glue_1.8.0              mnormt_2.1.1            prodlim_2024.06.25      xfun_0.52               mvnfast_0.2.8           numDeriv_2016.8-1.1     withr_3.0.2             combinat_0.0-8          fastmap_1.2.0          
 [64] boot_1.3-31             digest_0.6.37           timechange_0.3.0        R6_2.6.1                gridGraphics_0.5-1      estimability_1.5.1      colorspace_2.1-1        dichromat_2.0-0.1       generics_0.1.3         
 [73] recipes_1.2.1           clusterGeneration_1.3.8 class_7.3-22            httr_1.4.7              htmlwidgets_1.6.4       scatterplot3d_0.3-44    ModelMetrics_1.2.2.2    pkgconfig_2.0.3         gtable_0.3.6           
 [82] timeDate_4041.110       htmltools_0.5.8.1       dotCall64_1.2           multcompView_0.1-10     leaps_3.2               gower_1.0.2             reformulas_0.4.0        ggfun_0.1.8             knitr_1.50             
 [91] rstudioapi_0.17.1       tzdb_0.5.0              reshape2_1.4.4          coda_0.19-4.1           curl_6.2.2              nloptr_2.2.1            proxy_0.4-27            zoo_1.8-14              DEoptim_2.2-8          
[100] sjlabelled_1.2.0        KernSmooth_2.23-24      parallel_4.4.2          pillar_1.10.2           grid_4.4.2              vctrs_0.6.5             VGAM_1.1-13             xtable_1.8-4            cluster_2.1.6          
[109] evaluate_1.0.3          magick_2.8.6            mvtnorm_1.3-3           cli_3.6.4               compiler_4.4.2          rlang_1.1.6             future.apply_1.11.3     classInt_0.4-11         plyr_1.8.9             
[118] sjmisc_2.8.10           fs_1.6.6                viridisLite_0.4.2       munsell_0.5.1           optimParallel_1.0-2     lazyeval_0.2.2          sjstats_0.19.0          hms_1.1.3               patchwork_1.3.0        
[127] unmarked_1.5.0          future_1.40.0           rbibutils_2.3           igraph_2.1.4            gargle_1.5.2            phangorn_2.12.1         fastmatch_1.1-6         ggokabeito_0.1.0   

[gavinsimpson]

With apologies to @eric-pedersen, GH deleted both comments

Thanks for the report. I can't easily debug the issue without being able to see the data. Would you be willing to post a minimum reproducible sample with the phylogeny? Simulated data is fine, but it would need to have enough similarities to replicate the crash itself.

A couple things I would suggest checking:

• How large is your species penalty matrix (number of rows and columns)? This might be a memory issue due to trying to solve too large a matrix. If so, using bam instead of gam
• Does it still crash R when using a different family (e.g. gaussian()?
• Does it crash when the phylogenetic smoother is excluded? What about the other smooth term?
• Are all the species used to create the penalty available as levels in the subset of data used for the model? If not, does setting dropUnusedLevels=TRUE in the gam call solve the issue?

[mluerig]

How large is your species penalty matrix (number of rows and columns)? ->324 species
Does it still crash R when using a different family (e.g. gaussian()? -> gaussian works!
Does it crash when the phylogenetic smoother is excluded? What about the other smooth term? -> doesnt crash without
Are all the species used to create the penalty available as levels in the subset of data used for the model? -> all species are present

This is the actual data: data.csv
and this is the tree (rename to .tre) tree.txt

tree = read.tree(file="tree.tre")
data = fread(file="data.csv")

phylo_penalty <- mrf_penalty(tree)
data[, species := factor(species, levels=tree$tip.label)]

## this works
mod3p_ld1_f_v = gam(prob ~ 
                      s(LD1, k=5) +
                      s(species, bs="mrf", xt=list(penalty=phylo_penalty)), 
                    family = betar(link = "logit"),
                    method="REML",
                    data=data[subset=="ld1_f_d"]
)

## this crashes
mod3p_ld1_f_v = gam(prob ~ 
                      s(LD1, k=5) +
                      s(species, bs="mrf", xt=list(penalty=phylo_penalty)), 
                    family = betar(link = "logit"),
                    method="REML",
                    data=data[subset=="ld1_f_v"]
)

EDIT: so gaussian works, but it gives really awkward results for this and other traits which are on a scale between 0 and 1. betar seemed ideal...

[eric-pedersen]

Hi @mluerig,

I've re-created the issue on my end, but I am not 100% sure why it's crashing. I am pretty sure the issue is not with MRFtools per-se, though, but instead is some sort of memory overflow issue within gam.

Here's my general observations:

I see that each subset of the data you are fitting your model on has all 324 species in your data set, each one observed once, and you are using the default k for the mrf smoother, which in this case would be k=324 (one for each species). That means that each model actually has more coefficients than there are observations (1 intercept, 4 basis functions for LD1 (after dropping the null basis function), and 323 basis functions for the mrf smoother). While mgcv no longer returns an error when you have more parameters than data points (as it is possible to estimate this kind of model when at least part of the parameter space is penalized), I suspect that this is probably the cause of the issue, as this is not the kind of model that mgcv is really designed to estimate.

However, I found I can avoid the crash occurring by simply reducing the number of basis functions given to the mrf smoother such that the total number of coefficients is not larger than the data. This uses a low-rank approximation to the full MRF smoother, where the MRF smoother is approximated by a lower-dimensional smoother (see Simon Wood's 2017 book, section 5.8.1 for details on this approximation). What gets dropped in this approximation is the components of the penalty corresponding to the highest frequency differences (generally, fine-scale differences in values between groups of species very near the tips), which if you only have 324 observations per model, are the hardest basis functions to estimate accurately anyway.

Here's an example of this model, where I have reduced the phylogenic smoother to have 300 basis functions:

#this model shouldn't crash
mod3p_ld1_f_v = gam(prob ~ 
                      s(LD1, k=5) +
                      s(species, bs="mrf", xt=list(penalty=phylo_penalty),
                        k = 300), 
                    family = betar(link = "logit"),
                    method="REML",
                    data=data[subset=="ld1_f_v"])

summary(mod3p_ld1_f_v)

And this is the model output:

Family: Beta regression(26.195) 
Link function: logit 

Formula:
prob ~ s(LD1, k = 5) + s(species, bs = "mrf", xt = list(penalty = phylo_penalty), 
    k = 300)

Parametric coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -0.14200    0.02354  -6.032 1.62e-09 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
               edf  Ref.df Chi.sq p-value    
s(LD1)       3.966   3.986  227.3  <2e-16 ***
s(species) 172.207 299.000  758.6  <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =  0.881   Deviance explained = 94.1%
-REML = -178.98  Scale est. = 1         n = 324

Note that the phylogenic smoother only has an EDF of 172, indicating that using 300 basis functions might have been overkill anyway.

[mluerig]

Hey @eric-pedersen,

thanks for the fast response - with your help I got it to work on my end as well. I used k=N/2 and I get plausible results.

I checked the section from Simon Wood's book and ran the example. So, if I understand correctly, with the lower-dimensional mrf-smoother, I am approximating more coarse phylogenetic relationships (geographic in the example) rather than modelling all species at the tips explicitly, correct? Why then is the default a full rank smoother with all species - wouldn't this always give us a case where the model has more coefficients than there are observations?

[gavinsimpson]

@mluerig If you have $m$ species, the MRF basis will have $m - 1$ "basis functions" because the model will also have an intercept. I can't think of a situation where you'd have fewer observations $n$ (note I don't mean sampling locations here, I mean rows of data) than species because even if all you have is 1 observation per sampling location ($s$) you'd be stacking the $m$ vectors of species data into a vector of length $m \times s = n$, and hence $m \ll n$.

If we think about the spatial setting, where the MRF is most often encountered, if we have 1 observation per each of $l$ spatial locations (so we're in a purely spatial smoothing setup) and we use a full rank MRF basis, that basis would be of dimension $l - 1$ due to the identifiability constraint and hence the model would have as many coefficients as data. Assuming no other model terms are included (beyond the intercept / constant term) this model is identifiable because the MRF is penalized. But it wouldn't have more coefficients than data.

[gavinsimpson]

Regardless, this shouldn't crash R, so it might be worth contacting Simon Wood with a reproducible example as long as we're confident that this isn't something that we're doing wrong in MRFtools.

[mluerig]

@gavinsimpson thanks for the additional explanations.

In that regard, one important thing I just realized that it is possible to fit models with multiple observations per species, which I didn't think was possible (other statistical frameworks that account for phylogeny like mvMORPH or some phytools functions don't allow for that). That partly answers my follow-up question.

Moreover, this makes me think whether I should have fit the models that crashed together - in this particular case I would then have had 2 observations per species (same trait measured in different parts of the organism).

[eric-pedersen]

The multiple observations per species case is actually better when using mgcv; it's because the assumption is that each observation $i$ has some underlying mean $\mu_i$ that is predicted from the covariates (which here includes species identity), but the observed value also has idiosyncratic variation (following some distribution, here beta). As such, the assumption here is you have multiple observations from each species.

The other statistical tools you mentioned (to the best of my understanding) assume that you just have a single observation/outcome per species, and the distribution of those traits follows a statistical model with covariates plus a phylogenetic covariance matrix, which includes variation "at the tips" as it were. The standard phylogenetic covariance matrices used assumes that two individuals of the same species should be perfectly correlated (i.e. $Cor(y_1, y_2) = 1$ if $y_1$ and $y_2$ are the same species.

This is easiest to understand when thinking about Normally distributed outcomes. If we have outcomes $y_ij$ from species $j$ a covariate $x_i$ observed for each outcome and some phylogenetic correlation matrix $\Sigma(s)$, times some variance term $\sigma_{p}$

Then the standard phylogenetic regression model assumes that the distribution of $y$ is:

$$y \sim Normal(mean= \mu, Covariance = \sigma_p^2 \Sigma)$$ $$\mu_i = a+bx_i$$

Whereas in the GAM framework, the assumed model is:

$$y \sim Normal(mean= \mu, Covariance = \sigma_{obs}^2I)$$ $$\mu_i = a+bx_i + \beta_j$$ $$\beta \sim Normal(0, \sigma_p\Sigma)$$

Notice the additional variance parameter $\sigma_{obs}$, multiplied by the identity matrix $I$, and that in the second model, the phylogenetic covariance prior is on the $\beta$ coefficients instead of the observations themselves.

This means that the GAM model will generally struggle to converge if you only give it one observation per species and a full-rank covariance matrix (as it is difficult to tell observational noise apart from very fine-scale / high frequency variation in species similarities) but it can fit models with multiple observations per species, as it can model within-species variation.

[eric-pedersen]

Also, for the use-case you talked about, with the same trait measured in different locations: you could model this by including an interaction of LD1 with with body location (something like t2(LD1, location, bs = c("tp","re"), k = c(5,NA)), where the body location is measured as a random effect). That would give you one curve for each location.

Alternatively, if you want to pool information across locations to estimate an overall LD1 effect + location-specific deviations + a phylogenetic effect, this model should do it:

prob ~ s(LD1, k = 5) +                      # overall mean LD1 effect; the "global smoother" in terms of Pedersen et al. 2019 PeerJ
       s(LD1, location, bs = "sz") +        # models location-specific deviations of the LD1 curve from the overall mean. see  `?mgcv::smooth.construct.sz.smooth.spec` for more info on this smoother
       s(location, bs = "re") +             # random effect intercept for location-specific variation
       s(species, bs = "mrf", xt = list(penalty = phylo_penalty), k = 200) # phylogenetic smoother for species-specific intercepts

[mluerig]

thanks very much @eric-pedersen for this extensive description, this really helped my understanding of what is going on under the hood. overall an MRF term seems like a very powerful approach for cases where you want intraspecific variation to be part of your model - this gives me a lot of new ideas!

also, your suggestion for the model seems the way to go. only thing I did was to remove the random term for location, and make it a parametric effect

(I'll leave the comment open in case you want to follow up with Simon Wood about the crash)