Commit 23aac528 authored by Jerome Pasquier's avatar Jerome Pasquier
Browse files

updates longitudinal analyses

parent 5c9d9c8a
......@@ -14,7 +14,7 @@ library(writexl)
options(mc.cores = detectCores() - 1)
# Define subset
# k = 1 -> All observation
# k = 1 -> All observations
# k = 2 -> Inculude only observations until 16.5y
k <- 2
......@@ -33,6 +33,7 @@ dta <- read_xlsx("data-raw/Evolution-december_2021_AI_updated.xlsx",
names(dta) <- trimws(gsub("\\(.+\\)|\\r\\n", "", names(dta)))
# Subset
dta <- dta[dta$Age >= 14, ]
if (k == 2) {
dta <- dta[dta$Age <= 16.5, ]
}
......@@ -49,49 +50,36 @@ rm(dx)
Y <- c("TV", "T", "LH", "FSH", "IGF1", "AMH", "INB")
# Analyses
options(warn = 1)
R <- mclapply(setNames(Y, Y), function(y) {
sdta <- dta %>%
select(ID, DX, Age, one_of(y)) %>%
drop_na()
fml <- as.formula(paste(y, "~ DX * Age + (1 | ID)"))
drop_na() %>%
mutate(Age2 = Age^2, Age3 = Age^3)
#fml <- as.formula(paste(y, "~ DX * poly(Age, 3) + (1 | ID)"))
fml <- as.formula(paste(y, "~ DX * (Age + I(Age^2) + I(Age^3)) + (1 | ID)"))
ctrl <- lmerControl(optimizer ="Nelder_Mead")
DX_list <- c("CDGP", "Partial CHH", "Complete CHH")
DX_list <- setNames(DX_list, sub(" ", "_", DX_list))
ref_ages <- c(0, seq(14, c(20, 16.5)[k], .5)) %>% setNames(., .)
fits <- lapply(ref_ages, function(a) {
lapply(DX_list, function(dx) {
sdta <- sdta %>%
mutate(
DX = relevel(DX, dx),
Age = Age - a
)
fit <- do.call("lmer", list(formula = fml, data = quote(sdta),
control = quote(ctrl)))
tbl <- tidy(fit, conf.int = TRUE) %>%
select(term, estimate, std.error, conf.low, conf.high, p.value) %>%
mutate(
term = sub("^DX", "", term),
term = sub(":", " x ", term),
term = sub("sd__", "SD ", term),
)
if (a != 0) {
tbl <- tbl %>%
mutate(term = sub("Age", paste0("AgeC", a), term))
}
names(tbl)[1] <- paste0("term (ref: ", dx, ")")
list(fit = fit, tbl = tbl)
})
})
tbls <- lapply(fits, function(z) lapply(z, function(w) w$tbl))
fits <- lapply(fits, function(z) lapply(z, function(w) w$fit))
fit <- do.call("lmer", list(formula = fml, data = quote(sdta),
control = quote(ctrl)))
tbl <- tidy(fit, conf.int = TRUE) %>%
select(term, estimate, std.error, conf.low, conf.high, p.value) %>%
mutate(
term = sub("^DX", "", term),
term = sub(":", " x ", term),
term = sub("sd__", "SD ", term),
)
figs <- list()
dgp <- plot_model(fits[[1]][[1]], type = "diag")
dgp <- plot_model(fit, type = "diag")
dgp[[2]] <- dgp[[2]]$ID
figs$diag <- ggarrange(plotlist = dgp, nrow = 2, ncol = 2) %>%
annotate_figure(top = text_grob("Diagnostic plots",
face = "bold", size = 16)) %>%
suppressMessages()
figs$pred1 <- augment(fits[[1]][[1]]) %>%
figs$pred1 <- augment(fit) %>%
#mutate(
# p = `poly(Age, 3)`,
# Age = with(attr(p, "coefs"), alpha[1] + sqrt(norm2)[3] * p[, 1])
#) %>%
group_by(ID) %>%
filter(n() > 1) %>%
ggplot(aes(Age, !!sym(y))) +
......@@ -100,41 +88,30 @@ R <- mclapply(setNames(Y, Y), function(y) {
facet_wrap(~ID) +
labs(title = "Individual predictions") +
theme(legend.position = "bottom", legend.title = element_blank())
ndta <- do.call(rbind, lapply(names(fits[[1]]), function(dx) {
ggpredict(fits[[1]][[dx]], "Age") %>%
as_tibble() %>%
select(x, predicted, conf.low, conf.high) %>%
mutate(dx = dx)
})) %>%
mutate(dx = factor(dx, names(fits[[1]])))
ndta <- ggpredict(fit, c("Age [all]", "DX")) %>%
as_tibble() %>%
select(x, predicted, conf.low, conf.high, group)
figs$pred2 <- ggplot(ndta, aes(x = x, y = predicted)) +
geom_line(aes(colour = dx)) +
geom_ribbon(aes(ymin = conf.low, ymax = conf.high, fill = dx), alpha = 0.3,
show.legend = FALSE) +
geom_line(aes(colour = group)) +
geom_ribbon(aes(ymin = conf.low, ymax = conf.high, fill = group),
alpha = 0.3, show.legend = FALSE) +
labs(x = "Age", y = y, title = "Fixed effects predictions") +
theme(legend.position = "bottom", legend.title = element_blank())
list(fits = fits, tbls = tbls, figs = figs)
list(fit = fit, tbl = tbl, figs = figs)
})
# Any singular fit ?
b <- lapply(R, function(r) lapply(r$fits, function(z) lapply(z, isSingular)))
b <- sapply(R, function(r) isSingular(r$fit))
if (any(unlist(b))) warning("Singular fit")
rm(b)
# Export results - Tables
Z <- lapply(R, function(r) lapply(r$tbls, function(z) {
tbl <- lapply(z, cbind, NA) %>%
append(list(.name_repair = "minimal")) %>%
do.call(bind_cols, .)
names(tbl)[names(tbl) == "NA"] <- ""
return(tbl)
}))
for (y in names(Z)) {
for (y in names(R)) {
f <- paste(y, "regression_coefficients.xlsx", sep = "_")
f <- file.path(outdir, f)
write_xlsx(Z[[y]], f)
write_xlsx(R[[y]]$tbl, f)
}
rm(Z, f, y)
rm(f, y)
# Export results - Figures
for (y in names(R)) {
......
library(broom.mixed)
library(dplyr)
library(ggplot2)
library(lme4)
library(parallel)
library(readxl)
library(tidyr)
library(writexl)
options(mc.cores = detectCores() - 1)
# Set working directory
setwd("~/Projects/Consultations/Pitteloud Nelly (Puberté)")
# Output directory
outdir <-
paste0("results/longitudinal_analyses_", format(Sys.Date(), "%Y%m%d"))
outdir <- file.path(outdir, "nonlinear")
if (!dir.exists(outdir)) dir.create(outdir, recursive = TRUE)
# Load data
dta <- read_xlsx("data-raw/Evolution-december_2021_AI_updated.xlsx",
na = c("", "NA", "-"))
names(dta) <- trimws(gsub("\\(.+\\)|\\r\\n", "", names(dta)))
# Add DX
dta <- dta %>%
select(-DX) %>%
left_join(unique(na.omit(dta[c("ID", "DX")])), by = "ID") %>%
mutate(
DX = sub(", Kallmann", "", DX),
DX = factor(DX, c("CDGP", "Partial CHH", "Complete CHH")),
PCHH = as.numeric(DX == "Partial CHH"),
CCHH = as.numeric(DX == "Complete CHH")
)
# Outcomes
Y <- c("TV", "T", "LH", "FSH", "IGF1")
# Nonlinear model
Model <- function(Age, Asym, Asym2, Asym3, xmid, scal, PCHH, CCHH) {
(Asym + Asym2 * PCHH + Asym3 * CCHH) /
(1 + exp((xmid - Age) / scal))
}
ModelGradient <- deriv(
body(Model)[[2]],
namevec = c("Asym", "Asym2", "Asym3", "xmid", "scal"),
function.arg = Model
)
start_values <- list(
TV = c(Asym = 20, Asym2 = -10, Asym3 = -18, xmid = 16, scal = 2),
T = c(Asym = 20, Asym2 = -15, Asym3 = -18, xmid = 16, scal = 2),
LH = c(Asym = 5, Asym2 = -2.5, Asym3 = -4, xmid = 16, scal = 2),
FSH = c(Asym = 5, Asym2 = -2.5, Asym3 = -4, xmid = 16, scal = 2),
IGF1 = c(Asym = 600, Asym2 = -300, Asym3 = 0, xmid = 16, scal = 2)#,
#AMH = c(Asym = 500, Asym2 = 0, Asym3 = 0, xmid = 16, scal = 2),
#INB = c(Asym = 400, Asym2 = -250, Asym3 = -350, xmid = 16, scal = 2)
)
# Analyses
R <- mclapply(setNames(Y, Y), function(y) {
sdta <- dta %>%
select(ID, DX, PCHH, CCHH, Age, one_of(y)) %>%
drop_na()
fml <- ModelGradient(Age, Asym, Asym2, Asym3, xmid, scal, PCHH, CCHH) ~
Asym | ID
fml <- update(fml, as.formula(paste(y, "~ . ~ .")))
fit <- do.call("nlmer", list(formula = fml, data = quote(sdta),
start = start_values[[y]]))
tbl <- tidy(fit, conf.int = TRUE) %>%
select(term, estimate, std.error, conf.low, conf.high) %>%
mutate(
term = sub(":", " x ", term),
term = sub("sd__", "SD ", term)
)
figs <- list()
figs$pred1 <- augment(fit) %>%
mutate(
DX = c("CDGP", "Partial CHH", "Complete CHH")[1 + PCHH + 2 * CCHH]
) %>%
group_by(ID) %>%
filter(n() > 1) %>%
ggplot(aes(Age, !!sym(y))) +
geom_point() +
geom_line(aes(y = .fitted, colour = DX)) +
facet_wrap(~ID) +
labs(title = "Individual predictions") +
theme(legend.position = "bottom", legend.title = element_blank())
ndta <- expand_grid(
Age = seq(14, 20, 0.1),
DX = c("CDGP", "Partial CHH", "Complete CHH")
) %>% mutate(
PCHH = as.numeric(DX == "Partial CHH"),
CCHH = as.numeric(DX == "Complete CHH")
) %>% cbind(as.data.frame(t(setNames(tbl[["estimate"]], tbl[["term"]])))) %>%
rowwise() %>%
mutate(predicted = Model(Age, Asym, Asym2, Asym3, xmid, scal, PCHH, CCHH))
figs$pred2 <- ggplot(ndta, aes(x = Age, y = predicted)) +
geom_line(aes(colour = DX)) +
labs(x = "Age", y = y, title = "Fixed effects predictions") +
theme(legend.position = "bottom", legend.title = element_blank())
tbl <- tbl %>%
mutate(
term = sub("^Asym2", "AsymPartCHH", term),
term = sub("^Asym3", "AsymCompCHH", term)
)
list(fit = fit, tbl = tbl, figs = figs)
})
# Any singular fit ?
b <- sapply(R, function(r) isSingular(r$fit))
if (any(unlist(b))) warning("Singular fit")
rm(b)
# Export results - Tables
for (y in names(R)) {
f <- paste(y, "regression_coefficients.xlsx", sep = "_")
f <- file.path(outdir, f)
write_xlsx(R[[y]]$tbl, f)
}
rm(f, y)
# Export results - Figures
for (y in names(R)) {
tiff(filename = file.path(outdir, paste0(y, "_individual_predictions.tiff")),
height = 3600, width = 5400, res = 512, compression = "zip")
print(R[[y]]$figs$pred1)
dev.off()
tiff(filename = file.path(outdir,
paste0(y, "_fixed_effects_predictions.tiff")),
height = 3600, width = 5400, res = 1024, compression = "zip")
print(R[[y]]$figs$pred2)
dev.off()
}
rm(y)
# Session Info
sink(file.path(outdir, "sessionInfo.txt"))
print(sessionInfo(), locale = FALSE)
sink()
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment