Prediction of mild cognitive impairment progression using time-sensitive multimodal biomarkers.

2. ############### Script generated by Jonathan Gallego Rudolf ####################
5. ##### Progression to MCI analysis
8. ##### 1) Cox regression models
10. ### Survival analysis based on Cox regression proportional hazard models
11. ### to estimate the risk of progression to MCI and the added value of incorporating
12. ### neurophysiological, brain structure, plasma and PET biomarkers to clinical information
15. ### The script is used to fit Cox regression model and obtain survival curves
16. ### for each model. Likelihood ratio tests are used to assess whether
17. ### the model performed significantly better than a null model including no
18. ### features, and also compared to the reference clinical model
21. ### The following models were defined
22. ### Model 1: Clinical (age, sex, education, apoe)
23. ### Model 2: Clinical + MEG (initially relative alpha power, but was extended to other frequency bands)
24. ### Model 3: Clinical + MRI (hippocampal volume norm to TIV)
25. ### Model 4: Clinical + plasma (Ab42/40 ratio, p-tau217)
26. ### Model 5: Clinical + Ab PET uptake (neocortical)
27. ### Model 6: Clinical + Tau PET uptake (entorhinal)
28. ### Model 7: All features
31. ### Additional model comparisons are performed by calculating the difference in AIC
32. ### and the concordance index of the models
34. ### Risk scores estimated for each participants are used to stratify individuals
35. ### into low-, medium- and high-risk groups, to visualize whether the survival curves
36. ### changes as a function of a given biomarker
38. ### Linear models and independent t-test were implemented to assess the contribution
39. ### of different biomarkers by estimating their correlation with the risk scores
40. ### estimated derived from its corresponding model
42. ### Run models including time interactions to assess whether the hazard risk
43. ### associated with a given biomarker changes over time. Plot the change in
44. ### hazard ratio over time, calculating the CIs while accounting for the clinical covariates
45. ### Plot the change in the linear predictor (calculated for each participant) as a function of time, to show
46. ### the time interaction effect on the association between a biomarker and the risk of progression
48. ### Compare whether adding the neurophysiological activity time-varying effects
49. ### improved the performance of the models combining clinical information and plasma/PET biomarkers
51. ### Step AIC model feature selection regression analysis was performed on the model
52. ### including all biomarkers to select a parsimonious model including a set of features that would
53. ### maximize model performance while accounting for model complexity based on AIC values
55. #### XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX ###
61. #### XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX ###
63. ### 1. Prepare data and environment----
66. ##### Generate data_meg_pet_cogn spreadsheet (n=104) using the prepare_data script
69. ### Load required packages
70. library(scico)
71. library(survival)
72. library(survminer)
73. require(ggplot2)
74. require(nlme)
75. require(MuMIn)
76. library(lubridate)
77. library(ggsurvfit)
78. library(gtsummary)
79. library(tidycmprsk)
80. library(condsurv)
81. library(dplyr)
82. library(ggpubr)
83. library(see)
84. library(scales)
85. library(tidyr)
88. ### Set path to save data
89. path="C:/Users/jogar/Documents/PREVENT_AD/Multimodal_data/Figures/MCI_review"
91. ### Load Dk atlas labels
92. dk_labels=readLines('C:/Users/jogar/Documents/PREVENT_AD/Multimodal_data/PET/Nov2023/labels_dk.csv')
94. #### XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX ###
100. #### XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX ###
102. ### 2. Get dates and calculate time diff, scale variables, add lobes (EXTRA)----
104. ### Load dates for the scanning visits, plasma measures and MCI diagnosis
105. scan_dates=read.csv("C:/Users/jogar/Documents/PREVENT_AD/Multimodal_data/Spreadsheets/multimodal_scan_dates.csv",header=T)
107. ### Change format to dates (needed for computing time diff)
108. scan_dates=scan_dates[,c(seq(1,9))]
109. scan_dates$MEG_date=as.Date(as.character(scan_dates$MEG_date), format = "%Y%m%d")
110. scan_dates$MRI_date=as.Date(as.character(scan_dates$MRI_date), format = "%Y%m%d")
111. scan_dates$Ab_date=as.Date(as.character(scan_dates$Ab_date), format = "%Y%m%d")
112. scan_dates$Tau_date=as.Date(as.character(scan_dates$Tau_date), format = "%Y%m%d")
113. scan_dates$First_cogn=as.Date(as.character(scan_dates$First_cogn), format = "%Y%m%d")
114. scan_dates$MEG_cogn=as.Date(as.character(scan_dates$MEG_cogn), format = "%Y%m%d")
115. scan_dates$MCI_date=as.Date(as.character(scan_dates$MCI_date_july_24), format = "%Y%m%d")
116. scan_dates$Plasma_date=as.Date(as.character(scan_dates$Plasma_date), format = "%Y%m%d")
117. scan_dates$ids=data_meg_pet_cogn$pscid
120. ### Add dates to spreadsheet
121. data_meg_pet_cogn=cbind(data_meg_pet_cogn,scan_dates)
124. ### Generate function to compute difference in dates
125. compute_date_diff <- function(df, date_col1, date_col2, new_col_name) {
126. # Ensure date columns are in POSIXct format
127. date1 <- as.POSIXct(df[[date_col1]], tz = "UTC")
128. date2 <- as.POSIXct(df[[date_col2]], tz = "UTC")
129. # Compute difference in days
130. diff_days <- as.numeric(difftime(date1, date2, units = "days"))
131. # Add to dataframe
132. df[[new_col_name]] <- diff_days
133. return(df)
134. }
138. ### Function to compute days to cutoff date (July 31, 2024)
139. fill_missing_diff_days <- function(df, date_col, diff_col, cutoff = "2024-07-31", new_col = "diff_years") {
140. # Ensure cutoff is Date
141. cutoff_date <- as.Date(cutoff)
142. # Loop through rows where diff_col is NA
143. for (i in seq_len(nrow(df))) {
144. if (is.na(df[[diff_col]][i])) {
145. date_val <- df[[date_col]][i]
146. if (!is.na(date_val)) {
147. diff_days <- as.numeric(difftime(as.POSIXct(cutoff_date), as.POSIXct(date_val, tz = "UTC"), units = "days"))
148. df[[diff_col]][i] <- diff_days
149. }
150. }
151. }
152. # Convert to years
153. df[[new_col]] <- df[[diff_col]] / 365
154. return(df)
155. }
158. ### Calculate difference in dates
159. data_meg_pet_cogn=compute_date_diff(data_meg_pet_cogn, "MCI_date", "MEG_date", "diff_days_meg_mci")
162. ### Format variables of interest
164. ### Ensure MCI status, apoe, and sex are defined as factors as factors
165. data_meg_pet_cogn$mci_status=as.factor(data_meg_pet_cogn$mci_status)
166. data_meg_pet_cogn$apoe=as.factor(data_meg_pet_cogn$apoe)
167. data_meg_pet_cogn$sex=as.factor(data_meg_pet_cogn$sex)
168. data_meg_pet_cogn$sex=factor(data_meg_pet_cogn$sex, levels = c("Male", "Female"))
174. ### Z-score biomarker measurements
175. data_meg_pet_cogn$age_z=scale(data_meg_pet_cogn$age);
176. data_meg_pet_cogn$education_z=scale(data_meg_pet_cogn$education)
177. data_meg_pet_cogn$hippocampal_volume_z=scale(data_meg_pet_cogn$hippocampal_volume)
178. data_meg_pet_cogn$tau_meta_roi_meg_delta_z=scale(data_meg_pet_cogn$tau_meta_roi_meg_delta)
179. data_meg_pet_cogn$tau_meta_roi_meg_theta_z=scale(data_meg_pet_cogn$tau_meta_roi_meg_theta)
180. data_meg_pet_cogn$tau_meta_roi_meg_alpha_z=scale(data_meg_pet_cogn$tau_meta_roi_meg_alpha)
181. data_meg_pet_cogn$tau_meta_roi_meg_beta_z=scale(data_meg_pet_cogn$tau_meta_roi_meg_beta)
182. data_meg_pet_cogn$tau_meta_roi_meg_gamma1_z=scale(data_meg_pet_cogn$tau_meta_roi_meg_gamma1)
183. data_meg_pet_cogn$tau_meta_roi_meg_gamma2_z=scale(data_meg_pet_cogn$tau_meta_roi_meg_gamma2)
185. data_meg_pet_cogn$Ab42_40_ratio_z=scale(data_meg_pet_cogn$Ab42_40_ratio)
186. data_meg_pet_cogn$ptau_217_z=scale(data_meg_pet_cogn$ptau_217)
187. data_meg_pet_cogn$amyloid_index_z=scale(data_meg_pet_cogn$amyloid_index)
188. data_meg_pet_cogn$tau_entorhinal_z=scale(data_meg_pet_cogn$tau_entorhinal)
194. ### Generate new dataset removing the subjects that do not have plasma data
195. no_plasma=c("MTL0139","MTL0224")
196. data_meg_pet_cogn3=data_meg_pet_cogn %>% filter(!pscid %in% no_plasma)
199. ### Compute additional differences in dates
200. data_meg_pet_cogn3=compute_date_diff(data_meg_pet_cogn3, "MEG_date", "Plasma_date", "diff_days_meg_plasma")
201. data_meg_pet_cogn3=compute_date_diff(data_meg_pet_cogn3, "Ab_date", "Tau_date", "diff_days_ab_tau")
202. data_meg_pet_cogn3=compute_date_diff(data_meg_pet_cogn3, "MEG_date", "MRI_date", "diff_days_meg_mri")
203. data_meg_pet_cogn3=compute_date_diff(data_meg_pet_cogn3, "Ab_date", "MEG_date", "diff_days_ab_meg")
204. data_meg_pet_cogn3=compute_date_diff(data_meg_pet_cogn3, "MCI_date", "Plasma_date", "diff_days_plasma_mci")
206. #### XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX ###
212. #### XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX ###
214. ### 3. Load extended Ab and tau PET and plasma datasets----
216. ### Load full PET dataset
217. Ab_data_all_subs=read.csv('C:/Users/jogar/Documents/PREVENT_AD/Multimodal_data/PET/Nov2023/PAD_PET_NAV_suvr_space-anat_ref-cerebellumCortex_time-4070_Nov2023.tsv',header=T, sep="")
218. Tau_data_all_subs=read.csv('C:/Users/jogar/Documents/PREVENT_AD/Multimodal_data/PET/Nov2023/PAD_PET_TAU_suvr_space-anat_ref-infcereg_time-4070_Nov2023.tsv',header=T, sep="")
220. ### Separate BL and FU
221. Ab_data_bl_all_subs <- filter(Ab_data_all_subs, Visit_label_PET == "ses-01")
222. Ab_data_fu_all_subs <- filter(Ab_data_all_subs, Visit_label_PET == "ses-02")
223. Tau_data_bl_all_subs <- filter(Tau_data_all_subs, Visit_label_PET == "ses-01")
224. Tau_data_fu_all_subs <- filter(Tau_data_all_subs, Visit_label_PET == "ses-02")
227. ### Get subset of individuals with MEG (n=104)
228. Ab_data_bl_NN_subs=filter(Ab_data_bl_all_subs,PSCID%in%data_meg_pet_cogn$pscid)
229. Ab_data_fu_NN_subs=filter(Ab_data_fu_all_subs,PSCID%in%data_meg_pet_cogn$pscid)
230. Tau_data_bl_NN_subs=filter(Tau_data_bl_all_subs,PSCID%in%data_meg_pet_cogn$pscid)
231. Tau_data_fu_NN_subs=filter(Tau_data_fu_all_subs,PSCID%in%data_meg_pet_cogn$pscid)
233. ### Keep subjects with both Ab and tau FU (N=86)
234. Ab_data_bl_both=filter(Ab_data_bl_NN_subs,PSCID%in%Tau_data_fu_NN_subs$PSCID)
235. Ab_data_fu_both=filter(Ab_data_fu_NN_subs,PSCID%in%Tau_data_fu_NN_subs$PSCID)
236. Tau_data_bl_both=filter(Tau_data_bl_NN_subs,PSCID%in%Tau_data_fu_NN_subs$PSCID)
237. Tau_data_fu_both=Tau_data_fu_NN_subs
240. ### Load demographic information
241. demographics_387=read.csv("C:/Users/jogar/Documents/PREVENT_AD/Multimodal_data/Cognition/MCI_progression/new_consensus_and_petersen_july_2024/july_2024_updated_plasma/Demographics_387.csv",header=T)
242. demographics_ab_243=filter(demographics_387,pscid%in%Ab_data_bl_all_subs$PSCID)
243. demographics_tau_240=demographics_387 %>% filter(pscid %in% Tau_data_bl_all_subs$PSCID)
245. ### Get subset of individuals with BL Ab and Tau PET
246. Ab_data_bl_243=filter(Ab_data_bl_all_subs,PSCID%in%demographics_ab_243$pscid); amyloid_index_243=Ab_data_bl_243$thresh_amyl_index
247. Tau_data_bl_240=filter(Tau_data_bl_all_subs,PSCID%in%demographics_tau_240$pscid); tau_meta_roi_240=Tau_data_bl_240$thresh_meta_roi_bilat
251. ### Load MCI dates
252. MCI_dates_all=read.csv("C:/Users/jogar/Documents/PREVENT_AD/Multimodal_data/Cognition/MCI_progression/new_consensus_and_petersen_july_2024/mci_dates_july_24_all.csv",header=T, sep=","); colnames(MCI_dates_all)=c("pscid", "mci_status")
254. ### Subset MCI list for individuals with Ab PET (n=255)
255. MCI_dates_ab_255=data.frame(pscid=Ab_data_bl_all_subs$PSCID)
256. MCI_dates_ab_255$mci_status=as.integer(Ab_data_bl_all_subs$PSCID %in% MCI_dates_all$pscid)
257. MCI_dates_ab_255=merge(MCI_dates_ab_255, MCI_dates_all, by.x = "pscid", by.y = "pscid", all.x = TRUE); colnames(MCI_dates_ab_255)=c("pscid","mci_status","mci_date")
260. ### Subset MCI list for individuals with Ab PET (n=255)
261. MCI_dates_tau_252=data.frame(pscid=Tau_data_bl_all_subs$PSCID)
262. MCI_dates_tau_252$mci_status=as.integer(Tau_data_bl_all_subs$PSCID %in% MCI_dates_all$pscid)
263. MCI_dates_tau_252=merge(MCI_dates_tau_252, MCI_dates_all, by.x = "pscid", by.y = "pscid", all.x = TRUE); colnames(MCI_dates_tau_252)=c("pscid","mci_status","mci_date")
265. MCI_dates_ab_243=filter(MCI_dates_ab_255,pscid%in%demographics_ab_243$pscid);
266. MCI_dates_tau_240=filter(MCI_dates_tau_252,pscid%in%demographics_tau_240$pscid);
270. ### Load PET scans dates
271. PET_dates_all_subs=read.csv('C:/Users/jogar/Documents/PREVENT_AD/Multimodal_data/PET/Nov2023/PAD_PET_info_Nov2023.csv',header=T, sep=",")
272. PET_dates_all_subs$date_scan=as.Date(as.character(PET_dates_all_subs$date_scan), format = "%Y%m%d")
274. ### Get dates for BL and FU Ab and tau PET scans
275. Ab_dates_bl_all_subs=filter(PET_dates_all_subs, Visit_label_PET == "ses-01" & Tracer_PET == "NAV")
276. Ab_dates_fu_all_subs=filter(PET_dates_all_subs, Visit_label_PET == "ses-02" & Tracer_PET == "NAV")
277. Tau_dates_bl_all_subs=filter(PET_dates_all_subs, Visit_label_PET == "ses-01" & Tracer_PET == "TAU")
278. Tau_dates_fu_all_subs=filter(PET_dates_all_subs, Visit_label_PET == "ses-02" & Tracer_PET == "TAU")
281. ### Get subset of individuals with MEG (used for predictive model analysis)
282. Ab_dates_bl_NN_subs=filter(Ab_dates_bl_all_subs,PSCID%in%data_meg_pet_cogn$pscid)
283. Ab_dates_fu_NN_subs=filter(Ab_dates_fu_all_subs,PSCID%in%data_meg_pet_cogn$pscid)
284. Tau_dates_bl_NN_subs=filter(Tau_dates_bl_all_subs,PSCID%in%data_meg_pet_cogn$pscid)
285. Tau_dates_fu_NN_subs=filter(Tau_dates_fu_all_subs,PSCID%in%data_meg_pet_cogn$pscid)
287. ### Keep subjects with both Ab and tau FU (N=86; used for predictive model analysis)
288. Ab_dates_bl_both=filter(Ab_dates_bl_NN_subs,PSCID%in%Tau_dates_fu_NN_subs$PSCID)
289. Ab_dates_fu_both=filter(Ab_dates_fu_NN_subs,PSCID%in%Tau_dates_fu_NN_subs$PSCID)
290. Tau_dates_bl_both=filter(Tau_dates_bl_NN_subs,PSCID%in%Tau_dates_fu_NN_subs$PSCID)
291. Tau_dates_fu_both=Tau_dates_fu_NN_subs
294. Ab_dates_bl_243=filter(Ab_dates_bl_all_subs,PSCID%in%Ab_data_bl_243$PSCID)
295. Tau_dates_bl_240=filter(Tau_dates_bl_all_subs,PSCID%in%Tau_data_bl_240$PSCID)
299. ### Create df to run Cox regression models on full Ab sample
300. df_ab_243=data.frame(cbind(pscid=Ab_data_bl_243$PSCID, age=Ab_dates_bl_243$age_PET, sex=demographics_ab_243$Sex,
301. education=demographics_ab_243$Education, apoe=demographics_ab_243$apoe,
302. amyloid_index=amyloid_index_243,PET_date=Ab_dates_bl_243$Date_PET,
303. mci_status=MCI_dates_ab_243$mci_status,
304. mci_date=MCI_dates_ab_243$mci_date))
308. df_ab_243$mci_date=as.Date(as.character(df_ab_243$mci_date), format = "%Y%m%d")
309. df_ab_243$amyloid_index=as.double(df_ab_243$amyloid_index)
310. df_ab_243$amyloid_index_z=scale(df_ab_243$amyloid_index); df_ab_243$amyloid_index_z=df_ab_243$amyloid_index_z[,1]
313. ### Create df to run Cox regression models on full Ab sample
314. df_tau_240=data.frame(cbind(pscid=Tau_data_bl_240$PSCID, age=Tau_dates_bl_240$age_PET, sex=demographics_tau_240$Sex,
315. education=demographics_tau_240$Education, apoe=demographics_tau_240$apoe,
316. tau_meta_roi=tau_meta_roi_240,PET_date=Tau_dates_bl_240$Date_PET,
317. mci_status=MCI_dates_tau_240$mci_status,
318. mci_date=MCI_dates_tau_240$mci_date))
322. df_tau_240$mci_date=as.Date(as.character(df_tau_240$mci_date), format = "%Y%m%d")
323. df_tau_240$tau_meta_roi=as.double(df_tau_240$tau_meta_roi)
324. df_tau_240$tau_meta_roi_z=scale(df_tau_240$tau_meta_roi); df_tau_240$tau_meta_roi_z=df_tau_240$tau_meta_roi_z[,1]
327. ### Calculate entorhinal cortex SUVR
328. tau_entorhinal=vector()
329. for (i in 1:nrow(Tau_data_bl_240)){
330. tau_entorhinal_s=mean(Tau_data_bl_240$ctx.lh.entorhinal[i],Tau_data_bl_240$ctx.rh.entorhinal[i])
331. tau_entorhinal=rbind(tau_entorhinal,tau_entorhinal_s)
332. }
333. tau_entorhinal=data.frame(tau_entorhinal=tau_entorhinal)
335. df_tau_240$tau_entorhinal=as.double(tau_entorhinal$tau_entorhinal)
336. df_tau_240$tau_entorhinal_z=scale(df_tau_240$tau_entorhinal); df_tau_240$tau_entorhinal_z=df_tau_240$tau_entorhinal_z[,1]
341. ### Calculate difference in dates between Ab PET and MCI diagnosis
342. df_ab_243=compute_date_diff(df_ab_243, "mci_date", "PET_date", "diff_days_pet_mci")
346. ### Get MCI date for Cox regression models for full Ab sample
347. df_ab_243$diff_days_pet_mci_cox=df_ab_243$diff_days_pet_mci
348. df_ab_243 <- fill_missing_diff_days(df_ab_243, "PET_date", "diff_days_pet_mci_cox", cutoff = "2024-07-31", new_col = "diff_years_pet_mci_cox")
351. ### Transform age and education to z scores
352. df_ab_243$age=as.double(df_ab_243$age); df_ab_243$age_z=scale(df_ab_243$age); df_ab_243$age_z=df_ab_243$age_z[,1]
353. df_ab_243$education=as.double(df_ab_243$education); df_ab_243$education_z=scale(df_ab_243$education); df_ab_243$education_z=df_ab_243$education_z[,1]
355. df_ab_243$mci_status=as.double(df_ab_243$mci_status)+1
356. df_ab_243$mci_status_cox=factor(df_ab_243$mci_status, levels = c(1, 2), labels = c(0, 1)) # Define the factor with levels 0 and 1
357. df_ab_243$mci_status_cox=as.double(df_ab_243$mci_status_cox)
358. df_ab_243$mci_status_cox=df_ab_243$mci_status_cox-1
361. ### Remove subjects that had the PET after being diagnoses as MCI
362. df_ab_226=df_ab_243[df_ab_243$diff_years_pet_mci_cox >= 0, ]
364. ### Create df with only the additional subjects
365. df_ab_124=filter(df_ab_226,!pscid%in%data_meg_pet_cogn3$pscid)
369. ### Calculate difference in dates between Tau PET and MCI diagnosis
370. df_tau_240=compute_date_diff(df_tau_240, "mci_date", "PET_date", "diff_days_pet_mci")
373. ### Compute days to cutoff date (July 31, 2024) for non-progressors to run Cox regression models
374. df_tau_240$diff_days_pet_mci_cox=df_tau_240$diff_days_pet_mci
375. df_tau_240 <- fill_missing_diff_days(df_tau_240, "PET_date", "diff_days_pet_mci_cox", cutoff = "2024-07-31", new_col = "diff_years_pet_mci_cox")
378. ### Transform age and education to z scores
379. df_tau_240$age=as.double(df_tau_240$age); df_tau_240$age_z=scale(df_tau_240$age); df_tau_240$age_z=df_tau_240$age_z[,1]
380. df_tau_240$education=as.double(df_tau_240$education); df_tau_240$education_z=scale(df_tau_240$education); df_tau_240$education_z=df_tau_240$education_z[,1]
382. ### Add MCI vars for Cox
383. df_tau_240$mci_status=as.double(df_tau_240$mci_status)+1
384. df_tau_240$mci_status_cox=factor(df_tau_240$mci_status, levels = c(1, 2), labels = c(0, 1)) # Define the factor with levels 0 and 1
385. df_tau_240$mci_status_cox=as.double(df_tau_240$mci_status_cox)
386. df_tau_240$mci_status_cox=df_tau_240$mci_status_cox-1
388. ### Remove subjects that had the PET after being diagnoses as MCI
389. df_tau_224=df_tau_240[df_tau_240$diff_years_pet_mci_cox >= 0, ]
391. ### Create df with only the additional subjects
392. df_tau_122=filter(df_tau_224,!pscid%in%data_meg_pet_cogn3$pscid)
397. ### Load full plasma data spreadsheets
399. plasma_ab4240=read.csv("C:/Users/jogar/Documents/PREVENT_AD/Multimodal_data/Plasma/PREVENT-AD_n387_IPMS_totaltau_ptau_4plex.tsv",header=T, sep="");
400. plasma_ptau_217=read.csv("C:/Users/jogar/Documents/PREVENT_AD/Multimodal_data/Plasma/PREVENT-AD_n387_ptau217_UGOT.tsv",header=T, sep="");
403. plasma_ab4240$Candidate_Age=plasma_ab4240$Candidate_Age/12; plasma_ptau_217$Candidate_Age=plasma_ptau_217$Candidate_Age/12
405. plasma_ab4240_216=plasma_ab4240[!is.na(plasma_ab4240$AB_ratio_MS_UGOT), ]
406. plasma_ptau_217_216=filter(plasma_ptau_217,PSCID %in% plasma_ab4240_216$PSCID);
408. plasma_ptau_217_216=inner_join(plasma_ptau_217_216, plasma_ab4240_216, by = c("PSCID", "Study_visit_label"))
411. demographics_plasma_216=filter(demographics_387,pscid%in%plasma_ptau_217_216$PSCID)
415. ### Subset MCI list for individuals with Ab PET (n=255)
416. MCI_dates_plasma_216=data.frame(pscid=plasma_ptau_217_216$PSCID)
417. MCI_dates_plasma_216$mci_status=as.integer(plasma_ptau_217_216$PSCID %in% MCI_dates_all$pscid)
418. MCI_dates_plasma_216=merge(MCI_dates_plasma_216, MCI_dates_all, by.x = "pscid", by.y = "pscid", all.x = TRUE); colnames(MCI_dates_plasma_216)=c("pscid","mci_status","mci_date")
422. ### Create df to run Cox regression models on full Ab sample
423. df_plasma_216=data.frame(cbind(pscid=plasma_ab4240_216$PSCID, age=plasma_ab4240_216$Candidate_Age, sex=demographics_plasma_216$Sex,
424. education=demographics_plasma_216$Education, apoe=demographics_plasma_216$apoe,
425. plasma_date=plasma_ab4240_216$Date_taken,
426. mci_status=MCI_dates_plasma_216$mci_status,
427. mci_date=MCI_dates_plasma_216$mci_date,
428. Ab42_40_ratio=plasma_ab4240_216$AB_ratio_MS_UGOT,
429. ptau_217=plasma_ptau_217_216$UGOT_ptau217))
433. df_plasma_216$mci_date=as.Date(as.character(df_plasma_216$mci_date), format = "%Y%m%d")
438. ### Calculate difference in dates between plasma and MCI diagnosis
439. df_plasma_216=compute_date_diff(df_plasma_216, "mci_date", "plasma_date", "diff_days_plasma_mci")
441. ### Compute days to cutoff date (July 31, 2024) for non-progressors to run Cox regression models
442. df_plasma_216$diff_days_plasma_mci_cox=df_plasma_216$diff_days_plasma_mci
443. df_plasma_216 <- fill_missing_diff_days(df_plasma_216, "plasma_date", "diff_days_plasma_mci_cox", cutoff = "2024-07-31", new_col = "diff_years_plasma_mci_cox")
450. ### Transform age and education to z scores
451. df_plasma_216$age=as.double(df_plasma_216$age); df_plasma_216$age_z=scale(df_plasma_216$age); df_plasma_216$age_z=df_plasma_216$age_z[,1]
452. df_plasma_216$education=as.double(df_plasma_216$education); df_plasma_216$education_z=scale(df_plasma_216$education); df_plasma_216$education_z=df_plasma_216$education_z[,1]
453. df_plasma_216$Ab42_40_ratio=as.double(df_plasma_216$Ab42_40_ratio); df_plasma_216$Ab42_40_ratio_z=scale(df_plasma_216$Ab42_40_ratio); df_plasma_216$Ab42_40_ratio_z=df_plasma_216$Ab42_40_ratio_z[,1]
454. df_plasma_216$ptau_217=as.double(df_plasma_216$ptau_217); df_plasma_216$ptau_217_z=scale(df_plasma_216$ptau_217); df_plasma_216$ptau_217_z=df_plasma_216$ptau_217_z[,1]
457. df_plasma_216$mci_status=as.double(df_plasma_216$mci_status)+1
458. df_plasma_216$mci_status_cox=factor(df_plasma_216$mci_status, levels = c(1, 2), labels = c(0, 1)) # Define the factor with levels 0 and 1
459. df_plasma_216$mci_status_cox=as.double(df_plasma_216$mci_status_cox)
460. df_plasma_216$mci_status_cox=df_plasma_216$mci_status_cox-1
463. ### Remove subjects that had the PET after being diagnoses as MCI and missing ptau217 from the same visit as the Ab4240 ratio
464. df_plasma_211=df_plasma_216[df_plasma_216$diff_years_plasma_mci_cox >= 0, ];
465. df_plasma_211=df_plasma_211[!is.na(df_plasma_211$ptau_217), ]
467. ### Create df with only the additional subjects
468. df_plasma_109=filter(df_plasma_211,!pscid%in%data_meg_pet_cogn3$pscid)
473. df_plasma_211_mci <- df_plasma_211[df_plasma_211$mci_status == 2, ]
474. df_ab_226_mci <- df_ab_226[df_ab_226$mci_status == 2, ]
475. df_tau_224_mci <- df_tau_224[df_tau_224$mci_status == 2, ]
478. mean(df_plasma_211_mci$diff_years_plasma_mci_cox); sd(df_plasma_211_mci$diff_years_plasma_mci_cox)
479. min(df_plasma_211_mci$diff_years_plasma_mci_cox); max(df_plasma_211_mci$diff_years_plasma_mci_cox);
480. mean(df_ab_226_mci$diff_years_pet_mci_cox); sd(df_ab_226_mci$diff_years_pet_mci_cox)
481. min(df_ab_226_mci$diff_years_pet_mci_cox); max(df_ab_226_mci$diff_years_pet_mci_cox);
482. mean(df_tau_224_mci$diff_years_pet_mci_cox); sd(df_tau_224_mci$diff_years_pet_mci_cox)
483. min(df_tau_224_mci$diff_years_pet_mci_cox); max(df_tau_224_mci$diff_years_pet_mci_cox);
485. #### XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX ###
491. #### XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX ###
493. ### 4. Plot MCI proportion and compare progressors vs non-progressors----
495. ### Numbers updated up to July 2024
496. mci_chart <- data.frame(mci_status = c("CU", "MCI"), num_sub = c(71, 31), perc_sub = c(69.6, 30.4))
498. ### Create a basic bar chart
499. pie_chart <- ggplot(mci_chart, aes(x = "", y = num_sub, fill = mci_status)) +
500. geom_bar(stat = "identity") + coord_polar(theta = "y") +
501. scale_fill_manual(values = c("#ACD8A7","#E0B0FF")) + theme_void()
503. ### Print the pie chart
504. print(pie_chart)
507. ### Plot variability in the time bw biomarker collection and MCI diagnosis
508. ggplot(data_meg_pet_cogn3, aes(x = mci_status, y = diff_days_meg_mci, fill = mci_status)) +
509. geom_boxplot() + geom_jitter(width = 0.2, color = "black", alpha = 0.7) + # Adds jittered points
510. scale_fill_manual(values = c("white", "#E0B0FF")) + # Example colors
511. labs(title = "Boxplot with Modern Theme", x = "Group", y = "Value") +
512. theme_minimal(base_size = 15) + theme_modern()
514. #### XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX ###
520. #### XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX ###
522. ### 5. Group comparison between MCI progressors and non-progressors----
524. ### Calculate group means and sd for MCI progressors and non-progressors
525. aggregate(data_meg_pet_cogn3$tau_entorhinal, list(data_meg_pet_cogn3$mci_status), FUN=mean)
526. aggregate(data_meg_pet_cogn3$tau_entorhinal, list(data_meg_pet_cogn3$mci_status), FUN=sd)
528. ### Run Wilcox test
529. wilcox_age=wilcox.test(data_meg_pet_cogn3$age~data_meg_pet_cogn3$mci_status)
530. wilcox_education=wilcox.test(data_meg_pet_cogn3$education~data_meg_pet_cogn3$mci_status)
531. wilcox_tau_meta_roi_meg_delta=wilcox.test(data_meg_pet_cogn3$tau_meta_roi_meg_delta~data_meg_pet_cogn3$mci_status)
532. wilcox_tau_meta_roi_meg_theta=wilcox.test(data_meg_pet_cogn3$tau_meta_roi_meg_theta~data_meg_pet_cogn3$mci_status)
533. wilcox_tau_meta_roi_meg_alpha=wilcox.test(data_meg_pet_cogn3$tau_meta_roi_meg_alpha~data_meg_pet_cogn3$mci_status)
534. wilcox_tau_meta_roi_meg_beta=wilcox.test(data_meg_pet_cogn3$tau_meta_roi_meg_beta~data_meg_pet_cogn3$mci_status)
535. wilcox_tau_meta_roi_meg_gamma1=wilcox.test(data_meg_pet_cogn3$tau_meta_roi_meg_gamma1~data_meg_pet_cogn3$mci_status)
536. wilcox_hipp=wilcox.test(data_meg_pet_cogn3$hippocampal_volume~data_meg_pet_cogn3$mci_status)
537. wilcox_ab4240=wilcox.test(data_meg_pet_cogn3$Ab42_40_ratio~data_meg_pet_cogn3$mci_status)
538. wilcox_ptau217=wilcox.test(data_meg_pet_cogn3$ptau_217~data_meg_pet_cogn3$mci_status)
539. wilcox_ab=wilcox.test(data_meg_pet_cogn3$amyloid_index~data_meg_pet_cogn3$mci_status)
540. wilcox_tau=wilcox.test(data_meg_pet_cogn3$tau_entorhinal~data_meg_pet_cogn3$mci_status)
542. ### Check proportion of subjects for categorical variables (sex and APOE) and run chi squared
543. tt=table(data_meg_pet_cogn3$mci_status,data_meg_pet_cogn3$sex); tt; chi_sex=chisq.test(tt)
544. tt=table(data_meg_pet_cogn3$mci_status,data_meg_pet_cogn3$apoe); tt; chi_apoe=chisq.test(tt)
547. wil_chi_stats=c(wilcox_age$statistic, chi_sex$statistic, wilcox_education$statistic, chi_apoe$statistic,
548. wilcox_tau_meta_roi_meg_delta$statistic, wilcox_tau_meta_roi_meg_theta$statistic, wilcox_tau_meta_roi_meg_alpha$statistic,
549. wilcox_tau_meta_roi_meg_beta$statistic, wilcox_tau_meta_roi_meg_gamma1$statistic,
550. wilcox_hipp$statistic, wilcox_ab4240$statistic, wilcox_ptau217$statistic, wilcox_ab$statistic, wilcox_tau$statistic)
552. p_vals_wilcox=c(wilcox_age$p.value, chi_sex$p.value, wilcox_education$p.value, chi_apoe$p.value,
553. wilcox_tau_meta_roi_meg_delta$p.value, wilcox_tau_meta_roi_meg_theta$p.value, wilcox_tau_meta_roi_meg_alpha$p.value,
554. wilcox_tau_meta_roi_meg_beta$p.value, wilcox_tau_meta_roi_meg_gamma1$p.value,
555. wilcox_hipp$p.value, wilcox_ab4240$p.value, wilcox_ptau217$p.value, wilcox_ab$p.value, wilcox_tau$p.value)
557. p_vals_wilcox_fdr=p.adjust(p_vals_wilcox, method = "fdr")
558. sprintf("%.6f", p_vals_wilcox_fdr)
560. #### XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX ###
566. #### XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX ###
568. ### 6. MCI Cox regression survival analysis (no time interactions)----
571. ### 6.1 Prepare data----
573. ### Compute days to cutoff date (July 31, 2024) for non-progressors to run Cox regression models
574. data_meg_pet_cogn3$diff_days_meg_mci_cox=data_meg_pet_cogn3$diff_days_meg_mci
575. data_meg_pet_cogn3=fill_missing_diff_days(data_meg_pet_cogn3, "MEG_date", "diff_days_meg_mci_cox", cutoff = "2024-07-31", new_col = "diff_years_meg_mci_cox")
578. ### Define MCI list as 0s and 1s to run cox models
579. data_meg_pet_cogn3$mci_status_cox=factor(data_meg_pet_cogn3$mci_status, levels = c(1, 2), labels = c(0, 1)) # Define the factor with levels 0 and 1
580. data_meg_pet_cogn3$mci_status_cox=as.double(data_meg_pet_cogn3$mci_status_cox)
581. data_meg_pet_cogn3$mci_status_cox=data_meg_pet_cogn3$mci_status_cox-1
584. ### Compute days to cutoff date (July 31, 2024) for non-progressors to run Cox regression models
585. data_meg_pet_cogn3$diff_days_plasma_mci_cox=data_meg_pet_cogn3$diff_days_plasma_mci
586. data_meg_pet_cogn3=fill_missing_diff_days(data_meg_pet_cogn3, "Plasma_date", "diff_days_plasma_mci_cox", cutoff = "2024-07-31", new_col = "diff_years_plasma_mci_cox")
590. #### XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX ###
596. #### XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX ###
598. ### 6.2 Run Cox regression models----
600. ### Generate function that computes Cox regression models, get summary,
601. ### compute c-index, check PH assumption, fit survival curves and plot them
602. ### Only works for models that do not include a time interaction
605. run_cox_model <- function(formula, data, model_name) {
606. ### Fit Cox model
607. formula <- as.formula(formula)
608. model <- coxph(formula, data = data); model_summary <- summary(model)
609. ### Check proportional hazard assumption
610. if (model_name != "null_model") {
611. pha=cox.zph(model)}
612. ### Compute C-index
613. c_index <- model_summary$concordance
614. ### Fit survival curves
615. surv_fit <- survfit(model)
616. ### Log-rank test & p-value
617. lr_test <- model_summary$logtest[1]; p_val <- model_summary$logtest[3]
618. ### Plot survival curve
619. plot(surv_fit, xlab = "Time (years)", ylab = "Survival Probability", ylim = c(0, 1))
620. ### Save model summary to a text file
621. summary_filename <- paste0("Cox_model_", model_name, "_summary.txt")
622. #capture.output(model_summary, file = summary_filename)
623. ### Return results for null model (no PH assumption)
624. if (model_name == "null_model") {
625. return(list(model = model, summary = model_summary, c_index = c_index,
626. log_rank_test = lr_test, p_value = p_val))}
628. if (model_name != "null_model") {
629. return(list(model = model, summary = model_summary, ph_assumption=pha, c_index = c_index,
630. log_rank_test = lr_test, p_value = p_val))}
631. }
635. ### Run Cox regressions for each model (main models)
637. ### Clinical
638. results_clin = run_cox_model("Surv(diff_years_meg_mci_cox, mci_status_cox) ~ age_z+sex+education_z+apoe", data_meg_pet_cogn3, "clin")
640. ### MRI
641. results_clin_mri = run_cox_model("Surv(diff_years_meg_mci_cox, mci_status_cox) ~ age_z+sex+education_z+apoe+hippocampal_volume_z", data=data_meg_pet_cogn3, "clin_mri")
643. ### MEG
644. results_clin_meg_delta = run_cox_model("Surv(diff_years_meg_mci_cox, mci_status_cox) ~ age_z+sex+education_z+apoe+tau_meta_roi_meg_delta_z", data_meg_pet_cogn3, "clin_meg_delta")
645. results_clin_meg_theta = run_cox_model("Surv(diff_years_meg_mci_cox, mci_status_cox) ~ age_z+sex+education_z+apoe+tau_meta_roi_meg_theta_z", data_meg_pet_cogn3, "clin_meg_theta")
646. results_clin_meg_alpha = run_cox_model("Surv(diff_years_meg_mci_cox, mci_status_cox) ~ age_z+sex+education_z+apoe+tau_meta_roi_meg_alpha_z", data_meg_pet_cogn3, "clin_meg_alpha")
647. results_clin_meg_beta = run_cox_model("Surv(diff_years_meg_mci_cox, mci_status_cox) ~ age_z+sex+education_z+apoe+tau_meta_roi_meg_beta_z", data_meg_pet_cogn3, "clin_meg_beta")
648. results_clin_meg_gamma1 = run_cox_model("Surv(diff_years_meg_mci_cox, mci_status_cox) ~ age_z+sex+education_z+apoe+tau_meta_roi_meg_gamma1_z", data_meg_pet_cogn3, "clin_meg_gamma1")
649. results_clin_meg_gamma2 = run_cox_model("Surv(diff_years_meg_mci_cox, mci_status_cox) ~ age_z+sex+education_z+apoe+tau_meta_roi_meg_gamma2_z", data_meg_pet_cogn3, "clin_meg_gamma2")
650. results_clin_meg_alpha_gamma1 = run_cox_model("Surv(diff_years_meg_mci_cox, mci_status_cox) ~ age_z+sex+education_z+apoe+tau_meta_roi_meg_alpha_z+tau_meta_roi_meg_gamma1_z", data_meg_pet_cogn3, "clin_meg_alpha_gamma1")
652. ### Plasma
653. results_clin_plasma = run_cox_model("Surv(diff_years_meg_mci_cox, mci_status_cox) ~ age_z+sex+education_z+apoe+Ab42_40_ratio_z+ptau_217_z", data=data_meg_pet_cogn3, "clin_plasma")
655. ### PET
656. results_clin_ab = run_cox_model("Surv(diff_years_meg_mci_cox, mci_status_cox) ~ age_z+sex+education_z+apoe+amyloid_index_z", data=data_meg_pet_cogn3, "clin_ab")
657. results_clin_tau = run_cox_model("Surv(diff_years_meg_mci_cox, mci_status_cox) ~ age_z+sex+education_z+apoe+tau_entorhinal_z", data=data_meg_pet_cogn3, "clin_tau")
659. ### All biomarkers
660. results_all = run_cox_model("Surv(diff_years_meg_mci_cox, mci_status_cox) ~ age_z+sex+education_z+apoe+hippocampal_volume_z+tau_meta_roi_meg_alpha_z+tau_meta_roi_meg_gamma1_z+Ab42_40_ratio_z+ptau_217_z+amyloid_index_z+tau_entorhinal_z", data=data_meg_pet_cogn3, "all")
664. ### Results using full plasma (n=211) Ab (n=226) and Tau PET (n=224) samples
665. results_clin_211 = run_cox_model("Surv(diff_years_plasma_mci_cox, mci_status_cox) ~ age_z+sex+education_z+apoe", data=df_plasma_211, "clin_211")
666. results_clin_plasma_211 = run_cox_model("Surv(diff_years_plasma_mci_cox, mci_status_cox) ~ age_z+sex+education_z+apoe+Ab42_40_ratio_z+ptau_217_z", data=df_plasma_211, "clin_plasma_211")
667. results_clin_226 = run_cox_model("Surv(diff_years_pet_mci_cox, mci_status_cox) ~ age_z+sex+education_z+apoe", data=df_ab_226, "clin_226")
668. results_clin_ab_226 = run_cox_model("Surv(diff_years_pet_mci_cox, mci_status_cox) ~ age_z+sex+education_z+apoe+amyloid_index_z", data=df_ab_226, "clin_ab_226")
669. results_clin_224= run_cox_model("Surv(diff_years_pet_mci_cox, mci_status_cox) ~ age_z+sex+education_z+apoe", data=df_tau_224, "clin_224")
670. results_clin_tau_224= run_cox_model("Surv(diff_years_pet_mci_cox, mci_status_cox) ~ age_z+sex+education_z+apoe+tau_entorhinal_z", data=df_tau_224, "clin_tau_224")
676. ### Get LRT statistic and pvalues from the main models
677. ### Extract LRT statistic values (vs null model)
678. stats_lrt=c(results_clin$log_rank_test, results_clin_meg_alpha_gamma1$log_rank_test, results_clin_mri$log_rank_test, results_clin_plasma$log_rank_test,
679. results_clin_ab$log_rank_test, results_clin_tau$log_rank_test, results_all$log_rank_test)
681. ### Extract p-values and correct for multiple comparisons
682. p_vals_lrt=c(results_clin$p_value, results_clin_meg_alpha_gamma1$p_value, results_clin_mri$p_value, results_clin_plasma$p_value,
683. results_clin_ab$p_value, results_clin_tau$p_value, results_all$p_value)
685. p_vals_fdr_lrt=p.adjust(p_vals_lrt, method = "fdr")
686. sprintf("%.6f", p_vals_fdr_lrt)
692. ### Get LRT statistic and pvalues from the extended models
693. ### Extract LRT statistic values (vs null model)
694. stats_lrt=c(results_clin_plasma_211$log_rank_test, results_clin_ab_226$log_rank_test, results_clin_tau_224$log_rank_test)
696. ### Extract p-values and correct for multiple comparisons
697. p_vals_lrt=c(results_clin_plasma_211$p_value, results_clin_ab_226$p_value, results_clin_tau_224$p_value)
699. p_vals_fdr_lrt=p.adjust(p_vals_lrt, method = "fdr")
700. sprintf("%.6f", p_vals_fdr_lrt)
704. ### Generate function to compare two models using LRT
705. compare_models_lrt <- function(model1, model2) {
706. log_like_1 <- logLik(model1)
707. log_like_2 <- logLik(model2)
708. lr_test <- -2 * (as.numeric(log_like_1) - as.numeric(log_like_2))
709. df <- attr(log_like_2, "df") - attr(log_like_1, "df")
710. p_value <- pchisq(lr_test, df = df, lower.tail = FALSE)
711. return(list(
712. lr_test_statistic = lr_test,
713. df = df,
714. p_value = p_value
715. ))
716. }
721. ### Compare each model against the reference clinical model
722. clin_vs_clin_meg_alpha_gamma1 <- compare_models_lrt(results_clin$model, results_clin_meg_alpha_gamma1$model)
723. clin_vs_clin_mri <- compare_models_lrt(results_clin$model, results_clin_mri$model)
724. clin_vs_clin_plasma <- compare_models_lrt(results_clin$model, results_clin_plasma$model)
725. clin_vs_clin_ab <- compare_models_lrt(results_clin$model, results_clin_ab$model)
726. clin_vs_clin_tau <- compare_models_lrt(results_clin$model, results_clin_tau$model)
727. clin_vs_all <- compare_models_lrt(results_clin$model, results_all$model)
731. #
732. # ### Extract LRT statistic values (vs clinical model)
733. stats_lrt=c(clin_vs_clin_meg_alpha_gamma1$lr_test_statistic, clin_vs_clin_mri$lr_test_statistic, clin_vs_clin_plasma$lr_test_statistic,
734. clin_vs_clin_ab$lr_test_statistic, clin_vs_clin_tau$lr_test_statistic, clin_vs_all$lr_test_statistic)
736. ### Extract p-values and correct for multiple comparisons
737. p_vals_lrt=c(clin_vs_clin_meg_alpha_gamma1$p_value, clin_vs_clin_mri$p_value, clin_vs_clin_plasma$p_value,
738. clin_vs_clin_ab$p_value, clin_vs_clin_tau$p_value, clin_vs_all$p_value)
740. p_vals_fdr_lrt=p.adjust(p_vals_lrt, method = "fdr")
741. sprintf("%.6f", p_vals_fdr_lrt)
748. ### Comparisons of proteinopathy models using the extnded sample
749. clin_vs_clin_plasma_211 <- compare_models_lrt(results_clin_211$model, results_clin_plasma_211$model)
750. clin_vs_clin_ab_226 <- compare_models_lrt(results_clin_226$model, results_clin_ab_226$model)
751. clin_vs_clin_tau_224 <- compare_models_lrt(results_clin_224$model, results_clin_tau_224$model)
755. ### Extract LRT statistic values (vs clinical model)
756. stats_lrt=c(clin_vs_clin_plasma_211$lr_test_statistic, clin_vs_clin_ab_226$lr_test_statistic, clin_vs_clin_tau_224$lr_test_statistic)
758. ### Extract p-values and correct for multiple comparisons
759. p_vals_lrt=c(clin_vs_clin_plasma_211$p_value, clin_vs_clin_ab_226$p_value, clin_vs_clin_tau_224$p_value)
760. p_vals_fdr_lrt=p.adjust(p_vals_lrt, method = "fdr")
761. sprintf("%.6f", p_vals_fdr_lrt)
763. #### XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX ###
769. #### XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX ###
771. ### 6.3 Generate Kaplan-Meier survival curves----
773. ### SVG figures
774. for (i in seq(1,10)) {
776. ### Models
777. model_options = list(results_clin$model, results_clin_mri$model, results_clin_meg_alpha_gamma1$model,
778. results_clin_plasma$model, results_clin_ab$model, results_clin_tau$model,
779. results_all$model, results_clin_plasma_211$model,
780. results_clin_ab_226$model, results_clin_tau_224$model)
782. model_colors = c("gray5", "gray50","#F23F43","#FF7F50", "#1E98DD","#7600BC", "#00BA38",
783. "#FF7F50","#1E98DD","#7600BC")
785. model_names = c("clin","clin_mri","clin_meg","clin_plasma","clin_ab","clin_tau",
786. "all","clin_plasma_211","clin_ab_226","clin_tau_224")
788. ### Select model
789. model = model_options[[i]]
790. model_color = model_colors[i]
791. model_name = model_names[i]
793. ### Axis settings
794. if (i %in% seq(1, 7)) {
795. df_plot = data_meg_pet_cogn3
796. x_limits = c(0, 7.5)
797. y_limits = c(0.6, 1)
798. y_breaks = seq(0.6, 1, by = 0.1)
799. }
801. if (i %in% c(8)) {
802. df_plot = df_plasma_211
803. x_limits = c(0, 11.2)
804. y_limits = c(0.45, 1)
805. y_breaks = seq(0.5, 1, by = 0.1)
806. }
808. if (i %in% c(9)) {
809. df_plot = df_ab_226
810. x_limits = c(0, 7.5)
811. y_limits = c(0.45, 1)
812. y_breaks = seq(0.5, 1, by = 0.1)
813. }
815. if (i %in% c(10)) {
816. df_plot = df_tau_224
817. x_limits = c(0, 7.5)
818. y_limits = c(0.45, 1)
819. y_breaks = seq(0.5, 1, by = 0.1)
820. }
822. ### Generate base plot
823. p = ggsurvplot(
824. fit = survfit(model),
825. data = data_meg_pet_cogn3,
826. xlab = "Time (years)",
827. ylab = "Survival probability",
828. xlim = x_limits,
829. ylim = y_limits,
830. palette = model_color,
831. size = 0.3, # ???? thinner main survival curve
832. conf.int = TRUE,
833. conf.int.fill = model_color,
834. conf.int.alpha = 0.15,
835. conf.int.style = "ribbon",
836. surv.median.line = "none",
837. risk.table = FALSE,
838. legend = "none",
839. censor.size = 1, # ???? smaller + symbols
840. ggtheme = theme_classic(base_family = "Arial") +
841. theme(
842. axis.title = element_text(size = 8),
843. axis.text = element_text(size = 8),
844. axis.line = element_line(linewidth = 0.4),
845. axis.ticks = element_line(linewidth = 0.4),
846. plot.margin = margin(2, 2, 2, 2, "mm")
847. )
848. )
850. ### Refine plot (CI boundary lines slightly thinner)
851. g = p$plot +
852. geom_step(aes(y = upper), color = model_color, linewidth = 0.2, alpha = 1) +
853. geom_step(aes(y = lower), color = model_color, linewidth = 0.2, alpha = 1) +
854. scale_x_continuous(
855. limits = x_limits,
856. breaks = seq(0, ceiling(max(x_limits)), by = 1)
857. ) +
858. scale_y_continuous(
859. limits = y_limits,
860. breaks = y_breaks
861. )
863. ### Save as SVG (Illustrator-ready, editable text)
864. svg_filename <- file.path(path, paste0("cox_surv_curve_", model_name, "_panel.svg"))
866. svglite::svglite(
867. file = svg_filename,
868. width = 4.6 / 2.54, # cm ??? inches
869. height = 4.6 / 2.54
870. )
872. print(g)
873. dev.off()
875. ### Display in R
876. print(g)
877. }
879. #### XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX ###
885. #### XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX ###
887. ### 6.4 Plot survival curves stratified by risk group----
888. ##### Manually change features
891. library(survival)
892. library(survminer)
893. library(svglite)
894. library(grid)
896. ### Compute risk groups
897. risk_scores <- predict(results_all$model, type = "lp")
899. risk_group <- cut(
900. risk_scores,
901. breaks = quantile(risk_scores, probs = c(0, 1/3, 2/3, 1)),
902. labels = c("Low", "Medium", "High"),
903. include.lowest = TRUE
904. )
906. data_meg_pet_cogn3$risk_group <- risk_group
908. ### Fit survival model
909. fit <- survfit(
910. Surv(diff_years_meg_mci_cox, mci_status_cox) ~ risk_group,
911. data = data_meg_pet_cogn3
912. )
914. ### Generate survival plot + risk table
915. p <- ggsurvplot(
916. fit,
917. data = data_meg_pet_cogn3,
918. xlab = "Time (years)",
919. ylab = "Survival probability",
920. legend.labs = c("L", "M", "H"),
921. palette = c("#49beaa", "#ffbe0b", "#f27a7d"),
922. size = 0.4,
923. censor.size = 1,
924. conf.int = TRUE,
925. conf.int.alpha = 0.15,
926. conf.int.style = "ribbon",
927. risk.table = TRUE,
928. risk.table.col = "strata",
929. break.time.by = 1,
930. xlim = c(0, 7.5),
931. tables.theme = theme_classic(base_family = "Arial"),
932. ggtheme = theme_classic(base_family = "Arial"),
933. legend = "none"
934. )
936. ### Style MAIN plot
937. p$plot <- p$plot +
938. scale_x_continuous(limits = c(0, 7.5), breaks = seq(0, 7, by = 1)) +
939. scale_y_continuous(limits = c(0, 1), breaks = seq(0, 1, by = 0.2)) +
940. theme(
941. axis.title = element_text(size = 8),
942. axis.text = element_text(size = 8),
943. axis.line = element_line(linewidth = 0.4),
944. axis.ticks = element_line(linewidth = 0.4),
945. plot.margin = margin(2, 2, 2, 2, "mm")
946. )
948. ### Style RISK TABLE (clean)
949. # Set font sizes for numbers inside table
950. p$table <- p$table +
951. theme_classic(base_family = "Arial") +
952. theme(
953. axis.title = element_blank(),
954. axis.text.x = element_text(size = 8),
955. axis.text.y = element_text(size = 8), # row labels (Low, Medium, High)
956. axis.ticks = element_line(linewidth = 0.4),
957. axis.line = element_line(linewidth = 0.4),
958. plot.margin = margin(1, 2, 2, 2, "mm"),
959. )
961. # Remove "Number at risk" title
962. g <- ggplotGrob(p$table)
963. title_index <- which(g$layout$name == "title")
964. if(length(title_index) > 0) g$grobs[[title_index]] <- nullGrob()
965. p$table <- g # replace table with cleaned grob
967. ### SAVE FILES (SVG ??? Illustrator)
968. # Survival panel
969. svglite::svglite(
970. file = file.path(path, "survival_curve_all.svg"),
971. width = 4.6 / 2.54,
972. height = 4.6 / 2.54
973. )
974. print(p$plot)
975. dev.off()
977. # Risk table panel
978. svglite::svglite(
979. file = file.path(path, "risk_table_all.svg"),
980. width = 5.9 / 2.54,
981. height = 4.6 / 2.54
982. )
983. grid.draw(p$table) # ???? Use grid.draw for grob objects
984. dev.off()
986. ### Optional: print in R
987. print(p$plot)
988. grid.draw(p$table)
990. # ### Log-rank test
991. # survdiff(
992. # Surv(diff_years_meg_mci_cox, mci_status_cox) ~ risk_group,
993. # data = data_meg_pet_cogn3
994. # )
1001. create_time_table <- function(
1002. model,
1003. data,
1004. time_var,
1005. event_var,
1006. time_points = seq(0, 7, by = 1),
1007. output_name = "time_table",
1008. path = "."
1009. ) {
1013. # Ensure proper format
1014. data <- as.data.frame(data)
1016. ### Risk groups
1017. risk_scores <- predict(model, type = "lp")
1019. risk_group <- cut(
1020. risk_scores,
1021. breaks = quantile(risk_scores, probs = c(0, 1/3, 2/3, 1), na.rm = TRUE),
1022. labels = c("Low", "Medium", "High"),
1023. include.lowest = TRUE
1024. )
1026. data$risk_group <- factor(risk_group, levels = c("Low", "Medium", "High"))
1028. ### Survival fit
1029. surv_formula <- as.formula(
1030. paste0("Surv(", time_var, ", ", event_var, ") ~ risk_group")
1031. )
1033. fit <- survfit(surv_formula, data = data)
1035. ### Fixed time table
1036. fit_summary_fixed <- summary(fit, times = time_points)
1038. time_table_fixed <- data.frame(
1039. strata = fit_summary_fixed$strata,
1040. time = fit_summary_fixed$time,
1041. n_risk = fit_summary_fixed$n.risk,
1042. n_event = fit_summary_fixed$n.event,
1043. n_censor = fit_summary_fixed$n.censor
1044. )
1046. ### Cumulative values
1047. time_table_cum <- time_table_fixed %>%
1048. group_by(strata) %>%
1049. arrange(time) %>%
1050. mutate(
1051. cum_events = cumsum(n_event),
1052. cum_censor = cumsum(n_censor)
1053. ) %>%
1054. ungroup() %>%
1055. arrange(strata, time)
1057. ### Save CSV
1058. write.csv(
1059. time_table_cum,
1060. file = file.path(path, paste0(output_name, ".csv")),
1061. row.names = FALSE
1062. )
1064. ### Return table
1065. return(time_table_cum)
1066. }
1071. time_table_clin <- create_time_table(model = results_clin$model, data = data_meg_pet_cogn3,
1072. time_var = "diff_years_meg_mci_cox", event_var = "mci_status_cox", time_points = seq(0, 7, by = 1),
1073. output_name = "time_table_clin", path = path)
1075. time_table_clin_meg <- create_time_table(model = results_clin_meg_alpha_gamma1$model, data = data_meg_pet_cogn3,
1076. time_var = "diff_years_meg_mci_cox", event_var = "mci_status_cox", time_points = seq(0, 7, by = 1),
1077. output_name = "time_table_clin_meg", path = path)
1080. time_table_clin_mri <- create_time_table(model = results_clin_mri$model, data = data_meg_pet_cogn3,
1081. time_var = "diff_years_meg_mci_cox", event_var = "mci_status_cox", time_points = seq(0, 7, by = 1),
1082. output_name = "time_table_clin_mri", path = path)
1085. time_table_clin_plasma <- create_time_table(model = results_clin_plasma$model, data = data_meg_pet_cogn3,
1086. time_var = "diff_years_meg_mci_cox", event_var = "mci_status_cox", time_points = seq(0, 7, by = 1),
1087. output_name = "time_table_clin_plasma", path = path)
1090. time_table_clin_ab <- create_time_table(model = results_clin_ab$model, data = data_meg_pet_cogn3,
1091. time_var = "diff_years_meg_mci_cox", event_var = "mci_status_cox", time_points = seq(0, 7, by = 1),
1092. output_name = "time_table_clin_ab", path = path)
1095. time_table_clin_tau <- create_time_table(model = results_clin_tau$model, data = data_meg_pet_cogn3,
1096. time_var = "diff_years_meg_mci_cox", event_var = "mci_status_cox", time_points = seq(0, 7, by = 1),
1097. output_name = "time_table_clin_tau", path = path)
1100. time_table_all <- create_time_table(model = results_all$model, data = data_meg_pet_cogn3,
1101. time_var = "diff_years_meg_mci_cox", event_var = "mci_status_cox", time_points = seq(0, 7, by = 1),
1102. output_name = "time_table_all", path = path)
1105. time_table_clin_plasma_211 <- create_time_table(model = results_clin_plasma_211$model, data = df_plasma_211,
1106. time_var = "diff_years_plasma_mci_cox", event_var = "mci_status_cox", time_points = seq(0, 11, by = 1),
1107. output_name = "time_table_clin_plasma_211", path = path)
1110. time_table_clin_ab_226 <- create_time_table(model = results_clin_ab_226$model, data = df_ab_226,
1111. time_var = "diff_years_pet_mci_cox", event_var = "mci_status_cox", time_points = seq(0, 7, by = 1),
1112. output_name = "time_table_clin_ab_226", path = path)
1115. time_table_clin_tau_224 <- create_time_table(model = results_clin_tau_224$model, data = df_tau_224,
1116. time_var = "diff_years_pet_mci_cox", event_var = "mci_status_cox", time_points = seq(0, 7, by = 1),
1117. output_name = "time_table_clin_tau_224", path = path)
1120. #### XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX ###
1126. #### XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX ###
1128. ### Assess correlation between risk scores and biomarkers----
1130. ### Change labels for APOE
1131. data_meg_pet_cogn3$apoe <- factor(data_meg_pet_cogn3$apoe, levels = c(0, 1), labels = c("Non-carrier", "Carrier"))
1133. model_options=list(results_clin$model, results_clin_meg_alpha_gamma1$model, results_clin_mri$model, results_clin_plasma$model,
1134. results_clin_ab$model, results_clin_tau$model, results_all$model,
1135. results_clin_plasma_211$model, results_clin_ab_226$model, results_clin_tau_224$model)
1138. compare_risk_simple <- function(variable_name, variable_type, model_idx,
1139. model_options, data, y_label, custom_colors, path,integer_y=FALSE) {
1141. # Compute risk scores
1142. risk_scores <- predict(model_options[[model_idx]], type="lp")
1143. data$risk_scores <- risk_scores
1145. # Continuous variables
1146. if(variable_type == "continuous") {
1147. formula <- as.formula(paste("risk_scores ~", variable_name))
1148. lm_res <- summary(lm(formula, data=data))
1149. t_val <- lm_res$coefficients[variable_name, "t value"]
1150. p_val <- lm_res$coefficients[variable_name, "Pr(>|t|)"]
1152. # Scatterplot with regression line
1153. p_scatter <- ggplot(data, aes(x=.data[[variable_name]], y=risk_scores)) +
1154. geom_point(color=custom_colors[1], alpha=0.5, size=1.5, stroke=0.1) +
1155. geom_smooth(method="lm", se=TRUE, color=custom_colors[1], fill=custom_colors[1], alpha=0.15, linewidth=0.3) +
1156. labs(x=y_label, y="Relative risk score") +
1157. theme_classic(base_family="Arial") +
1158. theme(
1159. axis.title = element_text(size = 8),
1160. axis.text = element_text(size = 8),
1161. axis.line = element_line(linewidth = 0.4),
1162. axis.ticks = element_line(linewidth = 0.4),
1163. plot.margin = margin(2, 2, 2, 2, "mm")
1164. )
1166. if(integer_y) {
1167. p_scatter <- p_scatter +
1168. scale_y_continuous(
1169. breaks = seq(
1170. floor(min(data[[variable_name]], na.rm = TRUE)),
1171. ceiling(max(data[[variable_name]], na.rm = TRUE)),
1172. by = 2
1173. )
1174. )
1175. }
1178. svglite::svglite(file=file.path(path, paste0("scatter_risk_score_", variable_name, ".svg")),
1179. width=4.6/2.54, height=4.6/2.54)
1180. print(p_scatter); dev.off()
1182. return(list(lm_summary=lm_res, t_val=t_val, p_val=p_val))
1183. }
1185. # Categorical variables (sex, apoe)
1186. if(variable_type == "categorical") {
1187. formula <- as.formula(paste("risk_scores ~", variable_name))
1188. t_res <- t.test(formula, data=data)
1190. # Box + violin plot
1191. p_box <- ggplot(data, aes(x=.data[[variable_name]], y=risk_scores, fill=.data[[variable_name]])) +
1192. geom_violin(trim=FALSE, linewidth=0.3) +
1193. geom_boxplot(width=0.3, outlier.shape=NA, linewidth=0.3) +
1194. geom_jitter(width=0.15, color="black", alpha=0.25, size=1.5, stroke = 0.1) +
1195. scale_fill_manual(values=custom_colors) +
1196. labs(x=y_label, y="Relative risk score") +
1197. theme_classic(base_family="Arial") +
1198. theme(
1199. axis.title = element_text(size = 8),
1200. axis.text = element_text(size = 8),
1201. axis.line = element_line(linewidth = 0.4),
1202. axis.ticks = element_line(linewidth = 0.4),
1203. plot.margin = margin(2, 2, 2, 2, "mm"),
1204. legend.position="none")
1206. svglite::svglite(file=file.path(path, paste0("box_risk_score_", variable_name, ".svg")),
1207. width=4.6/2.54, height=4.6/2.54)
1208. print(p_box); dev.off()
1210. return(list(t_test=t_res))
1211. }
1212. }
1220. ### Define parameters for running the models
1222. variables = c("age_z", "sex","education_z", "apoe", "tau_meta_roi_meg_alpha_z","tau_meta_roi_meg_gamma1_z", "hippocampal_volume_z",
1223. "Ab42_40_ratio_z","ptau_217_z","amyloid_index_z","tau_entorhinal_z")
1226. y_labels=c("Age (z-score)", "Sex", "Education (z-score)", "APOE ??4 carrier status",
1227. "Alpha power (z-score)", "Gamma power (z-score)", "Hipp. volume (z-score)",
1228. "A??42/40 ratio (z-score)", "p-tau217 (z-score)", "Neocortical A?? (z-score)", "Entorhinal tau (z-score)")
1230. model=c("clin","clin_meg","clin_mri","clin_plasma","clin_ab","clin_tau","all","clin_plasma_211","clin_ab_226","clin_tau_224")
1232. custom_colors_con=c("#49beaa","#ffbe0b","#f27a7d")
1233. custom_colors_cat=c("#c8f0e7","#ffecb3","#fcdfe0","#49beaa","#ffbe0b","#f27a7d")
1235. model_colors=c("gray5", "#F23F43", "gray50", "#FF7F50", "#1E98DD","#7600BC", "#00BA38", "#FF7F50", "#1E98DD","#7600BC")
1238. My_Theme = theme(plot.title = element_text(hjust = 0.5, size = 16, family = "TT Arial"), axis.text.x = element_text(family = "TT Arial"),
1239. axis.text.y = element_text(family = "TT Arial"), axis.title.x = element_text(family = "TT Arial"),
1240. axis.title.y = element_text(family = "TT Arial"), legend.text=element_text(family = "TT Arial"),
1241. legend.title =element_text(family = "TT Arial"))
1251. # For sex (categorical)
1252. res_sex <- compare_risk_simple(variable_name = "sex", variable_type = "categorical", model_idx = 7,
1253. model_options = model_options, data = data_meg_pet_cogn3, y_label = "Sex",
1254. custom_colors = c("gray80","gray80"), path = path)
1256. # For sex (categorical)
1257. res_apoe <- compare_risk_simple(variable_name = "apoe", variable_type = "categorical", model_idx = 7,
1258. model_options = model_options, data = data_meg_pet_cogn3, y_label = "APOE e4 status",
1259. custom_colors = c("gray80","gray80"), path = path)
1264. # For age_z (continuous)
1265. res_age <- compare_risk_simple(variable_name = "age_z", variable_type = "continuous", model_idx = 7,
1266. model_options = model_options, data = data_meg_pet_cogn3$apoe, y_label = "Age (z-score)",
1267. custom_colors = c("gray5"), path = path)
1269. # For age_z (continuous)
1270. res_age <- compare_risk_simple(variable_name = "education_z", variable_type = "continuous", model_idx = 7,
1271. model_options = model_options, data = data_meg_pet_cogn3, y_label = "Education (z-score)",
1272. custom_colors = c("gray5"), path = path)
1274. # For age_z (continuous)
1275. res_age <- compare_risk_simple(variable_name = "tau_meta_roi_meg_alpha_z", variable_type = "continuous", model_idx = 7,
1276. model_options = model_options, data = data_meg_pet_cogn3, y_label = "Alpha power (z-score)",
1277. custom_colors = c("#F23F43"), path = path)
1279. # For age_z (continuous)
1280. res_age <- compare_risk_simple(variable_name = "tau_meta_roi_meg_gamma1_z", variable_type = "continuous", model_idx = 7,
1281. model_options = model_options, data = data_meg_pet_cogn3, y_label = "Gamma power (z-score)",
1282. custom_colors = c("#F23F43"), path = path)
1284. # For age_z (continuous)
1285. res_age <- compare_risk_simple(variable_name = "hippocampal_volume_z", variable_type = "continuous", model_idx = 7,
1286. model_options = model_options, data = data_meg_pet_cogn3, y_label = "Hipp. volume (z-score)",
1287. custom_colors = c("gray50"), path = path)
1289. # For age_z (continuous)
1290. res_age <- compare_risk_simple(variable_name = "Ab42_40_ratio_z", variable_type = "continuous", model_idx = 7,
1291. model_options = model_options, data = data_meg_pet_cogn3, y_label = "A??42/40 ratio (z-score)",
1292. custom_colors = c("#FF7F50"), path = path)
1294. # For age_z (continuous)
1295. res_age <- compare_risk_simple(variable_name = "ptau_217_z", variable_type = "continuous", model_idx = 7,
1296. model_options = model_options, data = data_meg_pet_cogn3, y_label = "p-tau217 (z-score)",
1297. custom_colors = c("#FF7F50"), path = path, integer_y = TRUE)
1299. # For age_z (continuous)
1300. res_age <- compare_risk_simple(variable_name = "amyloid_index_z", variable_type = "continuous", model_idx = 7,
1301. model_options = model_options, data = data_meg_pet_cogn3, y_label = "Neocortical A?? (z-score)",
1302. custom_colors = c("#1E98DD"), path = path)
1304. # For age_z (continuous)
1305. res_age <- compare_risk_simple(variable_name = "tau_entorhinal_z", variable_type = "continuous", model_idx = 7,
1306. model_options = model_options, data = data_meg_pet_cogn3, y_label = "Entorhinal tau (z-score)",
1307. custom_colors = c("#7600BC"), path = path)
1312. ### Extract t/t statistic values (risk scores)
1313. stats_t_all=c(results_risk_group_age_all$t_val_risk_score,results_risk_group_sex_all$t_value_risk_score,results_risk_group_edu_all$t_val_risk_score,
1314. results_risk_group_apoe_all$t_value_risk_score, results_risk_group_alpha_all$t_val_risk_score, results_risk_group_gamma_all$t_val_risk_score,
1315. results_risk_group_hipp_all$t_val_risk_score, results_risk_group_ab4240_all$t_val_risk_score,results_risk_group_ptau217_all$t_val_risk_score,
1316. results_risk_group_ab_all$t_val_risk_score,results_risk_group_tau_all$t_val_risk_score)
1318. ### Correct p-values for multiple comparisons
1319. p_vals_t_all=c(results_risk_group_age_all$p_value_risk_score,results_risk_group_sex_all$p_value_risk_score,results_risk_group_edu_all$p_value_risk_score,
1320. results_risk_group_apoe_all$p_value_risk_score,results_risk_group_alpha_all$p_value_risk_score, results_risk_group_gamma_all$p_value_risk_score,
1321. results_risk_group_hipp_all$p_value_risk_score, results_risk_group_ab4240_all$p_value_risk_score,results_risk_group_ptau217_all$p_value_risk_score,
1322. results_risk_group_ab_all$p_value_risk_score,results_risk_group_tau_all$p_value_risk_score)
1324. p_vals_fdr_t_all=p.adjust(p_vals_t_all, method = "fdr")
1325. sprintf("%.6f", p_vals_fdr_t_all)
1333. ### 8. MCI Cox regression survival analysis with time interactions----
1335. ### 8.1 Run Cox regressions including time interactions----
1337. ### Generate function to run a Cox model including an interaction term
1339. run_cox_model_time_int <- function(formula, data, model_name) {
1340. formula <- as.formula(formula)
1341. model <- coxph(formula, data = data, tt = function(x, time, ...) x * time)
1342. model_summary <- summary(model)
1343. c_index <- model_summary$concordance
1344. lr_test <- model_summary$logtest[1]
1345. p_val <- model_summary$logtest[3]
1346. summary_filename <- paste0("Cox_model_", model_name, "_summary.txt")
1347. #capture.output(model_summary, file = summary_filename)
1348. return(list(model = model, summary = model_summary, c_index = c_index, log_rank_test = lr_test, p_value = p_val))
1349. }
1353. ### Run Cox regressions including time interactions for each biomarker
1355. ### Age
1356. results_clin_age_nl=run_cox_model_time_int(formula = "Surv(diff_years_meg_mci_cox, mci_status_cox) ~ age_z + sex +
1357. education_z + apoe + tt(age_z)", data = data_meg_pet_cogn3, model_name = "age_tt")
1359. ### Education
1360. results_clin_edu_nl=run_cox_model_time_int(formula = "Surv(diff_years_meg_mci_cox, mci_status_cox) ~ age_z + sex +
1361. education_z + apoe + tt(education_z)", data = data_meg_pet_cogn3, model_name = "edu_tt")
1363. ### MEG delta
1364. results_clin_meg_delta_nl=run_cox_model_time_int(formula = "Surv(diff_years_meg_mci_cox, mci_status_cox) ~ age_z + sex +
1365. education_z + apoe + tau_meta_roi_meg_delta_z + tt(tau_meta_roi_meg_delta_z)", data = data_meg_pet_cogn3, model_name = "delta_tt")
1367. ### MEG theta
1368. results_clin_meg_theta_nl=run_cox_model_time_int(formula = "Surv(diff_years_meg_mci_cox, mci_status_cox) ~ age_z + sex +
1369. education_z + apoe + tau_meta_roi_meg_theta_z + tt(tau_meta_roi_meg_theta_z)", data = data_meg_pet_cogn3, model_name = "theta_tt")
1371. ### MEG alpha
1372. results_clin_meg_alpha_nl=run_cox_model_time_int(formula = "Surv(diff_years_meg_mci_cox, mci_status_cox) ~ age_z + sex +
1373. education_z + apoe + tau_meta_roi_meg_alpha_z + tt(tau_meta_roi_meg_alpha_z)", data = data_meg_pet_cogn3, model_name = "alpha_tt")
1375. ### MEG beta
1376. results_clin_meg_beta_nl=run_cox_model_time_int(formula = "Surv(diff_years_meg_mci_cox, mci_status_cox) ~ age_z + sex +
1377. education_z + apoe + tau_meta_roi_meg_beta_z + tt(tau_meta_roi_meg_beta_z)", data = data_meg_pet_cogn3, model_name = "beta_tt")
1379. ### MEG gamma1
1380. results_clin_meg_gamma1_nl=run_cox_model_time_int(formula = "Surv(diff_years_meg_mci_cox, mci_status_cox) ~ age_z + sex +
1381. education_z + apoe + tau_meta_roi_meg_gamma1_z + tt(tau_meta_roi_meg_gamma1_z)", data = data_meg_pet_cogn3, model_name = "gamma1_tt")
1383. ### MEG gamma2
1384. results_clin_meg_gamma2_nl=run_cox_model_time_int(formula = "Surv(diff_years_meg_mci_cox, mci_status_cox) ~ age_z + sex +
1385. education_z + apoe + tau_meta_roi_meg_gamma2_z + tt(tau_meta_roi_meg_gamma2_z)", data = data_meg_pet_cogn3, model_name = "gamma2_tt")
1387. ### MEG alpha and gamma1
1388. results_clin_meg_alpha_gamma1_nl=run_cox_model_time_int(formula = "Surv(diff_years_meg_mci_cox, mci_status_cox) ~ age_z + sex +
1389. education_z + apoe + tau_meta_roi_meg_alpha_z + tt(tau_meta_roi_meg_alpha_z) + tau_meta_roi_meg_gamma1_z + tt(tau_meta_roi_meg_gamma1_z)",
1390. data = data_meg_pet_cogn3, model_name = "alpha_gamma1_tt")
1392. ### MRI hipp vol
1393. results_clin_mri_nl=run_cox_model_time_int(formula = "Surv(diff_years_meg_mci_cox, mci_status_cox) ~ age_z + sex +
1394. education_z + apoe + hippocampal_volume_z + tt(hippocampal_volume_z)", data = data_meg_pet_cogn3, model_name = "hipp_tt")
1397. ### Plasma Ab42/40 ratio
1398. results_clin_Ab4240_nl=run_cox_model_time_int(formula = "Surv(diff_years_meg_mci_cox, mci_status_cox) ~ age_z + sex +
1399. education_z + apoe + Ab42_40_ratio_z + tt(Ab42_40_ratio_z)", data = data_meg_pet_cogn3, model_name = "Ab42_40_tt")
1401. ### Plasma p-tau217
1402. results_clin_ptau_217_nl=run_cox_model_time_int(formula = "Surv(diff_years_meg_mci_cox, mci_status_cox) ~ age_z + sex +
1403. education_z + apoe + ptau_217_z + tt(ptau_217_z)", data = data_meg_pet_cogn3, model_name = "ptau_217_tt")
1406. ### Ab PET
1407. results_clin_ab_nl=run_cox_model_time_int(formula = "Surv(diff_years_meg_mci_cox, mci_status_cox) ~ age_z + sex +
1408. education_z + apoe + amyloid_index_z + tt(amyloid_index_z)", data = data_meg_pet_cogn3, model_name = "amyloid_index_tt")
1410. ### Tau PET
1411. results_clin_tau_nl=run_cox_model_time_int(formula = "Surv(diff_years_meg_mci_cox, mci_status_cox) ~ age_z + sex +
1412. education_z + apoe + tau_entorhinal_z + tt(tau_entorhinal_z)", data = data_meg_pet_cogn3, model_name = "tau_entorhinal_tt")
1420. ### Run Cox regression models with time interaction for the extended samples
1422. ### Plasma Ab42/40 ratio (n=211)
1423. results_clin_Ab4240_211_nl=run_cox_model_time_int(formula = "Surv(diff_years_plasma_mci_cox, mci_status_cox) ~ age_z + sex +
1424. education_z + apoe + Ab42_40_ratio_z + tt(Ab42_40_ratio_z)", data = df_plasma_211, model_name = "ab4240_211_tt")
1426. ### Plasma p-tau217 (n=211)
1427. results_clin_ptau_217_211_nl=run_cox_model_time_int(formula = "Surv(diff_years_plasma_mci_cox, mci_status_cox) ~ age_z + sex +
1428. education_z + apoe + ptau_217_z + tt(ptau_217_z)", data = df_plasma_211, model_name = "ptau_217_211_tt")
1430. ### Ab PET (n=226)
1431. results_clin_ab_226_nl=run_cox_model_time_int(formula = "Surv(diff_years_pet_mci_cox, mci_status_cox) ~ age_z + sex +
1432. education_z + apoe + amyloid_index_z + tt(amyloid_index_z)", data = df_ab_226, model_name = "ab_226_tt")
1434. ### Tau PET (n=224)
1435. results_clin_tau_224_nl=run_cox_model_time_int(formula = "Surv(diff_years_pet_mci_cox, mci_status_cox) ~ age_z + sex +
1436. education_z + apoe + tau_entorhinal_z + tt(tau_entorhinal_z)", data = df_tau_224, model_name = "tau_224_tt")
1442. ### Combine MEG nl and proteinopathy
1444. ### MEG nl and plasma
1445. results_clin_plasma_meg_alpha_gamma1_nl=run_cox_model_time_int(formula = "Surv(diff_years_meg_mci_cox, mci_status_cox) ~ age_z + sex +
1446. education_z + apoe + Ab42_40_ratio_z + ptau_217_z + tau_meta_roi_meg_alpha_z + tt(tau_meta_roi_meg_alpha_z) + tau_meta_roi_meg_gamma1_z + tt(tau_meta_roi_meg_gamma1_z)",
1447. data = data_meg_pet_cogn3, model_name = "alpha_gamma1_tt")
1449. ### MEG nl and Ab PET
1450. results_clin_ab_meg_alpha_gamma1_nl=run_cox_model_time_int(formula = "Surv(diff_years_meg_mci_cox, mci_status_cox) ~ age_z + sex +
1451. education_z + apoe + tau_meta_roi_meg_alpha_z + tt(tau_meta_roi_meg_alpha_z) + amyloid_index_z + tt(amyloid_index_z) + tau_meta_roi_meg_gamma1_z + tt(tau_meta_roi_meg_gamma1_z)",
1452. data = data_meg_pet_cogn3, model_name = "alpha_gamma1_tt")
1454. ### MEG nl and tau PET
1455. results_clin_tau_meg_alpha_gamma1_nl=run_cox_model_time_int(formula = "Surv(diff_years_meg_mci_cox, mci_status_cox) ~ age_z + sex +
1456. education_z + apoe + tau_meta_roi_meg_alpha_z + tt(tau_meta_roi_meg_alpha_z) + tau_entorhinal_z + tau_meta_roi_meg_gamma1_z + tt(tau_meta_roi_meg_gamma1_z)",
1457. data = data_meg_pet_cogn3, model_name = "alpha_gamma1_tt")
1461. ### Comparisons of proteinopathy models including the gamma1 time interaction
1462. clin_plasma_vs_clin_plasma_meg_alpha_gamma1 <- compare_models_lrt(results_clin_plasma$model, results_clin_plasma_meg_alpha_gamma1_nl$model)
1463. clin_ab_vs_clin_ab_meg_alpha_gamma1 <- compare_models_lrt(results_clin_ab_nl$model, results_clin_ab_meg_alpha_gamma1_nl$model)
1464. clin_tau_vs_clin_tau_meg_alpha_gamma1 <- compare_models_lrt(results_clin_tau$model, results_clin_tau_meg_alpha_gamma1_nl$model)
1468. clin_plasma_vs_clin_plasma_meg_alpha_gamma1=anova(results_clin_plasma$model, results_clin_plasma_meg_alpha_gamma1_nl$model, test = "LRT")
1471. ### Extract p-values and correct for multiple comparisons
1472. stats_lrt=c(clin_plasma_vs_clin_plasma_meg_alpha_gamma1$lr_test_statistic,clin_ab_vs_clin_ab_meg_alpha_gamma1$lr_test_statistic,clin_tau_vs_clin_tau_meg_alpha_gamma1$lr_test_statistic)
1473. p_vals_lrt=c(clin_plasma_vs_clin_plasma_meg_alpha_gamma1$p_value,clin_ab_vs_clin_ab_meg_alpha_gamma1$p_value,clin_tau_vs_clin_tau_meg_alpha_gamma1$p_value)
1474. p_vals_fdr_lrt=p.adjust(p_vals_lrt, method = "fdr")
1475. sprintf("%.6f", p_vals_fdr_lrt)
1479. ### Get LRT statistic and pvalues from the models including a time interaction
1481. ### Extract LRT statistic values (vs null model)
1482. stats_lrt=c(results_clin_age_nl$log_rank_test, results_clin_education_nl$log_rank_test, results_clin_meg_nl$log_rank_test, results_clin_mri_nl$log_rank_test,
1483. results_clin_Ab4240_nl$log_rank_test, results_clin_ptau_217_nl$log_rank_test, results_clin_ab_nl$log_rank_test, results_clin_tau_nl$log_rank_test)
1485. ### Extract p-values and correct for multiple comparisons
1486. p_vals_lrt=c(results_clin_age_nl$p_value, results_clin_education$p_value, results_clin_meg_nl$p_value, results_clin_mri_nl$p_value,
1487. results_clin_Ab4240_nl$p_value, results_clin_ptau_217_nl$p_value, results_clin_ab_nl$p_value, results_clin_tau_nl$p_value)
1489. p_vals_fdr_lrt=p.adjust(p_vals_lrt, method = "fdr")
1492. ### XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX ###
1498. ### XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX ###
1500. ### 8.2 Plot interaction effect: Hazard ratio change over time----
1502. ### Generate function to plot change in Hazard ratio over time from the interaction coefficients
1505. library(ggplot2)
1506. library(svglite)
1508. plot_time_varying_hr_with_ci <- function(cox_model, var_name, time_max = 7, n_points = 100,
1509. color = "#ffbe0b", path = ".", ylim_fixed = "no") {
1511. # Extract coefficients
1512. coef_main <- coef(cox_model)[[var_name]]
1513. coef_tt <- coef(cox_model)[[paste0("tt(", var_name, ")")]]
1515. # Variance-covariance
1516. vcov_matrix <- vcov(cox_model)
1517. var_main <- vcov_matrix[var_name, var_name]
1518. var_tt <- vcov_matrix[paste0("tt(", var_name, ")"), paste0("tt(", var_name, ")")]
1519. cov_main_tt <- vcov_matrix[var_name, paste0("tt(", var_name, ")")]
1521. # Time sequence
1522. time_seq <- seq(0, time_max, length.out = n_points)
1524. # Linear predictor and variance
1525. lp_time <- coef_main + coef_tt * time_seq
1526. var_lp_time <- var_main + (time_seq^2) * var_tt + 2 * time_seq * cov_main_tt
1528. # Confidence intervals
1529. se_lp_time <- sqrt(var_lp_time)
1530. lower_lp <- lp_time - 1.96 * se_lp_time
1531. upper_lp <- lp_time + 1.96 * se_lp_time
1533. # Dataframe
1534. df_hr <- data.frame(
1535. Time = time_seq,
1536. HR = exp(lp_time),
1537. HR_lower = exp(lower_lp),
1538. HR_upper = exp(upper_lp)
1539. )
1541. # HR=1 crossing
1542. time_hr1 <- if (coef_tt != 0) -coef_main / coef_tt else NA
1543. y_max <- if (ylim_fixed == "yes") 2 else max(df_hr$HR_upper, na.rm = TRUE)
1544. y_seg_end <- y_max / 3
1546. # Create ggplot
1547. p <- ggplot(df_hr, aes(x = Time, y = HR)) +
1548. geom_line(color = color, linewidth = 0.6) +
1549. geom_ribbon(aes(ymin = HR_lower, ymax = HR_upper), fill = color, alpha = 0.2) +
1550. geom_hline(yintercept = 1, linetype = "dashed", color = "gray50", linewidth = 0.4) +
1551. labs(x = "Years from biomarker", y = "Hazard ratio") +
1552. theme_classic(base_family = "Arial") + # Use system font
1553. theme(
1554. axis.title = element_text(size = 8),
1555. axis.text = element_text(size = 8),
1556. axis.line = element_line(linewidth = 0.4),
1557. axis.ticks = element_line(linewidth = 0.4),
1558. plot.margin = margin(2, 2, 2, 2, "mm")
1559. )
1561. # Y-axis breaks
1562. if (ylim_fixed == "yes") {
1563. p <- p + scale_y_continuous(limits = c(0, 2), breaks = 0:2)
1564. } else {
1565. p <- p + scale_y_continuous(breaks = seq(0, ceiling(max(df_hr$HR_upper)), by = 1))
1566. }
1568. # HR=1 annotation
1569. if (!is.na(time_hr1) && time_hr1 >= 0 && time_hr1 <= time_max) {
1570. p <- p +
1571. annotate("segment",
1572. x = time_hr1, xend = time_hr1,
1573. y = 1, yend = 1 - (y_seg_end * 0.4),
1574. linetype = "dashed", color = "gray50", linewidth = 0.4) +
1575. annotate("text",
1576. x = time_hr1,
1577. y = 1 - (y_seg_end * 0.4),
1578. label = "HR=1",
1579. size = 2.2, # ??? 8 pt
1580. vjust = 1.2)
1581. }
1583. # Save as SVG (editable text in Illustrator)
1584. svg_filename <- paste0(path, "/", var_name, "_panel.svg")
1585. svglite::svglite(file = svg_filename, width = 4.6/2.54, height = 4.6/2.54) # cm ??? inches
1586. print(p)
1587. dev.off()
1589. return(p)
1590. }
1595. ## Generate plots for each biomarker
1596. plot_time_varying_hr_with_ci(results_clin_age_nl$model, "age_z", 7, 100, "#ffbe0b", path, "no")
1597. plot_time_varying_hr_with_ci(results_clin_edu_nl$model, "education_z", 7, 100, "#ffbe0b", path, "no")
1599. plot_time_varying_hr_with_ci(results_clin_meg_delta_nl$model, "tau_meta_roi_meg_delta_z", 7, 100, "#ffbe0b", path, "no")
1600. plot_time_varying_hr_with_ci(results_clin_meg_theta_nl$model, "tau_meta_roi_meg_theta_z", 7, 100, "#ffbe0b", path, "no")
1601. plot_time_varying_hr_with_ci(results_clin_meg_alpha_nl$model, "tau_meta_roi_meg_alpha_z", 7, 100, "#ffbe0b", path, "no")
1602. plot_time_varying_hr_with_ci(results_clin_meg_beta_nl$model, "tau_meta_roi_meg_beta_z", 7, 100, "#ffbe0b", path, "no")
1603. plot_time_varying_hr_with_ci(results_clin_meg_gamma1_nl$model, "tau_meta_roi_meg_gamma1_z", 7, 100, "#ffbe0b", path, "no")
1604. plot_time_varying_hr_with_ci(results_clin_meg_gamma2_nl$model, "tau_meta_roi_meg_gamma2_z", 7, 100, "#ffbe0b", path, "no")
1606. plot_time_varying_hr_with_ci(results_clin_mri_nl$model, "hippocampal_volume_z", 7, 100, "#ffbe0b", path, "no")
1607. plot_time_varying_hr_with_ci(results_clin_Ab4240_nl$model, "Ab42_40_ratio_z", 7, 100, "#ffbe0b", path,"yes")
1608. plot_time_varying_hr_with_ci(results_clin_ptau_217_nl$model, "ptau_217_z", 7, 100, "#ffbe0b", path, "no")
1609. plot_time_varying_hr_with_ci(results_clin_ab_nl$model, "amyloid_index_z", 7, 100, "#ffbe0b", path, "no")
1610. plot_time_varying_hr_with_ci(results_clin_tau_nl$model, "tau_entorhinal_z", 7, 100, "#ffbe0b", path, "no")
1613. ### Generate plots for the extended samples for plasma and PET
1614. plot_time_varying_hr_with_ci(results_clin_Ab4240_211_nl$model, "Ab42_40_ratio_z", 11, 100, "#ffbe0b", path,"no")
1615. plot_time_varying_hr_with_ci(results_clin_ptau_217_211_nl$model, "ptau_217_z", 11, 100, "#ffbe0b", path,"no")
1616. plot_time_varying_hr_with_ci(results_clin_ab_226_nl$model, "amyloid_index_z", 7, 100, "#ffbe0b", path, "no")
1617. plot_time_varying_hr_with_ci(results_clin_tau_224_nl$model, "tau_entorhinal_z", 7, 100, "#ffbe0b", path, "no")
1619. ### XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX ###
1625. ### XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX ###
1627. ### 8.3 Compute relative risk and plot linear predictor as a function of time----
1629. ### Function: Compute relative risks + Illustrator-friendly plotting
1630. compute_relative_risks <- function(data, biomarker, time_var, event_var, covariates = NULL,
1631. time_points = 1:7, id_var = "pscid", plot = TRUE, model_name = NULL,
1632. path = ".", include_ci = TRUE, overlay_baseline = FALSE,
1633. save_csv = FALSE, csv_filename = NULL, scale = "log") {
1635. # Check scale
1636. if (!scale %in% c("log", "exp")) stop("scale must be 'log' or 'exp'")
1638. # Build Cox formula
1639. all_covariates <- paste(c(covariates, biomarker, paste0("tt(", biomarker, ")")), collapse = " + ")
1640. cox_formula <- as.formula(paste0("Surv(", time_var, ", ", event_var, ") ~ ", all_covariates))
1642. # Fit Cox model with time interaction
1643. cox_model <- coxph(cox_formula, data = data, tt = function(x, time, ...) x * time)
1645. # Compute linear predictors
1646. linear_predictor <- predict(cox_model, type = "lp", newdata = data)
1648. coefs <- coef(cox_model)
1649. vcov_mat <- vcov(cox_model)
1651. biom_coef <- coefs[biomarker]
1652. interaction_name <- paste0("tt(", biomarker, ")")
1653. interaction_coef <- coefs[interaction_name]
1655. se_biom <- sqrt(vcov_mat[biomarker, biomarker])
1656. se_inter <- sqrt(vcov_mat[interaction_name, interaction_name])
1657. cov_biom_inter <- vcov_mat[biomarker, interaction_name]
1659. lp_no_main <- linear_predictor - data[[biomarker]] * biom_coef
1661. # Compute relative risks for each time point
1662. relative_risks <- lapply(time_points, function(t) {
1663. biom <- data[[biomarker]]
1664. lp <- lp_no_main + biom * (biom_coef + interaction_coef * t)
1666. if (include_ci) {
1667. se_lp <- sqrt((biom^2) * (se_biom^2 + t^2 * se_inter^2 + 2 * t * cov_biom_inter))
1668. lower <- lp - 1.96 * se_lp
1669. upper <- lp + 1.96 * se_lp
1671. if (scale == "exp") {
1672. return(data.frame(lp = exp(lp), lower = exp(lower), upper = exp(upper)))
1673. } else {
1674. return(data.frame(lp = lp, lower = lower, upper = upper))
1675. }
1676. } else {
1677. return(data.frame(lp = if (scale == "exp") exp(lp) else lp))
1678. }
1679. })
1681. # Combine results
1682. rr_df <- do.call(rbind, lapply(seq_along(time_points), function(i) {
1683. df <- relative_risks[[i]]
1684. df[[id_var]] <- data[[id_var]]
1685. df$biomarker_value <- data[[biomarker]]
1686. df$time <- time_points[i]
1687. return(df)
1688. }))
1690. rr_df <- rr_df[, c(id_var, "time", "biomarker_value", "lp", if (include_ci) c("lower", "upper") else NULL)]
1691. names(rr_df)[names(rr_df) == "lp"] <- "relative_risk"
1693. # Save CSV if requested
1694. if (save_csv && !is.null(csv_filename)) {
1695. write.csv(rr_df, file = file.path(path, csv_filename), row.names = FALSE)
1696. }
1698. # Plotting (Illustrator-friendly)
1699. if (plot && !is.null(model_name)) {
1701. # Colors
1702. parula_colors <- c("#352A87", "#3B52A1", "#3F7FBA", "#469BBA", "#58B89E", "#84CA79", "#D9D93A")
1703. if (length(time_points) != length(parula_colors)) stop("Number of time points must match parula_colors length")
1704. names(parula_colors) <- as.character(time_points)
1706. # Labels
1707. xlab <- switch(biomarker,
1708. "age_z" = "Age (z-score)",
1709. "education_z" = "Education (z-score)",
1710. "tau_meta_roi_meg_delta_z" = "Delta power (z-score)",
1711. "tau_meta_roi_meg_theta_z" = "Theta power (z-score)",
1712. "tau_meta_roi_meg_alpha_z" = "Alpha power (z-score)",
1713. "tau_meta_roi_meg_beta_z" = "Beta power (z-score)",
1714. "tau_meta_roi_meg_gamma1_z" = "Gamma power (z-score)",
1715. "tau_meta_roi_meg_gamma2_z" = "Gamma2 power (z-score)",
1716. "hippocampal_volume_z" = "Hipp. volume (z-score)",
1717. "Ab42_40_ratio_z" = "A??42/40 ratio (z-score)",
1718. "ptau_217_z" = "p-tau217 (z-score)",
1719. "amyloid_index_z" = "Neocortical A?? (z-score)",
1720. "tau_entorhinal_z" = "Entorhinal tau (z-score)",
1721. biomarker)
1723. ylab <- if (scale == "exp") "Relative Risk (HR)" else "Relative risk score"
1725. # Base plot
1726. p <- ggplot(rr_df, aes(x = biomarker_value, y = relative_risk, color = factor(time))) +
1727. geom_point(size = 0.6, alpha = 0.7) +
1728. geom_smooth(method = "lm", aes(group = time), se = FALSE, linewidth = 0.6) +
1729. scale_color_manual(values = parula_colors, name = "Time") +
1730. guides(color = guide_legend(nrow = 1))
1732. # CI ribbons
1733. if (include_ci) {
1734. p <- p + geom_ribbon(aes(ymin = lower, ymax = upper, fill = factor(time)), alpha = 0.15, color = NA) +
1735. scale_fill_manual(values = parula_colors, guide = "none")
1736. }
1738. # Adaptive integer y-axis
1739. p <- p + scale_y_continuous(
1740. breaks = function(x) {
1741. bks <- pretty(x)
1742. bks_int <- unique(round(bks))
1743. bks_int[bks_int >= min(x) & bks_int <= max(x)]
1744. },
1745. labels = scales::label_number(accuracy = 1)
1746. )
1748. # Color scale
1749. p <- p + scale_color_manual(values = parula_colors, name = "Time")
1751. # Theme
1752. p <- p +
1753. labs(x = xlab, y = ylab) +
1754. theme_classic(base_family = "Arial") +
1755. theme(
1756. axis.title = element_text(size = 8),
1757. axis.text = element_text(size = 8),
1758. axis.line = element_line(linewidth = 0.4),
1759. axis.ticks = element_line(linewidth = 0.4),
1760. legend.position = "none",
1761. legend.title = element_text(size = 8),
1762. legend.text = element_text(size = 7),
1763. legend.key.width = unit(0.4, "cm"),
1764. legend.spacing.x = unit(0.05, "cm"),
1765. legend.margin = margin(0, 0, 0, 0),
1766. legend.box.margin = margin(0, 0, 0, 0),
1767. plot.margin = margin(2, 2, 2, 2, "mm")
1768. )
1770. # Save as SVG for Illustrator
1771. svg_filename <- file.path(path, paste0(model_name, "_panel.svg"))
1772. svglite::svglite(file = svg_filename, width = 4.6/2.54, height = 4.6/2.54) # cm ??? inches
1773. print(p)
1774. dev.off()
1776. # Also show plot in R
1777. print(p)
1778. }
1780. return(rr_df)
1781. }
1785. rr_df <- compute_relative_risks(data = data_meg_pet_cogn3, biomarker = "age_z",
1786. time_var = "diff_years_meg_mci_cox", event_var = "mci_status_cox", covariates = c("sex", "education_z", "apoe"),
1787. time_points = 1:7, id_var = "pscid", model_name = "linear_predictor_age", path = path,
1788. include_ci = FALSE, scale = "log", save_csv = FALSE, csv_filename = "relative_risks_age.csv")
1790. rr_df <- compute_relative_risks(data = data_meg_pet_cogn3, biomarker = "education_z",
1791. time_var = "diff_years_meg_mci_cox", event_var = "mci_status_cox", covariates = c("age", "sex", "apoe"),
1792. time_points = 1:7, id_var = "pscid", model_name = "linear_predictor_education", path = path,
1793. include_ci = FALSE, scale = "log", save_csv = FALSE, csv_filename = "relative_risks_education.csv")
1795. rr_df <- compute_relative_risks(data = data_meg_pet_cogn3, biomarker = "tau_meta_roi_meg_delta_z",
1796. time_var = "diff_years_meg_mci_cox", event_var = "mci_status_cox", covariates = c("age", "sex", "education_z", "apoe"),
1797. time_points = 1:7, id_var = "pscid", model_name = "linear_predictor_meg_delta", path = path,
1798. include_ci = FALSE, scale = "log", save_csv = FALSE, csv_filename = "relative_risks_meg_delta.csv")
1800. rr_df <- compute_relative_risks(data = data_meg_pet_cogn3, biomarker = "tau_meta_roi_meg_theta_z",
1801. time_var = "diff_years_meg_mci_cox", event_var = "mci_status_cox", covariates = c("age", "sex", "education_z", "apoe"),
1802. time_points = 1:7, id_var = "pscid", model_name = "linear_predictor_meg_theta", path = path,
1803. include_ci = FALSE, scale = "log", save_csv = FALSE, csv_filename = "relative_risks_meg_theta.csv")
1805. rr_df <- compute_relative_risks(data = data_meg_pet_cogn3, biomarker = "tau_meta_roi_meg_alpha_z",
1806. time_var = "diff_years_meg_mci_cox", event_var = "mci_status_cox", covariates = c("age", "sex", "education_z", "apoe"),
1807. time_points = 1:7, id_var = "pscid", model_name = "linear_predictor_meg_alpha", path = path,
1808. include_ci = FALSE, scale = "log", save_csv = FALSE, csv_filename = "relative_risks_meg_alpha.csv")
1811. rr_df <- compute_relative_risks(data = data_meg_pet_cogn3, biomarker = "tau_meta_roi_meg_beta_z",
1812. time_var = "diff_years_meg_mci_cox", event_var = "mci_status_cox", covariates = c("age", "sex", "education_z", "apoe"),
1813. time_points = 1:7, id_var = "pscid", model_name = "linear_predictor_meg_beta", path = path,
1814. include_ci = FALSE, scale = "log", save_csv = FALSE, csv_filename = "relative_risks_meg_beta.csv")
1817. rr_df <- compute_relative_risks(data = data_meg_pet_cogn3, biomarker = "tau_meta_roi_meg_gamma1_z",
1818. time_var = "diff_years_meg_mci_cox", event_var = "mci_status_cox", covariates = c("age", "sex", "education_z", "apoe"),
1819. time_points = 1:7, id_var = "pscid", model_name = "linear_predictor_meg_gamma1", path = path,
1820. include_ci = FALSE, scale = "log", save_csv = FALSE, csv_filename = "relative_risks_meg_gamma1.csv")
1823. rr_df <- compute_relative_risks(data = data_meg_pet_cogn3, biomarker = "tau_meta_roi_meg_gamma2_z",
1824. time_var = "diff_years_meg_mci_cox", event_var = "mci_status_cox", covariates = c("age", "sex", "education_z", "apoe"),
1825. time_points = 1:7, id_var = "pscid", model_name = "linear_predictor_meg_gamma2", path = path,
1826. include_ci = FALSE, scale = "log", save_csv = FALSE, csv_filename = "relative_risks_meg_gamma2.csv")
1829. rr_df <- compute_relative_risks(data = data_meg_pet_cogn3, biomarker = "hippocampal_volume_z",
1830. time_var = "diff_years_meg_mci_cox", event_var = "mci_status_cox", covariates = c("age", "sex", "education_z", "apoe"),
1831. time_points = 1:7, id_var = "pscid", model_name = "linear_predictor_hipp_vol", path = path,
1832. include_ci = FALSE, scale = "log", save_csv = FALSE, csv_filename = "relative_risks_hipp_vol.csv")
1835. rr_df <- compute_relative_risks(data = data_meg_pet_cogn3, biomarker = "Ab42_40_ratio_z",
1836. time_var = "diff_years_meg_mci_cox", event_var = "mci_status_cox", covariates = c("age", "sex", "education_z", "apoe"),
1837. time_points = 1:7, id_var = "pscid", model_name = "linear_predictor_Ab4240", path = path,
1838. include_ci = FALSE, scale = "log", save_csv = FALSE, csv_filename = "relative_risks_Ab4240.csv")
1841. rr_df <- compute_relative_risks(data = data_meg_pet_cogn3, biomarker = "ptau_217_z",
1842. time_var = "diff_years_meg_mci_cox", event_var = "mci_status_cox", covariates = c("age", "sex", "education_z", "apoe"),
1843. time_points = 1:7, id_var = "pscid", model_name = "linear_predictor_ptau217", path = path,
1844. include_ci = FALSE, scale = "log", save_csv = FALSE, csv_filename = "relative_risks_ptau217.csv")
1847. rr_df <- compute_relative_risks(data = data_meg_pet_cogn3, biomarker = "amyloid_index_z",
1848. time_var = "diff_years_meg_mci_cox", event_var = "mci_status_cox", covariates = c("age", "sex", "education_z", "apoe"),
1849. time_points = 1:7, id_var = "pscid", model_name = "linear_predictor_amyloid_idx", path = path,
1850. include_ci = FALSE, scale = "log", save_csv = FALSE, csv_filename = "relative_risks_amyloid_idx.csv")
1853. rr_df <- compute_relative_risks(data = data_meg_pet_cogn3, biomarker = "tau_entorhinal_z",
1854. time_var = "diff_years_meg_mci_cox", event_var = "mci_status_cox", covariates = c("age", "sex", "education_z", "apoe"),
1855. time_points = 1:7, id_var = "pscid", model_name = "linear_predictor_tau_ent", path = path,
1856. include_ci = FALSE, scale = "log", save_csv = FALSE, csv_filename = "relative_risks_tau_ent.csv")
1860. rr_df <- compute_relative_risks(data = df_plasma_211, biomarker = "Ab42_40_ratio_z",
1861. time_var = "diff_years_plasma_mci_cox", event_var = "mci_status_cox", covariates = c("age", "sex", "education_z", "apoe"),
1862. time_points = 1:7, id_var = "pscid", model_name = "linear_predictor_Ab4240_211", path = path,
1863. include_ci = FALSE, scale = "log", save_csv = FALSE, csv_filename = "relative_risks_Ab4240_211.csv")
1866. rr_df <- compute_relative_risks(data = df_plasma_211, biomarker = "ptau_217_z",
1867. time_var = "diff_years_plasma_mci_cox", event_var = "mci_status_cox", covariates = c("age", "sex", "education_z", "apoe"),
1868. time_points = 1:7, id_var = "pscid", model_name = "linear_predictor_ptau217_211", path = path,
1869. include_ci = FALSE, scale = "log", save_csv = FALSE, csv_filename = "relative_risks_ptau217_211.csv")
1872. rr_df <- compute_relative_risks(data = df_ab_226, biomarker = "amyloid_index_z",
1873. time_var = "diff_years_pet_mci_cox", event_var = "mci_status_cox", covariates = c("age", "sex", "education_z", "apoe"),
1874. time_points = 1:7, id_var = "pscid", model_name = "linear_predictor_amyloid_idx_226", path = path,
1875. include_ci = FALSE, scale = "log", save_csv = FALSE, csv_filename = "relative_risks_amyloid_idx_226.csv")
1878. rr_df <- compute_relative_risks(data = df_tau_224, biomarker = "tau_entorhinal_z",
1879. time_var = "diff_years_pet_mci_cox", event_var = "mci_status_cox", covariates = c("age", "sex", "education_z", "apoe"),
1880. time_points = 1:7, id_var = "pscid", model_name = "linear_predictor_tau_ent_224", path = path,
1881. include_ci = FALSE, scale = "log", save_csv = FALSE, csv_filename = "relative_risks_tau_ent_224.csv")
1883. ### XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX ###
1889. ### XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX ###
1891. ### 9. Run stepwise Cox regression----
1893. ### including alpha and gamma1 and amyloid time interactions and hippocampal volume
1894. step_cox=step(coxph(Surv(diff_years_meg_mci_cox, mci_status_cox) ~ age_z + sex + education_z + apoe +
1895. hippocampal_volume_z + tau_meta_roi_meg_alpha_z + tt(tau_meta_roi_meg_alpha_z) + tau_meta_roi_meg_gamma1_z + tt(tau_meta_roi_meg_gamma1_z) +
1896. Ab42_40_ratio_z + ptau_217_z + amyloid_index_z + tt(amyloid_index_z) + tau_entorhinal_z, data=data_meg_pet_cogn3))
1900. ### Selected model
1901. final_cox <- coxph(Surv(diff_years_meg_mci_cox, mci_status_cox) ~ age_z + sex + hippocampal_volume_z + tau_meta_roi_meg_gamma1_z +
1902. tt(tau_meta_roi_meg_alpha_z) + amyloid_index_z + tt(amyloid_index_z) + tau_entorhinal_z, data = data_meg_pet_cogn3)
1904. ### XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX ###
1910. ### XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX ###
1912. ### Supplementary materials----
1914. ### 10. AIC-based model comparison----
1916. ### Set color palette for the diff matrix
1917. color_scale <- colorRampPalette(c("#71bfff","white", "#fc9c9c"))(100)
1919. ### Get AIC values
1920. AIC_1=AIC(results_clin$model); AIC_2=AIC(results_clin_meg_alpha_gamma1$model);
1921. AIC_3=AIC(results_clin_meg_alpha_nl$model); AIC_4=AIC(results_clin_meg_gamma1_nl$model);
1922. AIC_5=AIC(results_clin_mri$model); AIC_6=AIC(results_clin_plasma$model);
1923. AIC_7=AIC(results_clin_ab$model); AIC_8=AIC(results_clin_ab_nl$model);
1924. AIC_9=AIC(results_clin_tau$model); AIC_10=AIC(results_all$model);
1927. values_AIC <- c(AIC_1, AIC_2, AIC_3, AIC_4, AIC_5, AIC_6, AIC_7,AIC_8, AIC_9, AIC_10)
1931. ##### Calculate differences in AIC/BIC values between models
1932. diff_matrix <- outer(values_AIC, values_AIC, "-"); diag(diff_matrix) <- NA
1933. models=c("Ref","Ref MEG","Ref MEG A (t)","Ref MEG G (t)","Ref MRI","Ref plasma","Ref A??","Ref A?? (t)","Ref tau","All");
1934. colnames(diff_matrix)=models; rownames(diff_matrix)=models; print(diff_matrix)
1935. diff_matrix=diff_matrix*(-1)
1940. library(corrplot)
1941. library(svglite)
1943. plot_aic_matrix <- function(
1944. model_list,
1945. model_names,
1946. output_name = "AIC_matrix",
1947. path = "."
1948. ) {
1951. # Compute AIC values
1952. values_AIC <- sapply(model_list, AIC)
1954. # Difference matrix
1955. diff_matrix <- outer(values_AIC, values_AIC, "-")
1956. diff_matrix <- diff_matrix * (-1)
1958. colnames(diff_matrix) <- model_names
1959. rownames(diff_matrix) <- model_names
1961. diff_vals <- diff_matrix
1962. diag(diff_matrix) <- 0
1964. # Create star matrix (LOWER TRIANGLE ONLY)
1965. star_matrix <- matrix("", nrow = nrow(diff_vals), ncol = ncol(diff_vals))
1967. for(i in 1:nrow(diff_vals)) {
1968. for(j in 1:ncol(diff_vals)) {
1970. if(i > j) { # ???? only lower triangle
1971. val <- diff_vals[i, j]
1973. if(!is.na(val)) {
1974. if(abs(val) > 10) {
1975. star_matrix[i, j] <- "**"
1976. } else if(abs(val) > 2) {
1977. star_matrix[i, j] <- "*"
1978. }
1979. }
1980. }
1981. }
1982. }
1985. # Color scale
1986. color_scale <- colorRampPalette(c("#71bfff", "white", "#fc9c9c"))(100)
1988. # Save SVG (slightly larger ??? better scaling)
1989. svglite::svglite(
1990. file = file.path(path, paste0(output_name, ".svg")),
1991. width = 7 / 2.54,
1992. height = 7 / 2.54
1993. )
1995. # ???? remove ALL margins (outer + inner)
1996. par(mar = c(0, 0, 0, 0), mai = c(0, 0, 0, 0))
1998. # Plot matrix
1999. corrplot(
2000. corr = diff_matrix,
2001. method = "color",
2002. col = color_scale,
2003. is.corr = FALSE,
2004. type = "lower",
2005. diag = FALSE,
2006. outline = FALSE,
2007. addCoef.col = NULL,
2008. tl.cex = 0.6, # ~8 pt labels
2009. tl.col = "black",
2010. tl.srt = 45,
2011. tl.offset = 0.1, # ???? tighter ??? reduces top whitespace
2012. cl.pos = "n", # no legend
2013. bg = "white",
2014. addgrid.col = "black"
2015. )
2017. # Add stars (lower triangle only)
2018. n <- nrow(diff_matrix)
2020. for(i in 1:n) {
2021. for(j in 1:n) {
2023. if(i > j) { # ???? enforce lower triangle
2024. label <- star_matrix[i, j]
2026. if(label != "") {
2027. text(
2028. x = j,
2029. y = n - i + 1,
2030. labels = label,
2031. cex = 0.6 # ~8 pt
2032. )
2033. }
2034. }
2035. }
2036. }
2038. dev.off()
2040. return(diff_vals)
2041. }
2046. model_list <- list(
2047. results_clin$model,
2048. results_clin_meg_alpha_gamma1$model,
2049. results_clin_meg_alpha_nl$model,
2050. results_clin_meg_gamma1_nl$model,
2051. results_clin_mri$model,
2052. results_clin_plasma$model,
2053. results_clin_ab$model,
2054. results_clin_ab_nl$model,
2055. results_clin_tau$model,
2056. results_all$model
2057. )
2059. model_names <- c(
2060. "Ref","Ref MEG","Ref MEG A (t)","Ref MEG G (t)",
2061. "Ref MRI","Ref plasma","Ref A??","Ref A?? (t)",
2062. "Ref tau","All"
2063. )
2065. diff_matrix <- plot_aic_matrix(
2066. model_list = model_list,
2067. model_names = model_names,
2068. output_name = "AIC_difference_models",
2069. path = path
2070. )
2075. ### XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX ###
2081. ### XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX ###
2083. ### 11. Concordance index----
2085. ### Set plot colors
2086. # Define model labels and colors
2087. models <- c("Ref", "Ref MEG", "Ref MEG A (t)", "Ref MEG G (t)", "Ref MRI", "Ref plasma",
2088. "Ref A??", "Ref A?? (t)", "Ref tau", "All biomarkers")
2091. # Create dataframe with custom model labels
2092. df_c_index <- data.frame(
2093. Model = factor(models, levels = models), # Set desired order and labels
2094. C = c(results_clin$c_index[1], results_clin_meg_alpha_gamma1$c_index[1], results_clin_meg_alpha_nl$c_index[1], results_clin_meg_gamma1_nl$c_index[1],
2095. results_clin_mri$c_index[1], results_clin_plasma$c_index[1], results_clin_ab$c_index[1],
2096. results_clin_ab_nl$c_index[1], results_clin_tau$c_index[1], results_all$c_index[1]),
2097. seC = c(results_clin$c_index[2], results_clin_meg_alpha_gamma1$c_index[2], results_clin_meg_alpha_nl$c_index[2], results_clin_meg_gamma1_nl$c_index[2],
2098. results_clin_mri$c_index[2], results_clin_plasma$c_index[2], results_clin_ab$c_index[2],
2099. results_clin_ab_nl$c_index[2], results_clin_tau$c_index[2], results_all$c_index[2])
2100. )
2104. c_index_colors <- c("gray20", "#F23F43","#F23F43","#F23F43", "gray50", "#FF7F50",
2105. "#1E98DD","#1E98DD", "#7600BC", "#00BA38")
2107. plot_c_index <- function(
2108. df_c_index,
2109. colors,
2110. output_name = "c_index_models",
2111. path = "."
2112. ) {
2114. # Plot
2115. p <- ggplot(df_c_index, aes(x = Model, y = C, fill = Model)) +
2117. geom_bar(
2118. stat = "identity",
2119. color = "black",
2120. linewidth = 0.3 # ???? thin outline
2121. ) +
2123. geom_errorbar(
2124. aes(ymin = C - seC, ymax = C + seC),
2125. width = 0.15,
2126. linewidth = 0.3 # ???? thin error bars
2127. ) +
2129. scale_fill_manual(values = colors) +
2131. labs(
2132. x = NULL,
2133. y = "Concordance index"
2134. ) +
2136. coord_cartesian(ylim = c(0.5, 0.85)) +
2138. theme_classic(base_family = "Arial") +
2140. theme(
2141. legend.position = "none",
2143. axis.text.x = element_text(
2144. angle = 70,
2145. hjust = 1,
2146. size = 6
2147. ),
2148. axis.text.y = element_text(size = 8),
2150. axis.title.y = element_text(size = 8),
2151. axis.title.x = element_blank(),
2153. axis.line = element_line(linewidth = 0.4),
2154. axis.ticks = element_line(linewidth = 0.4),
2156. plot.margin = margin(2, 2, 2, 2, "mm") # ???? tight margins
2157. )
2159. # Save SVG
2160. svglite::svglite(
2161. file = file.path(path, paste0(output_name, ".svg")),
2162. width = 4.6 / 2.54, # ???? match your other panels
2163. height = 4.6 / 2.54
2164. )
2166. print(p)
2167. dev.off()
2169. return(p)
2170. }
2174. models <- c("Ref", "Ref MEG", "Ref MEG A (t)", "Ref MEG G (t)",
2175. "Ref MRI", "Ref plasma", "Ref A??", "Ref A?? (t)",
2176. "Ref tau", "All biomarkers")
2178. df_c_index <- data.frame(
2179. Model = factor(models, levels = models),
2180. C = c(
2181. results_clin$c_index[1],
2182. results_clin_meg_alpha_gamma1$c_index[1],
2183. results_clin_meg_alpha_nl$c_index[1],
2184. results_clin_meg_gamma1_nl$c_index[1],
2185. results_clin_mri$c_index[1],
2186. results_clin_plasma$c_index[1],
2187. results_clin_ab$c_index[1],
2188. results_clin_ab_nl$c_index[1],
2189. results_clin_tau$c_index[1],
2190. results_all$c_index[1]
2191. ),
2192. seC = c(
2193. results_clin$c_index[2],
2194. results_clin_meg_alpha_gamma1$c_index[2],
2195. results_clin_meg_alpha_nl$c_index[2],
2196. results_clin_meg_gamma1_nl$c_index[2],
2197. results_clin_mri$c_index[2],
2198. results_clin_plasma$c_index[2],
2199. results_clin_ab$c_index[2],
2200. results_clin_ab_nl$c_index[2],
2201. results_clin_tau$c_index[2],
2202. results_all$c_index[2]
2203. )
2204. )
2206. plot_c_index(
2207. df_c_index = df_c_index,
2208. colors = c_index_colors,
2209. output_name = "c_index_models",
2210. path = path
2211. )
2215. ### XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX ###
2221. ### XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX ###
2223. ### 12. Plot correlations coefficients between model features----
2225. ### Function to compute p-values for correlations
2226. cor_pvalues <- function(mat) {
2227. p_mat <- matrix(NA, ncol = ncol(mat), nrow = nrow(mat))
2228. colnames(p_mat) <- colnames(mat)
2229. rownames(p_mat) <- rownames(mat)
2231. for (i in 1:(ncol(mat) - 1)) {
2232. for (j in (i + 1):ncol(mat)) {
2233. test <- cor.test(mat[, i], mat[, j], use = "complete.obs")
2234. p_mat[i, j] <- test$p.value
2235. p_mat[j, i] <- test$p.value # Symmetric matrix
2236. }
2237. }
2238. diag(p_mat) <- NA # No p-values on the diagonal
2239. return(p_mat)
2240. }
2242. # Compute correlation matrix
2243. cor_matrix <- cor(data_meg_pet_cogn3[, c("age","education","tau_meta_roi_meg_alpha","tau_meta_roi_meg_gamma1","hippocampal_volume",
2244. "Ab42_40_ratio","ptau_217","amyloid_index",
2245. "tau_entorhinal")], use = "complete.obs")
2247. # Compute p-value matrix
2248. p_matrix <- cor_pvalues(data_meg_pet_cogn3[, c("age","education","tau_meta_roi_meg_alpha","tau_meta_roi_meg_gamma1","hippocampal_volume",
2249. "Ab42_40_ratio","ptau_217","amyloid_index",
2250. "tau_entorhinal")])
2251. p_matrix=p_matrix[c(seq(1,7)),seq(1,7)]
2255. # Labels
2256. labels <- c("Age", "Education", "Alpha power", "Gamma power", "Hipp. vol.",
2257. "A??42/40 ratio", "p-tau217", "Necortical A??", "Entorhinal tau")
2259. rownames(cor_matrix) <- labels
2260. colnames(cor_matrix) <- labels
2261. rownames(p_matrix) <- labels
2262. colnames(p_matrix) <- labels
2264. # Color scale for correlation matrix
2265. #color_scale_corr <- colorRampPalette(c("#3A97A3", "white", "#ff84e1"))(100)
2266. color_scale_corr <- colorRampPalette(c("#71bfff","white", "#fc9c9c"))(100)
2268. # Color scale for p-values (e.g., significant values darker)
2269. color_scale_pval <- colorRampPalette(c("white", "red"))(100)
2274. plot_correlation_matrix <- function(
2275. data,
2276. variables,
2277. labels = NULL,
2278. output_name = "correlation_matrix",
2279. path = "."
2280. ) {
2282. library(corrplot)
2283. library(svglite)
2285. # Subset data
2286. mat <- data[, variables]
2288. # Compute correlation matrix
2289. cor_matrix <- cor(mat, use = "complete.obs")
2291. # Labels
2292. if(!is.null(labels)) {
2293. rownames(cor_matrix) <- labels
2294. colnames(cor_matrix) <- labels
2295. }
2297. # Create star matrix (LOWER TRIANGLE ONLY)
2298. star_matrix <- matrix("", nrow = nrow(cor_matrix), ncol = ncol(cor_matrix))
2300. for(i in 1:nrow(cor_matrix)) {
2301. for(j in 1:ncol(cor_matrix)) {
2303. if(i > j) {
2304. val <- cor_matrix[i, j]
2306. if(!is.na(val)) {
2307. if(abs(val) > 0.5) {
2308. star_matrix[i, j] <- "**"
2309. } else if(abs(val) > 0.3) {
2310. star_matrix[i, j] <- "*"
2311. }
2312. }
2313. }
2314. }
2315. }
2317. # Color scale
2318. color_scale <- colorRampPalette(c("#71bfff","white","#fc9c9c"))(100)
2320. # Save SVG
2321. svglite::svglite(
2322. file = file.path(path, paste0(output_name, ".svg")),
2323. width = 7 / 2.54,
2324. height = 7 / 2.54
2325. )
2327. par(mar = c(0, 0, 0, 0), mai = c(0, 0, 0, 0))
2329. # Plot matrix (no numbers)
2330. corrplot(
2331. cor_matrix,
2332. method = "color",
2333. col = color_scale,
2334. is.corr = TRUE,
2335. type = "lower",
2336. diag = FALSE,
2337. outline = FALSE,
2338. addCoef.col = NULL,
2339. tl.cex = 0.6, # ~8 pt labels
2340. tl.col = "black",
2341. tl.srt = 45,
2342. tl.offset = 0.1,
2343. cl.pos = "b", # ???? color legend at bottom
2344. cl.cex = 0.6, # ~8 pt legend text
2345. cl.length = 5, # cleaner ticks
2346. cl.align.text = "c",
2347. bg = "white",
2348. addgrid.col = "black"
2349. )
2351. # Add stars (lower triangle only)
2352. n <- nrow(cor_matrix)
2354. for(i in 1:n) {
2355. for(j in 1:n) {
2357. if(i > j) {
2358. label <- star_matrix[i, j]
2360. if(label != "") {
2361. text(
2362. x = j,
2363. y = n - i + 1,
2364. labels = label,
2365. cex = 0.6 # ~8 pt
2366. )
2367. }
2368. }
2369. }
2370. }
2372. dev.off()
2374. return(cor_matrix)
2375. }
2381. vars <- c(
2382. "age","education","tau_meta_roi_meg_alpha","tau_meta_roi_meg_gamma1",
2383. "hippocampal_volume","Ab42_40_ratio","ptau_217",
2384. "amyloid_index","tau_entorhinal"
2385. )
2387. labels <- c(
2388. "Age", "Education", "Alpha power", "Gamma power",
2389. "Hipp. vol.", "A??42/40 ratio", "p-tau217",
2390. "Neocortical A??", "Entorhinal tau"
2391. )
2393. plot_correlation_matrix(
2394. data = data_meg_pet_cogn3,
2395. variables = vars,
2396. labels = labels,
2397. output_name = "corr_features_stars",
2398. path = path
2399. )

Cox_regression_analysis.R at commit 179a41c, no license · at the source