Author: Andrew A. Chen, chenandr@musc.edu
Maintainer: Andrew A. Chen, chenandr@musc.edu
License: MIT License
gamlss2charts includes additional functions based on gamlss2, with backwards compatibility for gamlss (WIP). Currently, our package contains predict_score for out-of-sample scoring (i.e. centile scoring) based on previous works (see references).
This package is still a work-in-progress and will be continually supported and updated to include additional functionalities.
The latest version can be installed via remotes by running the following code
# install.packages("remotes")
remotes::install_github("andy1764/gamlss2charts")
Then, you can load this package via
library(gamlss2charts)
A toy example for getting out-of-sample centile scores in the iris dataset is shown below
library(gamlss2)
ft <- gamlss2(Sepal.Length ~ Sepal.Width + Petal.Length + Petal.Width + Species | Species,
family = BCCG(), data = iris[1:100,])
predict_score(ft, iris[-(1:100),], adjust = FALSE, rm.term = "Species") |> hist()
predict_score(ft, iris[-(1:100),], rm.term = "Species") |> hist()
# with two test batches, fixed effect in adjustment fit
ft <- gamlss2(Sepal.Length ~ Sepal.Width + Petal.Length + Petal.Width,
family = BCCG(), data = iris[1:50,])
predict_score(ft, iris[-(1:50),], adjust = FALSE) |> hist()
predict_score(ft, iris[-(1:50),],
newformula = y ~ offset(mu) + Species | offset(sigma) + Species) |> hist()
Our package supports basic functionality for the gamlss package. Workarounds are still needed for certain model terms (including random effects, see below)
library(gamlss)
ft <- gamlss(Sepal.Length ~ Sepal.Width + Petal.Length + Petal.Width + Species, ~ Species,
family = BCCG(), data = iris[1:100,])
predict_score(ft, iris[-(1:100),], adjust = FALSE, rm.term = "Species") |> hist()
predict_score(ft, iris[-(1:100),], rm.term = "Species") |> hist()
# using random effects
testris <- iris
testris[-(1:100),]$Species <- "setosa" # set test batch to a batch in the train set
ft <- gamlss(Sepal.Length ~ Sepal.Width + Petal.Length + Petal.Width + random(Species), ~ random(Species),
family = BCCG(), data = testris[1:100,])
predict_score(ft, testris[-(1:100),], adjust = FALSE) |> hist()
predict_score(ft, testris[-(1:100),]) |> hist()
Our package supports the model format from the Bethlehem et al. (2022) charts. Below is a working example after cloning https://github.com/brainchart/Lifespan to ~/GitHub/Lifespan/:
library(gamlss2)
source("~/GitHub/Lifespan/102.gamlss-recode.r")
source("~/GitHub/Lifespan/301.functions.r")
FIT <- readRDS("~/GitHub/Lifespan/Share/OriginalModels/FIT_GMV.rds")
# generate age and sex
set.seed(8888)
newdat <- data.frame(age = rnorm(200, 50, 5),
sex = c(rep("Male", 100), rep("Female", 100)))
newdat$AgeTransformed <- log(newdat$age*365.25) # 2022 charts use log(age in days)
POP.CURVE.RAW <- newdat[,c("AgeTransformed", "sex")]
CURVE <- Apply.Param(NEWData=POP.CURVE.RAW, FITParam=FIT$param, Add.Derivative=FALSE)
params <- CURVE[c("mu.pop", "sigma.pop", "nu.pop")]
names(params) <- c("mu", "sigma", "nu")
# generate GMV (in units of 10,000 mm3) values from chart parameters
newdat$GMV <- rGGalt(200, mu = exp(params$mu), sigma = exp(params$sigma), nu = params$nu)
predict_score(list(new = params), newdata = newdat, feat = "GMV", family = GGalt) |> hist()
Please cite the following papers for out-of-sample centile scoring:
Dinga, R., Fraza, C. J., Bayer, J. M. M., Kia, S. M., Beckmann, C. F., & Marquand, A. F. (2021). Normative modeling of neuroimaging data using generalized additive models of location scale and shape (p. 2021.06.14.448106). bioRxiv. https://doi.org/10.1101/2021.06.14.448106
Bethlehem, R. A. I., Seidlitz, J., White, S. R., Vogel, J. W., Anderson, K. M., Adamson, C., Adler, S., Alexopoulos, G. S., Anagnostou, E., Areces-Gonzalez, A., Astle, D. E., Auyeung, B., Ayub, M., Bae, J., Ball, G., Baron-Cohen, S., Beare, R., Bedford, S. A., Benegal, V., … Alexander-Bloch, A. F. (2022). Brain charts for the human lifespan. Nature, 604(7906), Article 7906. https://doi.org/10.1038/s41586-022-04554-y