Skip to content

Instantly share code, notes, and snippets.

@WillForan
Last active July 21, 2022 20:10
Show Gist options
  • Save WillForan/4469423e225c75e0ed1c8a3c5813f5ac to your computer and use it in GitHub Desktop.
Save WillForan/4469423e225c75e0ed1c8a3c5813f5ac to your computer and use it in GitHub Desktop.
library(ggplot)
library(gamm4)
library(dplyr)
library(checkmate) # for expect_subset dataframe names test
met_res_for_roi <- function(MRSI_input, this_roi, met_name, minage=18, maxage=25) {
# have columns we need
expect_subset(c("roi", "age", "dateNumeric", met_name), names(MRSI_input))
roi_met <- MRSI_input %>%
filter(roi %in% this_roi, age >=minage, age <=maxage) %>%
mutate(met = .data[[met_name]], # map whatever metabolite we're using (met_name) onto the 'met' column
gamresids = gam(met ~ s(age, k=3), data=.) %>% residuals,
dateAdjust = gam(gamresids ~ s(dateNumeric, k=3), data=.) %>% predict,
met_adj = met - dateAdjust)
return(roi_met)
}
plot_met_adjusted <- function(roi_met) {
# have columns we need
expect_subset(c("age", "gamresids", "dateAdjust", "met", "met_adj", "visitnum", "subjID"),
names(roi_met))
qassert(roi_met$met_adj, "n+") # nonzero numeric vector in met_adj
cowplot::plot_grid(
ggplot(roi_met) + aes(x=age, y=gamresids) + geom_smooth(method="loess") + geom_point(),
ggplot(roi_met) + aes(x=vdate, y=dateAdjust) + geom_point() + geom_smooth(method="loess"),
ggplot(roi_met) + aes(x=vdate, y=met_adj, color=visitnum, group=visitnum) + geom_point() + geom_smooth(method="loess"),
ggplot(roi_met) + aes(x=age, y=met_adj) + geom_point() + stat_smooth(method='loess', se=F) + geom_path(aes(group = subjID)),
ggplot(roi_met) + aes(x=age, y=met) + geom_point() + stat_smooth(method='loess', se=F) + geom_path(aes(group = subjID))
)
}
test_data <- data.frame(subjID=1:32, roi=1, age=rep(10:25,2), gaba=randu[1:32,1], dateNumeric=sample(rep(1:8,4)),visitnum=1) %>% mutate(vdate=dateNumeric)
adj_df <- met_res_for_roi(test_data, 1, "gaba")
plot_met_adjusted(adj_df)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment