OSCR

Genetic drivers of progression in Alzheimer's disease are distinct from disease risk.

Code ↔ Paper

3 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 3 matches
  1. [1] § Methods › Genomic quality control and imputation ↔ Progression_Genetic_Modelling/Genetic QC/genetic_qc_consolidated.R, lines 1–59 · score 0.85 · Hardy Weinberg Equilibrium, minor allele frequency, missing genotyping, deviations, QC, quality
  2. [2] § Methods › Patient filtering ↔ Progression_Genetic_Modelling/Clinical Data Modelling and Filtering/comprehensive_analysis.qmd, lines 84–105 · score 0.70 · Pittsburgh compound, positive patients, CSF, AV45, PiB, filtered
  3. [3] § Methods › Mixed-effects modelling ↔ Progression_Genetic_Modelling/Clinical Data Modelling and Filtering/comprehensive_analysis.qmd, lines 150–186 · score 0.56 · random intercept, baseline MMSE, PC1, PC2, education, age

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The authors' code

Quarto · 1,071 lines · 36 KB · no license · 2 matches

  1. ---
  2. title: "Exploratory Data Analysis, Model Optimisation and Clinical Filtering"
  3. author:
  4. - "CCOHEN"
  5. - "MSHOAI"
  6. date: today
  7. format:
  8. html:
  9. theme: cosmo
  10. toc: true
  11. toc-depth: 3
  12. toc-location: left
  13. number-sections: true
  14. code-fold: true
  15. code-tools: true
  16. df-print: kable
  17. fig-width: 10
  18. fig-height: 6
  19. embed-resources: true
  20. execute:
  21. warning: false
  22. message: false
  23. cache: false
  24. ---
  26. # Introduction
  28. This document presents a comprehensive analysis of longitudinal cognitive decline in Alzheimer's Disease (AD) patients from the Alzheimer's Disease Neuroimaging Initiative (ADNI). The analysis includes:
  30. - Quality control and patient filtering based on amyloid-beta biomarkers
  31. - Exploratory data analysis and descriptive statistics
  32. - Mixed-effects model optimization and comparison
  33. - Integration with polygenic risk scores
  35. # Data Loading and Preparation
  37. ```{r setup, warning=FALSE, message=FALSE}
  38. # Load required libraries
  39. library(tidyverse)
  40. library(tseries)
  41. library(visdat)
  42. library(ggplot2)
  43. library(dplyr)
  44. library(gridExtra)
  45. library(Amelia)
  46. library(corrplot)
  47. library(knitr)
  48. library(lme4)
  49. library(lmerTest)
  50. library(performance)
  51. library(fitdistrplus)
  52. library(glmmTMB)
  54. # Set options for output formatting
  55. options(knitr.kable.NA = '', digits = 2)
  56. ```
  58. ```{r load_data}
  59. # Read ADNIMERGE clinical data
  60. adnimerge <- readRDS("./QCjune/clinicaladni/adnimerge.rds")
  62. # Select relevant columns
  63. dataset <- adnimerge[, c(4, 9, 10, 11, 114, 27, 69, 61, 8, 18, 110, 17, 109, 20, 105, 14, 15, 1)]
  64. colnames(dataset) <- c("Patient_ID", "Age", "Gender", "Education", "Months_since_bl",
  65. "MMSE", "MMSE.bl", "DX", "DX.bl", "AV45", "AV45.bl", "PIB", "PIB.bl",
  66. "CSF_AB", "CSF_AB.bl", "Marital", "APOE4", "RID")
  68. # Clean CSF Amyloid-beta values
  69. dataset$CSF_AB.bl <- gsub(">", "", dataset$CSF_AB.bl)
  70. dataset$CSF_AB.bl <- gsub("<", "", dataset$CSF_AB.bl)
  71. dataset$CSF_AB.bl <- as.numeric(dataset$CSF_AB.bl)
  72. dataset$CSF_AB <- gsub(">", "", dataset$CSF_AB)
  73. dataset$CSF_AB <- gsub("<", "", dataset$CSF_AB)
  74. dataset$CSF_AB <- as.numeric(dataset$CSF_AB)
  75. ```
  77. ## Initial Dataset Summary
  79. ```{r initial_summary}
  80. cat("Total unique patients:", length(unique(dataset$Patient_ID)), "\n")
  81. cat("Total datapoints:", nrow(dataset), "\n")
  82. ```
  84. # Patient Filtering Pipeline
  86. ## Step 1: Amyloid-Beta Positive Selection {#step1}
  88. Patients are classified as amyloid-beta positive if they meet at least one of the following criteria at baseline:
  90. - **AV45 (Florbetapir PET)**: ≥ 1.1
  91. - **PIB (Pittsburgh Compound B PET)**: ≥ 1.5
  92. - **CSF Amyloid-beta**: ≤ 550 pg/mL
  94. ```{r filter_amyloid}
  95. # Define thresholds
  96. th_amyl_AV45 <- 1.1
  97. th_amyl_PIB <- 1.5
  98. th_amyl_CSF <- 550
  100. # Filter for amyloid-beta positive patients
  101. df_amyloid <- subset(dataset, AV45.bl >= th_amyl_AV45 | PIB.bl >= th_amyl_PIB | CSF_AB.bl <= th_amyl_CSF)
  103. cat("Amyloid-beta positive patients:", length(unique(df_amyloid$Patient_ID)), "\n")
  104. cat("Number of datapoints:", nrow(df_amyloid), "\n")
  105. ```
  107. ### Amyloid-Beta Distribution
  109. ```{r amyloid_dist, fig.height=4}
  110. # Create binary amyloid status variable
  111. dataset$Amyloid_Status <- ifelse(
  112. dataset$AV45.bl >= th_amyl_AV45 | dataset$PIB.bl >= th_amyl_PIB | dataset$CSF_AB.bl <= th_amyl_CSF,
  113. "Positive", "Negative"
  114. )
  116. # Plot distribution
  117. ggplot(dataset %>% group_by(Patient_ID) %>% slice(1), aes(x = Amyloid_Status, fill = Amyloid_Status)) +
  118. geom_bar() +
  119. geom_text(stat='count', aes(label=..count..), vjust=-0.5) +
  120. labs(title = "Distribution of Amyloid-Beta Status",
  121. x = "Amyloid-Beta Status", y = "Number of Patients") +
  122. scale_fill_manual(values = c("Positive" = "#E74C3C", "Negative" = "#3498DB")) +
  123. theme_minimal() +
  124. theme(legend.position = "none", text = element_text(size = 12))
  125. ```
  127. # Model Optimization on Amyloid-Positive Cohort
  129. Before applying additional filtering, we test multiple mixed-effects models on the amyloid-positive cohort to identify the optimal model structure.
  131. ```{r prepare_model_data}
  132. # Prepare data for modeling (amyloid+ only, at this stage)
  133. df_model1 <- df_amyloid
  134. df_model1$Gender <- as.factor(df_model1$Gender)
  135. df_model1$Months_since_bl <- as.numeric(as.character(df_model1$Months_since_bl))
  137. # For this initial modeling, we need PC1 and PC2
  138. # Read PCA data
  139. pcs <- read.table("./QCjune/adnimerge_genomes/FINALADNI_nohapmap_pca.eigenvec")[, 2:4]
  140. colnames(pcs) <- c("PTID", "PC1", "PC2")
  142. # Merge with clinical data
  143. df_model1$PTID <- as.character(df_model1$Patient_ID)
  144. df_model1_pcs <- inner_join(df_model1, pcs, by = "PTID", multiple = "all")
  146. cat("Amyloid+ patients with genetic data for modeling:", length(unique(df_model1_pcs$Patient_ID)), "\n")
  147. cat("Datapoints for modeling:", nrow(df_model1_pcs), "\n")
  148. ```
  150. ## Model 1-9 Comparison
  152. We test nine hierarchical models of increasing complexity:
  154. ```{r model_comparison, results='hide'}
  155. # Model 1: Time only
  156. model1 <- lm(MMSE ~ Months_since_bl, data = df_model1_pcs)
  158. # Model 2: Time + Baseline age
  159. model2 <- lm(MMSE ~ Months_since_bl + Age, data = df_model1_pcs)
  161. # Model 3: Time + Current age
  162. model3 <- lm(MMSE ~ Months_since_bl + Age, data = df_model1_pcs)
  164. # Model 4: Time + Baseline MMSE + Baseline age
  165. model4 <- lm(MMSE ~ MMSE.bl + Months_since_bl + Age, data = df_model1_pcs)
  167. # Model 5: Random intercept
  168. model5 <- lmer(MMSE ~ MMSE.bl + Age + Months_since_bl + (1|Patient_ID),
  169. data = df_model1_pcs, REML = TRUE)
  171. # Model 6: Random intercept + Education
  172. model6 <- lmer(MMSE ~ MMSE.bl + Age + Months_since_bl + (1|Patient_ID) + Education,
  173. data = df_model1_pcs, REML = TRUE)
  175. # Model 7: Random intercept + Education + Gender
  176. model7 <- lmer(MMSE ~ MMSE.bl + Age + Months_since_bl + (1|Patient_ID) + Education + Gender,
  177. data = df_model1_pcs, REML = TRUE)
  179. # Model 8: Random intercept + Education + Gender + PCs
  180. model8 <- lmer(MMSE ~ MMSE.bl + Age + Months_since_bl + (1|Patient_ID) + Education + Gender + PC1 + PC2,
  181. data = df_model1_pcs, REML = TRUE)
  183. # Model 9: Random slope and intercept + Education + Gender + PCs
  184. model9 <- lmer(MMSE ~ MMSE.bl + Age + Months_since_bl + (Months_since_bl|Patient_ID) +
  185. Education + Gender + PC1 + PC2, data = df_model1_pcs, REML = TRUE)
  186. ```
  188. ```{r model_performance}
  189. # Compile model performance metrics
  190. models <- list(model1, model2, model3, model4, model5, model6, model7, model8, model9)
  191. model_names <- paste0("Model ", 1:9)
  193. # Get performance metrics manually
  194. aics <- sapply(models, AIC)
  195. bics <- sapply(models, BIC)
  196. r2_vals <- sapply(models, function(m) {
  197. r2_result <- tryCatch(r2(m), error = function(e) list(R2 = NA))
  198. if ("R2" %in% names(r2_result)) return(r2_result$R2)
  199. if ("R2_conditional" %in% names(r2_result)) return(r2_result$R2_conditional)
  200. return(NA)
  201. })
  203. performance_table <- data.frame(
  204. Model = model_names,
  205. AIC = aics,
  206. BIC = bics,
  207. R2 = r2_vals
  208. )
  210. kable(performance_table,
  211. caption = "Performance Comparison of Models 1-9 on Amyloid-Positive Cohort",
  212. digits = 2)
  213. ```
  215. ::: {.callout-note}
  216. ## Best Model Selection
  218. Based on AIC, BIC, and R² metrics, **Model 9** (random slopes and intercepts with full covariates) provides the best fit. This model will be used for subsequent analyses after full filtering.
  219. :::
  221. # Continued Patient Filtering
  223. ## Step 2: Remove Missing MMSE and Require ≥3 Datapoints
  225. ```{r filter_mmse}
  226. df <- df_amyloid
  228. cat("Patients before filtering:", length(unique(df$Patient_ID)), "\n")
  229. cat("Datapoints before filtering:", nrow(df), "\n\n")
  231. # Remove rows with missing MMSE
  232. df <- df[-which(is.na(df$MMSE)), ]
  234. # Require at least 3 datapoints per patient
  235. Keep <- df %>%
  236. group_by(Patient_ID) %>%
  237. summarize(datapoints = n()) %>%
  238. filter(datapoints >= 3)
  239. df <- df[df$Patient_ID %in% Keep$Patient_ID, ]
  241. cat("Patients after filtering:", length(unique(df$Patient_ID)), "\n")
  242. cat("Datapoints after filtering:", nrow(df), "\n")
  243. ```
  245. ## Step 3: Remove "Always Cognitively Normal" Patients
  247. ```{r filter_always_cn}
  248. # Identify patients who are CN at baseline
  249. PTID_CNbl <- unique(df[df$Months_since_bl == 0 & df$DX == "CN", ]$Patient_ID)
  250. PTID_CNbl_df <- df[df$Patient_ID %in% PTID_CNbl, ]
  252. # Find those who progress to MCI or Dementia
  253. PTID_CNthenMCI <- unique(PTID_CNbl_df[which(PTID_CNbl_df$DX == "MCI" | PTID_CNbl_df$DX == "Dementia"), ]$Patient_ID)
  255. # Remove always-CN patients
  256. PTID_alwaysCN <- PTID_CNbl[-which(PTID_CNbl %in% PTID_CNthenMCI)]
  257. df_clean <- df[-which(df$Patient_ID %in% PTID_alwaysCN), ]
  259. cat("Patients not always cognitively normal:", length(unique(df_clean$Patient_ID)), "\n")
  260. cat("Datapoints:", nrow(df_clean), "\n")
  261. ```
  263. ## Step 4: Remove Patients with Minimum MMSE ≥ 28
  265. ```{r filter_min_mmse}
  266. # Calculate minimum MMSE per patient
  267. minMMSEperpatient <- df_clean %>%
  268. group_by(Patient_ID) %>%
  269. summarise(minMMSE = min(MMSE))
  271. min28 <- minMMSEperpatient$Patient_ID[which(minMMSEperpatient$minMMSE >= 28)]
  272. df_clean <- df_clean[!df_clean$Patient_ID %in% min28, ]
  274. cat("Patients with minimum MMSE < 28:", length(unique(df_clean$Patient_ID)), "\n")
  275. cat("Datapoints:", nrow(df_clean), "\n")
  276. ```
  278. ## Step 5: Remove Patients with Last MMSE 28-30 OR Last DX = CN
  280. ```{r filter_last_status}
  281. # Calculate last MMSE and DX per patient
  282. lastMMSE <- c()
  283. lastDX <- c()
  284. for (i in unique(df_clean$Patient_ID)) {
  285. d <- subset(df_clean, Patient_ID == i)
  286. d_sort <- d[order(d$Months_since_bl, decreasing = FALSE), ]
  287. lastMMSE <- c(lastMMSE, d_sort$MMSE[length(d_sort$MMSE)])
  288. lastDX <- c(lastDX, d_sort$DX[length(d_sort$DX)])
  289. }
  291. # Remove patients with last MMSE 28-30
  292. df_lastMMSE30 <- subset(df_clean, Patient_ID %in% unique(df_clean$Patient_ID)[which(lastMMSE == 30 | lastMMSE == 29 | lastMMSE == 28)])
  293. df_clean <- df_clean[!df_clean$Patient_ID %in% df_lastMMSE30$Patient_ID, ]
  295. # Remove patients with last DX = CN
  296. df_lastDXCN <- subset(df_clean, Patient_ID %in% unique(df_clean$Patient_ID)[which(lastDX == "CN")])
  297. df_clean <- df_clean[!df_clean$Patient_ID %in% df_lastDXCN$Patient_ID, ]
  299. cat("Patients without ceiling effects:", length(unique(df_clean$Patient_ID)), "\n")
  300. cat("Datapoints:", nrow(df_clean), "\n")
  301. ```
  303. ## Step 6: Remove Datapoints Before MMSE = 30
  305. For patients who achieved MMSE of 30 at any point, we remove all observations before the last occurrence of MMSE = 30. This focuses the analysis on the decline phase.
  307. ```{r filter_before_30}
  308. before_count <- nrow(df_clean)
  309. df_clean$dpid <- seq(1, nrow(df_clean))
  310. rm_rows <- c()
  312. for (i in unique(df_clean$Patient_ID)) {
  313. d <- subset(df_clean, Patient_ID == i)
  314. d_sort <- d[order(d$Months_since_bl, decreasing = FALSE), ]
  315. instances <- c()
  317. for (n in 1:(length(d_sort$MMSE) - 1)) {
  318. if (d_sort$MMSE[n] == 30) {
  319. instances <- c(instances, n)
  320. }
  321. }
  323. if (length(instances) > 0) {
  324. rm <- d_sort$dpid[1:(max(instances) - 1)]
  325. rm_rows <- c(rm_rows, rm)
  326. }
  327. }
  329. df_clean <- df_clean[!df_clean$dpid %in% rm_rows, ]
  330. df_clean$dpid <- NULL
  332. cat("Datapoints removed:", before_count - nrow(df_clean), "\n")
  333. cat("Patients remaining:", length(unique(df_clean$Patient_ID)), "\n")
  334. cat("Datapoints remaining:", nrow(df_clean), "\n")
  335. ```
  337. ## Step 7: Keep Only Last 5 Observations Per Patient
  339. This critical step focuses the analysis on recent disease trajectory.
  341. ```{r filter_last_5}
  342. before_count <- nrow(df_clean)
  343. df_clean$dpid <- seq(1, nrow(df_clean))
  344. rm_rows <- c()
  346. for (i in unique(df_clean$Patient_ID)) {
  347. d <- subset(df_clean, Patient_ID == i)
  348. d_sort <- d[order(d$Months_since_bl, decreasing = FALSE), ]
  350. if (length(d$dpid) > 5) {
  351. rm <- d_sort$dpid[1:(length(d_sort$dpid) - 5)]
  352. rm_rows <- c(rm_rows, rm)
  353. }
  354. }
  356. df_clean <- df_clean[!df_clean$dpid %in% rm_rows, ]
  357. df_clean$dpid <- NULL
  359. cat("Datapoints removed:", before_count - nrow(df_clean), "\n")
  360. cat("Datapoints remaining:", nrow(df_clean), "\n")
  361. ```
  363. ## Step 8: Require ≥3 Datapoints (After Pruning)
  365. ```{r filter_3pts_final}
  366. Keep <- df_clean %>%
  367. group_by(Patient_ID) %>%
  368. summarize(datapoints = n()) %>%
  369. filter(datapoints >= 3)
  370. df_clean <- df_clean[df_clean$Patient_ID %in% Keep$Patient_ID, ]
  372. cat("Final clinically filtered patients:", length(unique(df_clean$Patient_ID)), "\n")
  373. cat("Final datapoints:", nrow(df_clean), "\n")
  374. ```
  376. ::: {.callout-important}
  377. ## Filtering Summary
  379. Starting with **`r length(unique(dataset$Patient_ID))`** patients, the filtering pipeline identified **`r length(unique(df_clean$Patient_ID))`** patients with:
  381. - Confirmed amyloid-beta pathology
  382. - Sufficient longitudinal data (3-5 recent observations)
  383. - Evidence of cognitive decline
  384. :::
  386. # Descriptive Statistics
  388. ## Missing Data Visualization
  390. ```{r missing_data, fig.height=8}
  391. # Visualize missing data patterns
  392. missmap(df_clean[, c("MMSE", "Age", "Education", "Gender", "DX", "AV45", "PIB", "CSF_AB", "APOE4")],
  393. main = "Missing Data Map - Filtered Cohort",
  394. col = c("lightcoral", "lightblue"),
  395. legend = TRUE,
  396. x.cex = 0.8,
  397. y.cex = 0.6)
  398. ```
  400. ## Baseline Characteristics
  402. ```{r baseline_char}
  403. # Get baseline characteristics (one row per patient)
  404. df_baseline <- df_clean %>%
  405. group_by(Patient_ID) %>%
  406. arrange(Months_since_bl) %>%
  407. slice(1) %>%
  408. ungroup()
  410. # Summary statistics
  411. baseline_stats <- data.frame(
  412. Variable = c("Age", "Education", "MMSE"),
  413. Mean = c(mean(df_baseline$Age, na.rm = TRUE),
  414. mean(df_baseline$Education, na.rm = TRUE),
  415. mean(df_baseline$MMSE, na.rm = TRUE)),
  416. SD = c(sd(df_baseline$Age, na.rm = TRUE),
  417. sd(df_baseline$Education, na.rm = TRUE),
  418. sd(df_baseline$MMSE, na.rm = TRUE)),
  419. Median = c(median(df_baseline$Age, na.rm = TRUE),
  420. median(df_baseline$Education, na.rm = TRUE),
  421. median(df_baseline$MMSE, na.rm = TRUE)),
  422. Min = c(min(df_baseline$Age, na.rm = TRUE),
  423. min(df_baseline$Education, na.rm = TRUE),
  424. min(df_baseline$MMSE, na.rm = TRUE)),
  425. Max = c(max(df_baseline$Age, na.rm = TRUE),
  426. max(df_baseline$Education, na.rm = TRUE),
  427. max(df_baseline$MMSE, na.rm = TRUE))
  428. )
  430. kable(baseline_stats, caption = "Baseline Characteristics (Continuous Variables)", digits = 2)
  432. # Categorical variables
  433. cat_summary <- data.frame(
  434. Variable = c("Gender (Male)", "Gender (Female)",
  435. "DX.bl (AD)", "DX.bl (MCI)", "DX.bl (CN)", "DX.bl (Dementia)"),
  436. N = c(sum(df_baseline$Gender == "Male", na.rm = TRUE),
  437. sum(df_baseline$Gender == "Female", na.rm = TRUE),
  438. sum(df_baseline$DX.bl == "AD", na.rm = TRUE),
  439. sum(df_baseline$DX.bl == "MCI" | df_baseline$DX.bl == "EMCI" | df_baseline$DX.bl == "LMCI", na.rm = TRUE),
  440. sum(df_baseline$DX.bl == "CN", na.rm = TRUE),
  441. sum(df_baseline$DX.bl == "Dementia", na.rm = TRUE)),
  442. Percentage = c(100 * sum(df_baseline$Gender == "Male", na.rm = TRUE) / nrow(df_baseline),
  443. 100 * sum(df_baseline$Gender == "Female", na.rm = TRUE) / nrow(df_baseline),
  444. 100 * sum(df_baseline$DX.bl == "AD", na.rm = TRUE) / nrow(df_baseline),
  445. 100 * sum(df_baseline$DX.bl == "MCI" | df_baseline$DX.bl == "EMCI" | df_baseline$DX.bl == "LMCI", na.rm = TRUE) / nrow(df_baseline),
  446. 100 * sum(df_baseline$DX.bl == "CN", na.rm = TRUE) / nrow(df_baseline),
  447. 100 * sum(df_baseline$DX.bl == "Dementia", na.rm = TRUE) / nrow(df_baseline))
  448. )
  450. kable(cat_summary, caption = "Baseline Characteristics (Categorical Variables)", digits = 2)
  451. ```
  453. ## Age Distribution by Diagnosis and Gender
  455. ```{r age_distributions, fig.height=8}
  456. # Plot age distribution by baseline diagnosis
  457. p1 <- ggplot(df_baseline, aes(x = Age, fill = DX.bl)) +
  458. geom_histogram(binwidth = 5, alpha = 0.7) +
  459. labs(title = "Age Distribution by Baseline Diagnosis",
  460. x = "Age (years)", y = "Count", fill = "Baseline DX") +
  461. theme_minimal() +
  462. theme(text = element_text(size = 11))
  464. # Plot age distribution by gender
  465. p2 <- ggplot(df_baseline, aes(x = Age, fill = Gender)) +
  466. geom_histogram(binwidth = 5, alpha = 0.7, position = "identity") +
  467. labs(title = "Age Distribution by Gender",
  468. x = "Age (years)", y = "Count", fill = "Gender") +
  469. scale_fill_manual(values = c("Male" = "#3498DB", "Female" = "#E74C3C")) +
  470. theme_minimal() +
  471. theme(text = element_text(size = 11))
  473. # Boxplot of age by gender with t-test
  474. t_test_result <- t.test(Age ~ Gender, data = df_baseline)
  475. p3 <- ggplot(df_baseline, aes(x = Gender, y = Age, fill = Gender)) +
  476. geom_boxplot(alpha = 0.7) +
  477. labs(title = paste0("Age by Gender (p = ", format(t_test_result$p.value, digits = 3, scientific = TRUE), ")"),
  478. x = "Gender", y = "Age (years)") +
  479. scale_fill_manual(values = c("Male" = "#3498DB", "Female" = "#E74C3C")) +
  480. theme_minimal() +
  481. theme(legend.position = "none", text = element_text(size = 11))
  483. # MMSE trajectory over time by gender
  484. p4 <- ggplot(df_clean, aes(x = Months_since_bl, y = MMSE, color = Gender)) +
  485. geom_point(alpha = 0.3, size = 0.8) +
  486. geom_smooth(method = "lm", se = TRUE, size = 1.2) +
  487. labs(title = "MMSE Trajectory Over Time by Gender",
  488. x = "Months Since Baseline", y = "MMSE Score", color = "Gender") +
  489. scale_color_manual(values = c("Male" = "#3498DB", "Female" = "#E74C3C")) +
  490. theme_minimal() +
  491. theme(text = element_text(size = 11))
  493. # Arrange plots
  494. grid.arrange(p1, p2, p3, p4, ncol = 2)
  495. ```
  497. ## Age and Progression Analysis
  499. ### Progression Rate by Age Group
  501. ```{r age_progression, fig.height=6}
  502. # Calculate progression slopes by age groups
  503. agerange <- c(50, 60, 65, 70, 75, 80, 85, 95)
  504. slopes <- c()
  505. age_group <- c()
  507. for (i in 1:(length(agerange) - 1)) {
  508. low <- agerange[i]
  509. high <- agerange[i + 1]
  510. dfa <- df_clean[which(df_clean$Age >= low & df_clean$Age < high), ]
  512. for (id in unique(dfa$Patient_ID)) {
  513. d <- dfa[which(dfa$Patient_ID == id), ]
  514. if (nrow(d) >= 2) {
  515. slopes <- c(slopes, coef(lm(MMSE ~ Months_since_bl, data = d))[2])
  516. age_group <- c(age_group, paste(low, "to", high))
  517. }
  518. }
  519. }
  521. age_prog_df <- data.frame(Age_Group = age_group, Progression_Rate = slopes)
  523. # Box plot
  524. ggplot(age_prog_df, aes(x = Age_Group, y = Progression_Rate)) +
  525. geom_boxplot(fill = "#3498DB", alpha = 0.7) +
  526. geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  527. labs(title = "Cognitive Decline Rate by Age Group",
  528. x = "Age Group (years)", y = "MMSE Change per Month",
  529. caption = paste("ANOVA p-value:", format(anova(lm(Progression_Rate ~ Age_Group, data = age_prog_df))$Pr[1], digits = 3))) +
  530. theme_minimal() +
  531. theme(axis.text.x = element_text(angle = 45, hjust = 1), text = element_text(size = 11))
  532. ```
  534. ### Progression by Age of Onset
  536. ```{r age_onset, fig.height=4}
  537. # Analyze progression by age of onset (for dementia patients)
  538. df_dementia <- df_clean[which(df_clean$DX == "Dementia"), ]
  539. slopes_dem <- c()
  540. onset_age <- c()
  542. for (i in unique(df_dementia$Patient_ID)) {
  543. d <- subset(df_dementia, Patient_ID == i)
  544. if (nrow(d) >= 2) {
  545. slopes_dem <- c(slopes_dem, coef(lm(MMSE ~ Months_since_bl, data = d))[2])
  546. onset_age <- c(onset_age, min(d$Age))
  547. }
  548. }
  550. onset_df <- data.frame(Onset_Age = onset_age, Progression_Rate = slopes_dem)
  551. onset_df$Onset_Category <- cut(onset_df$Onset_Age,
  552. breaks = c(-Inf, 65, 75, Inf),
  553. labels = c("EOAD (<65)", "MOAD (65-75)", "LOAD (>75)"))
  555. # Plot
  556. p1 <- ggplot(onset_df, aes(x = Onset_Age, y = Progression_Rate)) +
  557. geom_point(alpha = 0.6, color = "#E74C3C") +
  558. geom_smooth(method = "lm", se = TRUE, color = "#3498DB") +
  559. labs(title = "Progression Rate vs Age of Onset",
  560. x = "Age of Onset (years)", y = "MMSE Change per Month") +
  561. theme_minimal()
  563. p2 <- ggplot(onset_df, aes(x = Onset_Category, y = Progression_Rate, fill = Onset_Category)) +
  564. geom_boxplot(alpha = 0.7) +
  565. labs(title = "Progression by Onset Category",
  566. x = "Age of Onset Category", y = "MMSE Change per Month") +
  567. scale_fill_brewer(palette = "Set2") +
  568. theme_minimal() +
  569. theme(legend.position = "none")
  571. grid.arrange(p1, p2, ncol = 2)
  573. # Statistical test
  574. lm_onset <- lm(Progression_Rate ~ Onset_Age, data = onset_df)
  575. cat("\nLinear regression p-value:", format(summary(lm_onset)$coefficients[2, 4], digits = 3), "\n")
  576. ```
  578. # Final Model on Filtered Cohort
  580. ## Data Preparation
  582. ```{r prepare_final_model}
  583. # Calculate baseline values
  584. df_clean <- df_clean[order(df_clean$Patient_ID, df_clean$Months_since_bl), ]
  586. Age_bl <- c()
  587. MMSE_bl <- c()
  588. DX_bl <- c()
  589. AV45_bl <- c()
  590. PIB_bl <- c()
  591. CSF_AB_bl <- c()
  592. prograte <- c()
  593. visits <- c()
  595. for (i in unique(df_clean$Patient_ID)) {
  596. d <- subset(df_clean, Patient_ID == i)
  597. datapoints <- nrow(d)
  598. visits <- c(visits, rep(datapoints, datapoints))
  600. # Calculate progression rate
  601. if (datapoints >= 2) {
  602. prograte <- c(prograte, rep(coef(lm(MMSE ~ Months_since_bl, data = d))[2], datapoints))
  603. } else {
  604. prograte <- c(prograte, rep(NA, datapoints))
  605. }
  607. Age_bl <- c(Age_bl, rep(min(d$Age), datapoints))
  608. MMSE_bl <- c(MMSE_bl, rep(d$MMSE[1], datapoints))
  609. DX_bl <- c(DX_bl, rep(d$DX[1], datapoints))
  610. AV45_bl <- c(AV45_bl, rep(d$AV45[1], datapoints))
  611. PIB_bl <- c(PIB_bl, rep(d$PIB[1], datapoints))
  612. CSF_AB_bl <- c(CSF_AB_bl, rep(d$CSF_AB[1], datapoints))
  613. }
  615. df_clean$Age_bl <- Age_bl
  616. df_clean$MMSE_bl <- MMSE_bl
  617. df_clean$DX_bl <- DX_bl
  618. df_clean$AV45_bl_calc <- AV45_bl
  619. df_clean$PIB_bl_calc <- PIB_bl
  620. df_clean$CSF_AB_bl_calc <- CSF_AB_bl
  621. df_clean$prograte <- prograte
  622. df_clean$visits <- visits
  624. # Remove old baseline columns
  625. df_clean$MMSE.bl <- NULL
  626. df_clean$CSF_AB.bl <- NULL
  627. df_clean$PIB.bl <- NULL
  628. df_clean$AV45.bl <- NULL
  629. df_clean$DX.bl <- NULL
  631. # Convert variables
  632. df_clean$Months_since_bl <- as.numeric(as.character(df_clean$Months_since_bl))
  633. df_clean$Gender <- as.factor(df_clean$Gender)
  634. df_clean$PTID <- as.character(df_clean$Patient_ID)
  636. # Merge with PCA data
  637. pcs <- read.table("./QCjune/adnimerge_genomes/FINALADNI_nohapmap_pca.eigenvec")[, 2:4]
  638. colnames(pcs) <- c("PTID", "PC1", "PC2")
  639. df_final <- inner_join(df_clean, pcs, by = "PTID", multiple = "all")
  641. cat("Patients with clinical and genetic data:", length(unique(df_final$Patient_ID)), "\n")
  642. cat("Final datapoints for modeling:", nrow(df_final), "\n")
  643. ```
  645. ## Optimal Model: Random Slopes and Intercepts
  647. Based on the model comparison in @step1, we fit Model 9 to the filtered cohort:
  649. ```{r final_model}
  650. final_model <- lmer(MMSE ~ MMSE_bl + Age_bl + Months_since_bl + (Months_since_bl|Patient_ID) +
  651. Education + Gender + PC1 + PC2,
  652. data = df_final, REML = TRUE)
  654. summary(final_model)
  655. ```
  657. ## Model Diagnostics
  659. ```{r model_diagnostics, fig.height=8}
  660. # Shapiro-Wilk test for normality of residuals
  661. shapiro_result <- shapiro.test(resid(final_model))
  662. cat("Shapiro-Wilk test for normality:\n")
  663. cat(" W =", format(shapiro_result$statistic, digits = 4), "\n")
  664. cat(" p-value =", format(shapiro_result$p.value, digits = 4), "\n\n")
  666. # Homogeneity of variances test
  667. hom_var_test <- function(m) {
  668. p <- coef(summary(lm(abs(resid(m)) ~ fitted(m))))[2, 4]
  669. return(p)
  670. }
  672. hom_p <- hom_var_test(final_model)
  673. cat("Homogeneity of variance test:\n")
  674. cat(" p-value =", format(hom_p, digits = 4), "\n")
  675. cat(" ", ifelse(hom_p < 0.05, "❌ Assumption violated", "✓ Assumption met"), "\n\n")
  677. # Model performance
  678. model_perf <- performance::r2(final_model)
  679. cat("Model Performance:\n")
  680. cat(" R² Conditional:", format(model_perf$R2_conditional, digits = 3), "\n")
  681. cat(" R² Marginal:", format(model_perf$R2_marginal, digits = 3), "\n")
  682. cat(" AIC:", format(AIC(final_model), digits = 1), "\n")
  683. cat(" BIC:", format(BIC(final_model), digits = 1), "\n")
  684. ```
  686. ### Random Effects Visualization
  688. ```{r random_effects, fig.height=6}
  689. # Extract random effects
  690. ranef_data <- as.data.frame(ranef(final_model)$Patient_ID)
  691. ranef_data$Patient_ID <- rownames(ranef_data)
  692. colnames(ranef_data)[1:2] <- c("Intercept", "Slope")
  694. # Merge with baseline diagnosis
  695. group_info <- df_final %>%
  696. dplyr::select(Patient_ID, DX_bl) %>%
  697. dplyr::distinct()
  699. ranef_with_group <- merge(ranef_data, group_info, by = "Patient_ID")
  701. # Plot random slopes by baseline diagnosis
  702. ggplot(ranef_with_group, aes(x = Slope, fill = DX_bl)) +
  703. geom_density(alpha = 0.6) +
  704. geom_vline(xintercept = 0, linetype = "dashed", color = "black", size = 1) +
  705. labs(title = "Distribution of Individual Progression Rates by Baseline Diagnosis",
  706. x = "Random Slope (MMSE change per month)",
  707. y = "Density",
  708. fill = "Baseline DX") +
  709. scale_fill_brewer(palette = "Set2") +
  710. theme_minimal() +
  711. theme(text = element_text(size = 12))
  712. ```
  714. # Alternative Model Specifications
  716. ## Testing MMSE_bl vs DX_bl as Baseline Predictor
  718. ```{r baseline_predictor_comparison, results='hide'}
  719. # Model with MMSE_bl (may have convergence issues)
  720. model_mmse <- glmmTMB(MMSE ~ MMSE_bl + Age_bl + Months_since_bl * APOE4 +
  721. (1 + Months_since_bl|Patient_ID) + Education + Gender + PC1 + PC2,
  722. data = df_final, REML = TRUE)
  724. # Model with MMSE_bl and BFGS optimizer
  725. model_mmse_bfgs <- glmmTMB(MMSE ~ MMSE_bl + Age_bl + Months_since_bl * APOE4 +
  726. (1 + Months_since_bl|Patient_ID) + Education + Gender + PC1 + PC2,
  727. data = df_final, REML = TRUE,
  728. control = glmmTMBControl(optimizer = optim, optArgs = list(method = "BFGS")))
  730. # Model with DX_bl
  731. model_dx <- glmmTMB(MMSE ~ DX_bl + Age_bl + Months_since_bl * APOE4 +
  732. (1 + Months_since_bl|Patient_ID) + Education + Gender + PC1 + PC2,
  733. data = df_final, REML = TRUE)
  735. # Model with DX_bl and BFGS optimizer
  736. model_dx_bfgs <- glmmTMB(MMSE ~ DX_bl + Age_bl + Months_since_bl * APOE4 +
  737. (1 + Months_since_bl|Patient_ID) + Education + Gender + PC1 + PC2,
  738. data = df_final, REML = TRUE,
  739. control = glmmTMBControl(optimizer = optim, optArgs = list(method = "BFGS")))
  740. ```
  742. ```{r baseline_comparison}
  743. models_baseline <- list(model_mmse, model_mmse_bfgs, model_dx, model_dx_bfgs)
  744. model_names_baseline <- c("MMSE_bl (default)", "MMSE_bl (BFGS)", "DX_bl (default)", "DX_bl (BFGS)")
  746. # Performance metrics
  747. aics <- sapply(models_baseline, AIC)
  748. r2cond <- sapply(models_baseline, function(m) r2(m)$R2_conditional)
  749. r2marg <- sapply(models_baseline, function(m) r2(m)$R2_marginal)
  750. shap_w <- sapply(models_baseline, function(m) shapiro.test(resid(m))$statistic)
  751. shap_p <- sapply(models_baseline, function(m) shapiro.test(resid(m))$p.value)
  752. homosced <- sapply(models_baseline, function(m) hom_var_test(m))
  754. perform_baseline <- data.frame(
  755. Model = model_names_baseline,
  756. AIC = aics,
  757. R2_Conditional = r2cond,
  758. R2_Marginal = r2marg,
  759. Shapiro_W = shap_w,
  760. Shapiro_p = shap_p,
  761. Homoscedasticity_p = homosced
  762. )
  764. kable(perform_baseline, caption = "Model Comparison: MMSE_bl vs DX_bl as Baseline Predictor", digits = 2)
  765. ```
  767. # Polygenic Score Analysis
  769. ## Polygenic Hazard Score (PHS)
  771. ```{r phs_analysis}
  772. # Read PHS data
  773. PHS <- read.csv("./QCjune/clinicaladni/DESIKANLAB_18Jun2024.csv")
  774. df_final$RID <- as.integer(df_final$RID)
  775. df_PHS <- inner_join(df_final, PHS, by = "RID", multiple = "all")
  777. cat("Patients with PHS data:", length(unique(df_PHS$Patient_ID)), "\n")
  778. cat("Datapoints:", nrow(df_PHS), "\n")
  780. # Model with PHS
  781. model_PHS <- glmmTMB(MMSE ~ DX_bl + Age_bl + PHS * Months_since_bl +
  782. (1 + Months_since_bl|Patient_ID) + Education + Gender + PC1 + PC2,
  783. data = df_PHS, REML = TRUE)
  785. # Extract coefficients
  786. phs_coef <- summary(model_PHS)$coef$cond
  787. kable(phs_coef, caption = "PHS Model Coefficients", digits = 4)
  788. ```
  790. ## Polygenic Risk Scores (PRS)
  792. ```{r prs_analysis}
  793. # Read PRS data
  794. PRS <- read.csv("./QCjune/clinicaladni/PRS_andrealtmann.csv")[, 1:5]
  795. df_PRS <- inner_join(df_final, PRS, by = "RID", multiple = "all")
  797. cat("Patients with PRS data:", length(unique(df_PRS$Patient_ID)), "\n")
  798. cat("Datapoints:", nrow(df_PRS), "\n")
  800. # Models with different PRS
  801. model_PRS1 <- glmmTMB(MMSE ~ DX_bl + Age_bl + PRS1 * Months_since_bl +
  802. (1 + Months_since_bl|Patient_ID) + Education + Gender + PC1 + PC2,
  803. data = df_PRS, REML = TRUE)
  805. model_PRS2 <- glmmTMB(MMSE ~ DX_bl + Age_bl + PRS2 * Months_since_bl +
  806. (1 + Months_since_bl|Patient_ID) + Education + Gender + PC1 + PC2,
  807. data = df_PRS, REML = TRUE)
  809. model_PRScs <- glmmTMB(MMSE ~ DX_bl + Age_bl + PRScs_auto * Months_since_bl +
  810. (1 + Months_since_bl|Patient_ID) + Education + Gender + PC1 + PC2,
  811. data = df_PRS, REML = TRUE)
  813. # Summary table
  814. prs_results <- data.frame(
  815. Model = c("PRS1", "PRS2", "PRScs_auto"),
  816. AIC = c(AIC(model_PRS1), AIC(model_PRS2), AIC(model_PRScs)),
  817. PRS_Interaction_Beta = c(
  818. summary(model_PRS1)$coef$cond["PRS1:Months_since_bl", "Estimate"],
  819. summary(model_PRS2)$coef$cond["PRS2:Months_since_bl", "Estimate"],
  820. summary(model_PRScs)$coef$cond["PRScs_auto:Months_since_bl", "Estimate"]
  821. ),
  822. PRS_Interaction_p = c(
  823. summary(model_PRS1)$coef$cond["PRS1:Months_since_bl", "Pr(>|z|)"],
  824. summary(model_PRS2)$coef$cond["PRS2:Months_since_bl", "Pr(>|z|)"],
  825. summary(model_PRScs)$coef$cond["PRScs_auto:Months_since_bl", "Pr(>|z|)"]
  826. )
  827. )
  829. kable(prs_results, caption = "Polygenic Risk Score Models - Interaction Effects", digits = 4)
  830. ```
  832. # Advanced Model Testing
  834. ## Data Normalization
  836. Testing whether data transformation improves model fit:
  838. ```{r normalization}
  839. # Reflect and transform MMSE_bl (beta-distributed)
  840. MMSE_bl_refl <- max(df_final$MMSE_bl + 1) - df_final$MMSE_bl
  841. MMSE_bl_norm <- sqrt(MMSE_bl_refl)
  842. df_final$MMSE_bl_norm <- MMSE_bl_norm
  844. # Reflect and transform MMSE
  845. MMSE_refl <- max(df_final$MMSE + 1) - df_final$MMSE
  846. MMSE_norm <- sqrt(MMSE_refl)
  847. df_final$MMSE_norm <- MMSE_norm
  849. # Normalized model
  850. model_norm <- lmer(MMSE_norm ~ MMSE_bl_norm + Age_bl + Months_since_bl +
  851. (Months_since_bl|Patient_ID) + Education + Gender + PC1 + PC2,
  852. data = df_final, REML = TRUE)
  853. ```
  855. ## Data Scaling
  857. ```{r scaling}
  858. # Scale numeric variables
  859. df_scaled <- df_final
  860. numeric_cols <- sapply(df_scaled, is.numeric)
  862. for (col in names(df_scaled)[numeric_cols]) {
  863. if (col != "Patient_ID" && col != "RID") {
  864. df_scaled[[col]] <- round((df_scaled[[col]] - mean(df_scaled[[col]], na.rm = TRUE)) /
  865. sd(df_scaled[[col]], na.rm = TRUE), 3)
  866. }
  867. }
  869. # Scaled model
  870. model_scaled <- lmer(MMSE ~ MMSE_bl + Age_bl + Months_since_bl +
  871. (Months_since_bl|Patient_ID) + Education + Gender + PC1 + PC2,
  872. data = df_scaled, REML = TRUE)
  874. # Normalized and scaled model
  875. model_norm_scaled <- lmer(MMSE_norm ~ MMSE_bl_norm + Age_bl + Months_since_bl +
  876. (Months_since_bl|Patient_ID) + Education + Gender + PC1 + PC2,
  877. data = df_scaled, REML = TRUE)
  878. ```
  880. ## Comprehensive Model Comparison
  882. ```{r model_comparison_all}
  883. models_all <- list(final_model, model_scaled, model_norm, model_norm_scaled)
  884. model_names_all <- c("Original", "Scaled", "Normalized", "Normalized + Scaled")
  886. aics <- sapply(models_all, AIC)
  887. bics <- sapply(models_all, BIC)
  888. r2cond <- sapply(models_all, function(m) r2(m)$R2_conditional)
  889. r2marg <- sapply(models_all, function(m) r2(m)$R2_marginal)
  890. shap_w <- sapply(models_all, function(m) shapiro.test(resid(m))$statistic)
  891. shap_p <- sapply(models_all, function(m) shapiro.test(resid(m))$p.value)
  892. homosced <- sapply(models_all, function(m) hom_var_test(m))
  894. comparison_table <- data.frame(
  895. Model = model_names_all,
  896. AIC = aics,
  897. BIC = bics,
  898. R2_Conditional = r2cond,
  899. R2_Marginal = r2marg,
  900. Shapiro_W = shap_w,
  901. Shapiro_p = shap_p,
  902. Homoscedasticity_p = homosced
  903. )
  905. kable(comparison_table,
  906. caption = "Comprehensive Model Comparison: Transformation and Scaling Effects",
  907. digits = 4)
  908. ```
  910. ::: {.callout-tip}
  911. ## Model Selection Summary
  913. Lower AIC/BIC values indicate better model fit. Higher R² values indicate more variance explained. The Shapiro-Wilk test assesses normality of residuals (p > 0.05 desired). Homoscedasticity test assesses equal variance (p > 0.05 desired).
  914. :::
  916. # Patient Summary Table
  918. ## Generate Comprehensive Summary
  920. ```{r patient_summary_table}
  921. # Read FAM file for genotype status
  922. fam_data <- read.table("./QCjune/adnimerge_genomes/adnimerge.fam", header = FALSE)
  923. fam_ptids <- as.character(fam_data[, 2])
  925. # Get list of filtered patients
  926. filtered_patients <- unique(df_clean$Patient_ID)
  928. # Create summary for ALL patients in original dataset
  929. unique_patients <- unique(dataset$Patient_ID)
  931. summary_list <- list()
  933. for (patient in unique_patients) {
  934. patient_data <- dataset[dataset$Patient_ID == patient, ]
  935. patient_data <- patient_data[order(patient_data$Months_since_bl), ]
  937. # Basic info
  938. ptid <- as.character(patient)
  939. n_obs <- nrow(patient_data)
  940. n_mmse_obs <- sum(!is.na(patient_data$MMSE))
  942. # Baseline diagnosis
  943. dx_baseline <- as.character(patient_data$DX.bl[1])
  944. if (is.na(dx_baseline) || dx_baseline == "") dx_baseline <- "NA"
  946. # Last diagnosis
  947. max_idx <- which.max(patient_data$Months_since_bl)
  948. dx_last <- as.character(patient_data$DX[max_idx])
  949. if (is.na(dx_last) || dx_last == "") dx_last <- "NA"
  951. # Last MMSE
  952. last_mmse <- patient_data$MMSE[max_idx]
  953. if (is.na(last_mmse)) last_mmse <- "NA" else last_mmse <- as.character(last_mmse)
  955. # In FAM file
  956. in_fam <- ifelse(any(grepl(ptid, fam_ptids, fixed = TRUE)), "Yes", "No")
  958. # Amyloid-beta status
  959. av45_bl <- patient_data$AV45.bl[1]
  960. pib_bl <- patient_data$PIB.bl[1]
  961. csf_ab_bl <- patient_data$CSF_AB.bl[1]
  963. abeta_pos <- FALSE
  964. if (!is.na(av45_bl) && av45_bl >= th_amyl_AV45) abeta_pos <- TRUE
  965. if (!is.na(pib_bl) && pib_bl >= th_amyl_PIB) abeta_pos <- TRUE
  966. if (!is.na(csf_ab_bl) && csf_ab_bl <= th_amyl_CSF) abeta_pos <- TRUE
  967. abeta_status <- ifelse(abeta_pos, "AbetaPos", "AbetaNeg")
  969. # DX trajectory
  970. dx_trajectory <- sapply(patient_data$DX, function(x) {
  971. if (is.na(x) || x == "") return("NA")
  972. return(as.character(x))
  973. })
  974. dx_trajectory_str <- paste(dx_trajectory, collapse = "_")
  976. # Always CN
  977. dx_values_no_na <- patient_data$DX[!is.na(patient_data$DX) & patient_data$DX != ""]
  978. always_cn <- ""
  979. if (length(dx_values_no_na) > 0 && all(dx_values_no_na == "CN")) {
  980. always_cn <- "alwaysCN"
  981. }
  983. # Survives filtering
  984. survives <- ifelse(ptid %in% filtered_patients, "Yes", "No")
  986. summary_list[[patient]] <- data.frame(
  987. PTID = ptid,
  988. N_Observations = n_obs,
  989. N_MMSE_Observations = n_mmse_obs,
  990. DX_Baseline = dx_baseline,
  991. DX_Last = dx_last,
  992. Last_MMSE = last_mmse,
  993. In_FAM_File = in_fam,
  994. AbetaPos = abeta_status,
  995. DX_Trajectory = dx_trajectory_str,
  996. AlwaysCN = always_cn,
  997. Survives_Filtering = survives,
  998. stringsAsFactors = FALSE
  999. )
  1000. }
  1002. # Combine into data frame
  1003. summary_table <- bind_rows(summary_list)
  1004. summary_table <- summary_table[order(summary_table$PTID), ]
  1006. # Save table
  1007. write.csv(summary_table, "./patient_summary_table_complete.csv", row.names = FALSE, quote = TRUE)
  1009. cat("Total patients in summary table:", nrow(summary_table), "\n")
  1010. cat("Patients surviving filtering:", sum(summary_table$Survives_Filtering == "Yes"), "\n")
  1011. cat("Patients not surviving filtering:", sum(summary_table$Survives_Filtering == "No"), "\n")
  1012. ```
  1014. ## Summary Table Preview
  1016. ```{r summary_preview}
  1017. kable(head(summary_table, 20),
  1018. caption = "Patient Summary Table (First 20 Patients)",
  1019. digits = 2)
  1020. ```
  1022. ## Filtering Statistics
  1024. ```{r filtering_stats}
  1025. filter_stats <- data.frame(
  1026. Category = c("Total Patients", "Amyloid-Beta Positive", "After MMSE Filtering",
  1027. "After Clinical Filtering", "With Genotype Data", "Final Cohort"),
  1028. Count = c(
  1029. length(unique(dataset$Patient_ID)),
  1030. length(unique(df_amyloid$Patient_ID)),
  1031. length(unique(df$Patient_ID)),
  1032. length(unique(df_clean$Patient_ID)),
  1033. length(unique(df_final$Patient_ID)),
  1034. length(unique(df_final$Patient_ID))
  1035. ),
  1036. Percentage = c(
  1037. 100,
  1038. 100 * length(unique(df_amyloid$Patient_ID)) / length(unique(dataset$Patient_ID)),
  1039. 100 * length(unique(df$Patient_ID)) / length(unique(dataset$Patient_ID)),
  1040. 100 * length(unique(df_clean$Patient_ID)) / length(unique(dataset$Patient_ID)),
  1041. 100 * length(unique(df_final$Patient_ID)) / length(unique(dataset$Patient_ID)),
  1042. 100 * length(unique(df_final$Patient_ID)) / length(unique(dataset$Patient_ID))
  1043. )
  1044. )
  1046. kable(filter_stats, caption = "Filtering Pipeline Statistics", digits = 2)
  1047. ```
  1049. # Conclusions
  1051. This comprehensive analysis identified **`r length(unique(df_final$Patient_ID))` patients** with:
  1053. - Confirmed amyloid-beta pathology
  1054. - Sufficient longitudinal cognitive data (3-5 observations)
  1055. - Evidence of cognitive decline
  1056. - High-quality genotype data
  1058. The optimal model (random slopes and intercepts with full covariates) achieved:
  1060. - **R² Conditional**: `r format(r2(final_model)$R2_conditional, digits = 3)`
  1061. - **R² Marginal**: `r format(r2(final_model)$R2_marginal, digits = 3)`
  1062. - **AIC**: `r format(AIC(final_model), digits = 1)`
  1064. Polygenic scores (PHS and PRS) were successfully integrated and show no effect on progression.
  1066. # Session Information
  1068. ```{r session_info}
  1069. sessionInfo()
  1070. ```

comprehensive_analysis.qmd at commit 0bd2f20, no license · at the source

Overview

Authors: Celeste E. Cohen1,2, Shane Fernandez3,4, Umran Yaman1,5, Ahmad R. Ehyaei6, Eleftheria Kodosaki1,5, Aydan Askarova7,8, Tenielle Porter3,4, Eleanor O’Brien3,4, Australian Imaging Biomarkers and Lifestyle Study, Alzheimer’s Disease Neuroimaging Initiative, Paul Maruff9,10, Alexi Nott7,8, John A. Hardy1,5,11, Simon M. Laws3,4, Dervis A. Salih1,5, Maryam Shoai1,5
ORCID iDs: Umran Yaman, Eleftheria Kodosaki, John A. Hardy, Maryam Shoai
  1. Department of Neurodegenerative Disease, UCL Queen Square Institute of Neurology,London, UK
  2. Wellcome Sanger Institute,Cambridge, UK
  3. Centre for Precision Health, Edith Cowan University,Joondalup, WA Australia
  4. Collaborative Genomics and Translation Group, School of Medical and Health Sciences, Edith Cowan University,Joondalup, WA Australia
  5. Dementia Research Institute, University College London,London, UK
  6. Max-Planck-Institute for Intelligent Systems,Tubingen, Germany
  7. UK Dementia Research Institute, Imperial College London,London, UK
  8. Department of Brain Sciences, Imperial College London,London, UK
  9. Florey Department of Neuroscience and Mental Health, University of Melbourne,Melbourne, Australia
  10. CogState Ltd,Melbourne, Australia
  11. Reta Lila Weston Institute of Neurological Studies,London, UK
Institutions: UCL Queen Square Institute of Neurology (United Kingdom); Wellcome Sanger Institute (United Kingdom); Edith Cowan University (Australia); University College London (United Kingdom); Max Planck Institute for Intelligent Systems (Germany); Imperial College London (United Kingdom); The University of Melbourne (Australia)
Journal: Alzheimer's research & therapy, volume 18, issue 1, article 142
Dates: received 3 November 2025; accepted 24 March 2026; published online 25 April 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1186/s13195-026-02036-1 · PMID 42035106 · PMCID PMC13248444 · OpenAlex W4416378515
Open access: gold, a free copy (OpenAlex)
Preprint: doi:10.1101/2025.11.17.25340247
Status: code verified
Categories: genetics / omics (modality), human (organism), Alzheimer's / dementia (population), clinical / translational (subfield)
Methods: Preprocessing, Statistics, Smoothing, state filtering, decompositions
Keywords: Progression, Alzheimer’s disease, Genome-wide association studies, APOE‑ε4, Polygenic risk score
MeSH: Alzheimer Disease*, Disease Progression*, Genetic Predisposition to Disease*, Aged, Aged, 80 and over, Apolipoprotein E4, Female, Genetic Risk Score, Genome-Wide Association Study, Humans, Longitudinal Studies, Male, Polymorphism, Single Nucleotide, Risk Factors (* major topic)
Topic: Alzheimer's disease research and treatments (Physiology, Medicine), according to OpenAlex
Funding: Imperial College London (President’s PhD Scholarships award); UK Dementia Research Institute through UK DRI Ltd (UKDRI-5208); Vivensa Foundation (AISRPG2305\26); Alzheimer’s Association (ADSF-24-1345198-C); Fidelity Foundation; Dolby Family Ventures; National Health and Medical Research Council (GNT1161706, GNT1191535); Fidelity Foundation, United States
Citations: not cited yet (Europe PMC); 103 references in the paper

Abstract

Background: Recent trials in Alzheimer’s disease (AD) demonstrate encouraging outcomes. These trials target risk mechanisms identified through genetic analysis whilst directly aiming to reduce progression rates. Evidence from other neurodegenerative diseases suggests the genetics of progression is distinct from risk of disease. To expand these initial successes and improve clinical outcomes further we need to understand genetics of progression of disease. These can be deduced through rigorous analysis of meticulously phenotyped longitudinal cohorts. In this study we first looked at known genetic drivers of risk, namely polygenic risk scores for AD and APOE‑ε4, to assess their role in progression. This was then extended to a genome wide association analysis to identify the role of other genetic variants in progression of AD.

Methods: A total of 387 individuals with genetic data, amyloid positivity, and in active decline (ADNI (n = 222) and AIBL(n = 165)) were used to perform generalised mixed effects linear model genome wide association studies of longitudinal cognitive decline as measured by mini mental state examination (MMSE). The resulting summary statistics were subjected to functional annotation, and colocalisation analyses.

Results: Established AD risk factors, including APOE‑ε4 dosage and polygenic risk scores, were not associated with disease progression in amyloid positive individuals who are actively declining. A mixed effects GWAS meta-analysis revealed one genome-wide significant locus on chromosome 22 (rs78369883) and several nominally significant loci linked with AD progression. Functional annotation, finemapping, and colocalisation analyses implicated genes primarily involved in immune response, neurodegeneration (including tau pathology), brain resilience, and neurogenesis. These progression-related genes were significantly enriched in neuronal-interferon-microglial signalling pathways and normal homeostatic processes of neuronal networks, with specific enrichment in dopaminergic and inhibitory neuronal populations.

Conclusion: These findings enhance our understanding of the biological underpinnings of AD progression, opening new avenues for therapeutic intervention.

Reproduced under the paper's license (CC BY), from the paper cited above.

Repository

Its files are read in the Code ↔ Paper reader above, with 3 matches between paragraphs and lines of code.

MaryamShoai/Codes

License: none: the authors keep all their rights
State: the link answers, verified on 30 September 2026
Evidence: files inventoried
Commit: 0bd2f20200b3fcd0d8403aada2da72be3393bb58, 16 October 2025
Languages: R (5), Quarto (1)
Size: 19 files, 6 scripts
Software Heritage: not archived
Found in: “Data availability”
Holds: README, 1 notebook
Not found: license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: tidyverse (4 files), data.table (3 files), ggplot2 (3 files), glmmTMB (2 files), cowplot (1 file), easystats (1 file), emmeans (1 file), lme4 (1 file), lmerTest (1 file), metafor (1 file), nlme (1 file), reticulate (1 file)
Availability: 1 check, the latest on 30 September 2026: the link answers
  • 30 September 2026: the link answers
7 files
At the source: github.com/MaryamShoai/Codes

The paper's code and data availability statement is in the Data section.

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

  • 1 repository of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 6 scripts, each with its path and the digest of its content;
  • 3 matches between paragraphs of the paper and lines of the code (method lexical-v1);
  • neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

No dataset and no data link were found in the paper.

Data availability

The codes used for running the mixed effect models and the full resulting summary statistics are available from: https://github.com/MaryamShoai/Codes The Alzheimer's Disease Neuroimaging Initiative (ADNI) and the Australian Imaging, Biomarkers and Lifestyle (AIBL) study of aging are publicly available and can be accessed upon request. ADNI data are available via the LONI Image and Data Archive (IDA LONI). AIBL data are available via expression of interest to AIBL data access committee.

Reproduced under the paper's license (CC BY), from the paper cited above.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 1, 30 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 16 authors, 5 keywords, 14 MeSH terms, 8 funders, 103 references.

Cite

This paper

Cohen, C. E., Fernandez, S., Yaman, U., Ehyaei, A. R., Kodosaki, E., Askarova, A., Porter, T., O’Brien, E., Australian Imaging Biomarkers and Lifestyle Study, Alzheimer’s Disease Neuroimaging Initiative, Maruff, P., Nott, A., Hardy, J. A., Laws, S. M., A. Salih, D., & Shoai, M. (2026). Genetic drivers of progression in Alzheimer's disease are distinct from disease risk. Alzheimer's research & therapy, 18(1), 142. https://doi.org/10.1186/s13195-026-02036-1

BibTeX

@article{cohen2026genetic,
author = {Cohen, Celeste E. and Fernandez, Shane and Yaman, Umran and Ehyaei, Ahmad R. and Kodosaki, Eleftheria and Askarova, Aydan and Porter, Tenielle and O’Brien, Eleanor and {Australian Imaging Biomarkers and Lifestyle Study} and {Alzheimer’s Disease Neuroimaging Initiative} and Maruff, Paul and Nott, Alexi and Hardy, John A. and Laws, Simon M. and A. Salih, Dervis and Shoai, Maryam},
title = {{Genetic drivers of progression in Alzheimer's disease are distinct from disease risk}},
journal = {Alzheimer's research \& therapy},
year = {2026},
month = apr,
volume = {18},
number = {1},
pages = {142},
publisher = {BMC},
issn = {1758-9193},
doi = {10.1186/s13195-026-02036-1},
url = {https://doi.org/10.1186/s13195-026-02036-1},
pmid = {42035106},
pmcid = {PMC13248444}
}

RIS

TY - JOUR
AU - Cohen, Celeste E.
AU - Fernandez, Shane
AU - Yaman, Umran
AU - Ehyaei, Ahmad R.
AU - Kodosaki, Eleftheria
AU - Askarova, Aydan
AU - Porter, Tenielle
AU - O’Brien, Eleanor
AU - Australian Imaging Biomarkers and Lifestyle Study
AU - Alzheimer’s Disease Neuroimaging Initiative
AU - Maruff, Paul
AU - Nott, Alexi
AU - Hardy, John A.
AU - Laws, Simon M.
AU - A. Salih, Dervis
AU - Shoai, Maryam
TI - Genetic drivers of progression in Alzheimer's disease are distinct from disease risk
T2 - Alzheimer's research & therapy
J2 - Alzheimers Res Ther
PY - 2026
DA - 2026/04/25
VL - 18
IS - 1
SP - 142
SN - 1758-9193
PB - BMC
DO - 10.1186/s13195-026-02036-1
UR - https://doi.org/10.1186/s13195-026-02036-1
LA - en
ER -

CSL-JSON

{
"id": "10.1186/s13195-026-02036-1",
"type": "article-journal",
"title": "Genetic drivers of progression in Alzheimer's disease are distinct from disease risk",
"container-title": "Alzheimer's research & therapy",
"author": [
{
"family": "Cohen",
"given": "Celeste E."
},
{
"family": "Fernandez",
"given": "Shane"
},
{
"family": "Yaman",
"given": "Umran"
},
{
"family": "Ehyaei",
"given": "Ahmad R."
},
{
"family": "Kodosaki",
"given": "Eleftheria"
},
{
"family": "Askarova",
"given": "Aydan"
},
{
"family": "Porter",
"given": "Tenielle"
},
{
"family": "O’Brien",
"given": "Eleanor"
},
{
"literal": "Australian Imaging Biomarkers and Lifestyle Study"
},
{
"literal": "Alzheimer’s Disease Neuroimaging Initiative"
},
{
"family": "Maruff",
"given": "Paul"
},
{
"family": "Nott",
"given": "Alexi"
},
{
"family": "Hardy",
"given": "John A."
},
{
"family": "Laws",
"given": "Simon M."
},
{
"family": "A. Salih",
"given": "Dervis"
},
{
"family": "Shoai",
"given": "Maryam"
}
],
"container-title-short": "Alzheimers Res Ther",
"volume": "18",
"issue": "1",
"page": "142",
"DOI": "10.1186/s13195-026-02036-1",
"PMID": "42035106",
"PMCID": "PMC13248444",
"ISSN": "1758-9193",
"publisher": "BMC",
"URL": "https://doi.org/10.1186/s13195-026-02036-1",
"language": "en",
"issued": {
"date-parts": [
[
2026,
4,
25
]
]
}
}

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.1038/s41380-026-03686-1 [code]
Early oligodendrocyte dysfunction signature in Alzheimer's disease: Insights from DNA methylomics and transcriptomics.
Journal: Molecular psychiatry
In common: data.table, ggplot2, tidyverse, Alzheimer's / dementia, genetics / omics, 5 references, 2 authors
[2] doi:10.1038/s41562-026-02486-5 [code]
Genome-wide association studies of infant and toddler temperament in European and multi-ancestry populations.
Journal: Nature human behaviour
In common: metafor, nlme, data.table, 2 other tools, genetics / omics, 7 references
[3] doi:10.1016/j.celrep.2026.117505 [code]
Impaired spatial coding and neuronal hyperactivity in the medial entorhinal cortex of aged APP knock-in mice.
Journal: Cell reports
In common: glmmTMB, metafor, nlme, 6 other tools, Alzheimer's / dementia
[4] doi:10.1371/journal.pgen.1012170 [code]
Genome wide association study meta-analysis of neuropathologic lesions of Alzheimer's disease and related dementias in a multi-site autopsy cohort.
Journal: PLoS genetics
In common: data.table, ggplot2, tidyverse, Alzheimer's / dementia, clinical / translational, genetics / omics, 8 references
[5] doi:10.1002/hbm.70605 [code]
BrainEnrich: Revealing Biological Insights for Imaging-Derived Phenotypes Through Transcriptomic Enrichment.
Journal: Human brain mapping
In common: nlme, reticulate, easystats, 5 other tools, genetics / omics, 1 reference
[6] doi:10.1038/s41467-026-71682-8 [code]
GWAS meta-analysis of cerebrospinal fluid Alzheimer's biomarkers reveals loci regulating lipids, brain volume and autophagy.
Journal: Nature communications
In common: data.table, ggplot2, tidyverse, Alzheimer's / dementia, genetics / omics, 7 references
[7] doi:10.1038/s41467-026-73865-9 [code]
Histamine shapes the neurocomputational dynamics of human learning.
Journal: Nature communications
In common: metafor, easystats, emmeans, 6 other tools
[8] doi:10.1073/pnas.2606871123 [code]
Oxytocin modulates the neurocomputational mechanisms engaged in learning rank relationships in social networks.
Journal: Proceedings of the National Academy of Sciences of the United States of America
In common: nlme, easystats, emmeans, 6 other tools
[9] doi:10.1186/s40168-026-02342-8 [code]
Impacts of host genetics on gut microbiome composition in Alzheimer's disease.
Journal: Microbiome
In common: metafor, cowplot, data.table, 2 other tools, Alzheimer's / dementia, genetics / omics, 4 references
[10] doi:10.1038/s41467-026-74753-y [code]
A human-specific microRNA controls the timing of excitatory synaptogenesis.
Journal: Nature communications
In common: reticulate, easystats, emmeans, 6 other tools

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.