Sleep increases firing rate modulation during interictal epileptiform discharges in mesial temporal structures.

1. install.packages("circular")
2. install.packages("ggplot2")
3. install.packages("units")
4. install.packages("reshape2")
5. install.packages("circlize")
6. install.packages("ggthemes")
7. install.packages("lemon")
8. install.packages("egg")
9. install.packages("readxl")
10. install.packages('Rcpp')
11. install.packages("equatiomatic")
12. install.packages("kableExtra")
13. install.packages("CircStats")
14. install.packages("emmeans")
15. install.packages("lmerTest")
16. install.packages("devtools")
17. install.packages("sjPlot")
18. install.packages("ggpubr")
19. install.packages("viridis")
20. install.packages("latex2exp")
22. #if(!require(devtools)) install.packages("devtools")
23. #devtools::install_github("kassambara/ggpubr")
24. #setTimeLimit(100000); setSessionTimeLimit(10000)
25. #devtools::install_github("strengejacke/sjPlot")
27. # ggpval
29. # install.packages("spiralize")
30. # library(spiralize)
31. library(ggplot2)
32. #library("cowplot")
33. #library("gridExtra")
34. library(plyr)
35. library(reshape2)
36. library(RColorBrewer)
37. #library("ggthemes")
38. #library(lemon)
39. # library("sjPlot")
40. library(emmeans)
41. library(gridExtra)
42. library(gtable)
43. library(grid)
44. library(egg)
45. library(ggpubr)
46. require(dplyr)
47. library(xtable)
48. library("readxl")
49. library(Rcpp)
50. library(equatiomatic)
51. # library(CircStats)
52. library(circular)
53. library(kableExtra)
54. library("sjPlot") # plot_model
55. library(viridis)
56. library(latex2exp)
59. #####################
60. ## Support function #
61. #####################
63. # replace subsequent fields with NAN for visualization purposes
64. cleanf <- function(x){
65. oldx <- c(FALSE, x[-1]==x[-length(x)]) # is the value equal to the previous?
66. res <- x
67. res[oldx] <- NA
68. res}
71. ###############################
72. # Latex table: Clinical table #
73. ###############################
75. data_clinical <- read.csv("D:/Dropbox/Apps/Overleaf/Hspike/tables/clinical.csv")
76. data_clinical$ID <- NULL
77. data_clinical$Label <- NULL
78. colnames(data_clinical) <- c("Patient","Sex","Age","Onset","Type","SOZ","MRI","PET","SPECT","Medication","Implantation","Pre-implantation surgery" )
79. kbl(data_clinical, "latex", booktabs = T, label = "clinical",
80. caption = "Clinical summary")%>%
81. kable_styling(latex_options = c("scale_down"))%>%
82. kable_styling(latex_options = c("HOLD_position"))%>%
83. column_spec(1, width = "3em") %>%
84. column_spec(2, width = "1em") %>%
85. column_spec(3, width = "1em") %>%
86. column_spec(4, width = "2em") %>%
87. column_spec(5, width = "5em") %>%
88. column_spec(6, width = "5em") %>%
89. column_spec(7, width = "5em") %>%
90. column_spec(8, width = "5em") %>%
91. column_spec(9, width = "5em") %>%
92. column_spec(10, width = "5em") %>%
93. column_spec(11, width = "5em") %>%
94. collapse_rows(columns = 1) %>%
95. footnote(general_title = "",
96. footnote_as_chunk = TRUE,
97. threeparttable = TRUE,
98. escape = FALSE,
99. general = c("
100. AED: Antiepileptic drugs,
101. CBZ: Carbamazepine,
102. ESL: Eslicarabazepine,
103. FBTCS: Focal to Bilateral Tonic-Clonic Seizures,
104. FIAS: Focal Impaired Awareness Seizure,
105. FSWLA: Focal seizures without loss of awareness,
106. IEDs: Interictal Epileptiform discharges,
107. LCS: Lacosamide,
108. LTG: Lamotrigine,
109. MRI: Indications from Magnetic Resonance Imaging,
110. OXC: Oxicarbazepine,
111. PER: Perampanel,
112. PMG: Polymicrogyria,
113. PNH: Periventricular Nodular Heterotopia,
114. SNH: Subcortical Nodular Heterotopia,
115. SOZ: Seizure Onset Zone,
116. TPM: Topiramate
117. VPA: Valproic Acid,
118. ZNG: Zonisamide."))%>%
119. save_kable("D:/Dropbox/Apps/Overleaf/Hspike/tables/clinical.tex")
122. ###############################################
123. # Latex table: electrode anatomical locations #
124. ###############################################
126. options(knitr.kable.NA = '')
128. macro <- read.csv("D:/Dropbox/Apps/Overleaf/Hspike/tables/macro_anatomical.csv", fileEncoding = 'UTF-8-BOM')
129. macro$ID <- NULL
130. macro$Label <- NULL
131. # clean.cols <- c("Patient")
132. # macro[clean.cols] <- lapply(macro[clean.cols], cleanf)
134. micro <- read.csv("D:/Dropbox/Apps/Overleaf/Hspike/tables/micro_anatomical.csv", fileEncoding = 'UTF-8-BOM')
135. micro$ID <- NULL
136. micro$Label <- NULL
137. # clean.cols <- c("Patient")
138. # micro[clean.cols] <- lapply(micro[clean.cols], cleanf)
140. locations <- bind_rows(macro,micro)
142. kbl(locations, "latex", booktabs = T, linesep = "", label = 'anatomical',
143. caption = "Anatomical locations of macro and micro electrodes")%>%
144. kable_styling(latex_options = c("HOLD_position"))%>%
145. pack_rows("Macro contacts", 1, 14) %>% # latex_gap_space = "2em"
146. pack_rows("Micro electrodes", 15, 27) %>%
147. row_spec(1, extra_latex_after = "\\cline{2-5}") %>%
148. row_spec(3, extra_latex_after = "\\cline{2-5}") %>%
149. row_spec(5, extra_latex_after = "\\cline{2-5}") %>%
150. row_spec(7, extra_latex_after = "\\cline{2-5}") %>%
151. row_spec(9, extra_latex_after = "\\cline{2-5}") %>%
152. row_spec(11, extra_latex_after = "\\cline{2-5}") %>%
153. row_spec(13, extra_latex_after = "\\cline{2-5}") %>%
154. row_spec(15, extra_latex_after = "\\cline{2-5}") %>%
155. row_spec(16, extra_latex_after = "\\cline{2-5}") %>%
156. row_spec(18, extra_latex_after = "\\cline{2-5}") %>%
157. row_spec(20, extra_latex_after = "\\cline{2-5}") %>%
158. row_spec(21, extra_latex_after = "\\cline{2-5}") %>%
159. row_spec(23, extra_latex_after = "\\cline{2-5}") %>%
160. row_spec(24, extra_latex_after = "\\cline{2-5}") %>%
161. row_spec(25, extra_latex_after = "\\cline{2-5}") %>%
163. save_kable("D:/Dropbox/Apps/Overleaf/Hspike/tables/anatomical.tex")
166. #########################
167. # Detection performance #
168. #########################
170. options(knitr.kable.NA = '')
171. data <- read.csv(file="Z:/analyses/stephen.whitmarsh/data/hspike/performance.txt", sep=',', header=TRUE, dec='.', na.strings = " ")
172. data$Patient[9] = "\\textit{Mean}"
173. data$Patient[10] = "\\textit{Std.}"
174. data[,2] = round(data[,2],digits=0)
175. data[,3] = round(data[,3],digits=0)
176. data[,4] = round(data[,4],digits=1)
177. data[,5] = round(data[,5],digits=1)
178. data[,6] = round(data[,6],digits=0)
179. data[,7] = round(data[,7],digits=1)
180. #data[,8] = round(data[,8],digits=1)
181. #data[,9] = round(data[,9],digits=1)
182. #data[,10] = round(data[,10],digits=0)
183. #data[,11] = round(data[,11],digits=1)
186. kbl(data, "latex", booktabs = T, linesep = "", label = 'performance', escape = FALSE,
187. col.names = c("Patient","24 hrs.", "24 hrs.","Hit (\\%)","FA (\\%)", "Total", "Total hrs."),
188. # col.names = c("Patient","24hrs", "24hrs","Hit (\\%)","FA (\\%)", "Total","24hrs","Hit (\\%)","FA (\\%)", "Total", "Total hrs."),
189. caption = "Automatic IED detection performance")%>%
190. kable_styling(latex_options = c("HOLD_position"))%>%
191. # kable_styling(latex_options = c("scale_down"))%>%
192. row_spec(8, hline_after = TRUE)%>%
193. add_header_above(c(" " = 1, "Visual" = 1, "Automatic detection" = 4)) %>%
194. # add_header_above(c(" " = 1, "Visual" = 1, "All templates" = 4, "Selected templates" = 4)) %>%
195. save_kable("D:/Dropbox/Apps/Overleaf/Hspike/tables/performance.tex")
198. ########################################
199. # Latex table: Electrode locations MNI #
200. ########################################
202. options(knitr.kable.NA = '')
203. data <- read.csv(file="//lexport/iss02.charpier/analyses/stephen.whitmarsh/data/hspike/MNI_table.txt", sep=',', header=TRUE, dec='.', na.strings = " ")
204. data <- data[!data$color == 0, ] # only make table of used contacts
205. data <- data[, c('patient','electrode','contact','X','Y','Z')]
206. colnames(data) = c('Patient','Electrode','Contact', 'X', 'Y', 'Z')
207. rownames(data) <- NULL
208. clean.cols <- c("Patient","Electrode")
209. data[clean.cols] <- lapply(data[clean.cols], cleanf)
211. kbl(data, "latex", booktabs = T, linesep = "", label = 'MNI',
212. caption = "Anatomical locations of micro electrodes") %>%
213. kable_styling(font_size = 6) %>%
214. kable_styling(latex_options = c("HOLD_position"))%>%
215. row_spec(5, extra_latex_after = "\\cline{1-6}") %>%
216. row_spec(10, extra_latex_after = "\\cline{2-6}") %>%
217. row_spec(14, extra_latex_after = "\\cline{1-6}") %>%
218. row_spec(19, extra_latex_after = "\\cline{1-6}") %>%
219. row_spec(24, extra_latex_after = "\\cline{2-6}") %>%
220. row_spec(29, extra_latex_after = "\\cline{1-6}") %>%
221. row_spec(34, extra_latex_after = "\\cline{2-6}") %>%
222. row_spec(39, extra_latex_after = "\\cline{1-6}") %>%
223. row_spec(44, extra_latex_after = "\\cline{2-6}") %>%
224. row_spec(49, extra_latex_after = "\\cline{1-6}") %>%
225. row_spec(54, extra_latex_after = "\\cline{1-6}") %>%
226. row_spec(59, extra_latex_after = "\\cline{2-6}") %>%
227. save_kable("D:/Dropbox/Apps/Overleaf/Hspike/tables/MNI.tex")
230. ###################################
231. # Latex table: Time in sleepstage #
232. ###################################
234. options(knitr.kable.NA = '')
236. # normalize by time spend in sleep stages
237. data <- read.csv(file="Z:/analyses/stephen.whitmarsh/data/hspike/hypnogram_duration.txt", sep=',', header=TRUE, dec='.', na.strings = " ")
238. data <- data[, c("patient", "part", "TOTAL", "pPHASE_3", "pPHASE_2", "pPHASE_1", "pREM", "pWASO")]
239. colnames(data) = c("Patient", "Night","Total (hrs.)", "S3", "S2", "S1", "REM", "WASO")
240. data <- round(data, digits=1)
242. df <- c("\\textit{Mean}", NA, round(mean(data$Total), digits = 1), round(mean(data$S3), digits = 1), round(mean(data$S2), digits = 1), round(mean(data$S1), digits = 1), round(mean(data$REM), digits = 1), round(mean(data$WASO), digits = 1))
243. data <- rbind(data, df)
244. # df <- c("\\textit{Std.}", NA, round(mean(data$Total), digits = 1), round(sd(data$S3), digits = 0), round(sd(data$S2), digits = 0), round(sd(data$S1), digits = 0), round(sd(data$REM), digits = 0), round(sd(data$WASO), digits = 0))
245. # data <- rbind(data, df)
247. clean.cols <- c("Patient")
248. data[clean.cols] <- lapply(data[clean.cols], cleanf)
250. library(dplyr)
251. kbl(data, "latex", booktabs = T, linesep = "", label = 'stageduration', escape = FALSE,
252. caption = "Time spend in sleep stages.", digits=2) %>%
253. kable_styling(latex_options = c("HOLD_position"))%>%
254. add_header_above(c(" " = 3, "Sleep stage (%)" = 5)) %>%
255. row_spec(3, extra_latex_after = "\\cline{2-8}") %>%
256. row_spec(6, extra_latex_after = "\\cline{2-8}") %>%
257. row_spec(9, extra_latex_after = "\\cline{2-8}") %>%
258. row_spec(12, extra_latex_after = "\\cline{2-8}") %>%
259. row_spec(15, extra_latex_after = "\\cline{2-8}") %>%
260. row_spec(18, extra_latex_after = "\\cline{2-8}") %>%
261. row_spec(21, extra_latex_after = "\\cline{2-8}") %>%
262. row_spec(24, hline_after = TRUE) %>%
263. save_kable("D:/Dropbox/Apps/Overleaf/Neuronal mechanisms underlying the modulation of interictal epileptic activity during sleep/tables/stageduration.tex")
266. ####################################
267. # Latex table: number of SUA / MUA #
268. ####################################
269. #
270. # options(knitr.kable.NA = '')
271. #
272. # data_MUA <- read.csv(file="//lexport/iss02.charpier/analyses/stephen.whitmarsh/data/hspike/DataMUASUA.txt", sep=',', header=TRUE, dec='.', na.strings = " ")
273. # data_MUA <- data_MUA[, c("PatientNr", "Part", "nrSUA", "nrMUA")]
274. # colnames(data_MUA) = c("Patient", "Night","SUA", "MUA")
275. # df <- c("\\textit{Sum}", NA, sum(data_MUA$SUA), sum(data_MUA$MUA))
276. # data_MUA <- rbind(data_MUA, df)
277. #
278. # clean.cols <- c("Patient")
279. # data_MUA[clean.cols] <- lapply(data_MUA[clean.cols], cleanf)
280. #
281. # kbl(data_MUA, "latex", booktabs = T, linesep = "", label = 'SUAMUA', escape = FALSE,
282. # caption = "Number of putatively isolated single units (SUA) and number of multiunits (MUA)")%>%
283. # kable_styling(latex_options = c("HOLD_position"))%>%
284. # row_spec(3, extra_latex_after = "\\cline{2-4}") %>%
285. # row_spec(6, extra_latex_after = "\\cline{2-4}") %>%
286. # row_spec(9, extra_latex_after = "\\cline{2-4}") %>%
287. # row_spec(12, extra_latex_after = "\\cline{2-4}") %>%
288. # row_spec(15, extra_latex_after = "\\cline{2-4}") %>%
289. # row_spec(18, extra_latex_after = "\\cline{2-4}") %>%
290. # row_spec(21, extra_latex_after = "\\cline{2-4}") %>%
291. # row_spec(24, hline_after = TRUE)%>%
292. # save_kable("D:/Dropbox/Apps/Overleaf/Hspike/tables/SUAMUA.tex")
293. #
295. #### ADD RESPONSIVE UNITS #######
296. library(tidyr)
298. # prepare data
299. data_psth <- read.csv(file="Z:/analyses/stephen.whitmarsh/data/hspike/psth_table.txt", sep=',', header=TRUE, dec='.', na.strings = " ")
300. data_psth$hyplabel <- factor(data_psth$hyplabel, levels = c("Pre", "Post", "REM", "Wake", "S1", "S2", "S3"))
301. data_psth$Type <- factor(data_psth$SUA)
303. data_sel <- setNames(aggregate(data_psth$responsive, by = c(list(data_psth$Patient, data_psth$part, data_psth$unit, data_psth$Type)), mean), c("Patient", "part", "unit", "Type", "responsive"))
304. data_sel$responsive = as.integer(data_sel$responsive > 0)
306. t = data_sel %>% count(Patient, part, responsive, Type)
307. t2 <- pivot_wider(t, names_from = "Type", names_prefix = "SUA", values_from = "n")
308. t3 <- pivot_wider(t2, names_from = "responsive", names_prefix = "responsive", values_from = c("SUA1", "SUA0"))
309. t3[is.na(t3)] <- 0
311. df <- c("\\textit{Sum}", NA, sum(t3$SUA1_responsive1), sum(t3$SUA1_responsive0), sum(t3$SUA0_responsive1), sum(t3$SUA0_responsive0))
313. data_MUA <- rbind(t3, df)
315. colnames(data_MUA) = c("Patient", "Night","Responsive", "Unresponsive", "Responsive", "Unresponsive")
317. clean.cols <- c("Patient")
318. data_MUA[clean.cols] <- lapply(data_MUA[clean.cols], cleanf)
320. kbl(data_MUA, "latex", booktabs = T, linesep = "", label = 'SUAMUA', escape = FALSE,
321. caption = "Number of responsive or unresponsive putatively isolated single units (SUA) and multiunits (MUA)")%>%
322. kable_styling(latex_options = c("HOLD_position"))%>%
323. add_header_above(c(" ", " ", "SUA" = 2, "MUA" = 2)) %>%
324. row_spec(3, extra_latex_after = "\\cline{1-6}") %>%
325. row_spec(6, extra_latex_after = "\\cline{1-6}") %>%
326. row_spec(9, extra_latex_after = "\\cline{1-6}") %>%
327. row_spec(12, extra_latex_after = "\\cline{1-6}") %>%
328. row_spec(15, extra_latex_after = "\\cline{1-6}") %>%
329. row_spec(18, extra_latex_after = "\\cline{1-6}") %>%
330. row_spec(21, extra_latex_after = "\\cline{1-6}") %>%
331. row_spec(24, hline_after = TRUE) %>%
332. save_kable("D:/Dropbox/Apps/Overleaf/Neuronal mechanisms underlying the modulation of interictal epileptic activity during sleep/tables/SUAMUA.tex")
335. #######################
336. # LFP power circadian #
337. #######################
339. # prepare data
340. data_power <- read.csv(file="Z:/analyses/stephen.whitmarsh/data/hspike/power_table_long.txt", sep=',', header=TRUE, dec='.', na.strings = " ")
341. data_power$hyplabel <- factor(data_power$hyplabel, ordered = TRUE, levels = c("NO_SCORE", "REM", "AWAKE", "PHASE_1", "PHASE_2", "PHASE_3"))
342. # data_power$band <- factor(data_power$band, ordered = TRUE, levels = c("delta", "theta", "alpha", "beta", "delta_div_alpha"))
343. data_power$band <- factor(data_power$band, ordered = TRUE, levels = c("Delta1", "Delta2"))
344. data_power$Patient <- factor(data_power$patient, levels = c(8:1))
345. data_power$part <- factor(data_power$part)
348. # for finding the median for plotting
349. # data_power$bin <- as.integer(cut(data_power$minute, seq(0, 60*24, by = 1)))
350. # data_binned_fine <- setNames(aggregate(data_power$power, c(list(data_power$Patient), list(data_power$bin), list(data_power$band)), mean), c("Patient", "bin", "band", "power"))
351. # data_binned_fine <- as.data.frame(data_binned_fine %>% group_by(Patient, band) %>% mutate(Npower = (power-min(power))/max(power-min(power))))
352. # ggplot(data=data_binned_fine[data_binned_fine$band == "Delta1" & data_binned_fine$Patient == 7, ], aes(x = bin, y = power, fill=Patient, col=Patient)) + geom_smooth()
355. # bin for polar representation
356. data_power$bin <- as.integer(cut(data_power$minute, seq(0, 24*60, by = 60)))
358. # duplicate midnight to connect in figure
359. temp <- subset(data_power, bin == 24)
360. temp$bin <- 0
361. data_power <- bind_rows(data_power, temp)
363. data_binned <- setNames(aggregate(data_power$power, c(list(data_power$Patient), list(data_power$bin), list(data_power$band)), mean), c("Patient", "bin", "band", "power"))
364. data_binned <- as.data.frame(data_binned %>% group_by(Patient, band) %>% mutate(Npower = (power-min(power))/max(power-min(power))))
367. # data_binned$rad <- data_binned$bin / 24 * pi
368. # c <- circular(control, units = "degrees", template = "geographics")
369. # data_binned$rad
370. # d <- density.circular(data_binned$bin[data_binned$patient == 1], bw = 50)
372. ###
373. # data_binned <- as.data.frame(data_binned %>% group_by(Patient, band) %>% mutate(Npower = (mean(power)-power)/sd(power)))
375. ####
376. # getCentroid <- function(x, width = 1) {
377. # A <- x * width # area of each bar
378. # xc <- seq(width/2, length(x), 1) # x coordinates of center of bars
379. # yc <- x/2 # y coordinatey
380. #
381. # cx <- sum(xc * A) / sum(A)
382. # cy <- sum(yc * A) / sum(A)
383. # return(list(x = cx, y = cy))
384. # }
385. # points(getCentroid(x), col = 'red', pch = 19)
386. ####
388. # plot
389. data_binned$title1 = "Slow waves (0.1-2.5 Hz)"
390. data_binned$title2 = "Delta (2.5-4 Hz)"
391. # data_binned$title2 = "Theta (5-7Hz)"
392. # data_binned$title3 = "Alpha (8-14Hz)"
393. # data_binned$title4 = "Delta (1-4Hz) / Alpha (8-14Hz)"
395. polarpowerplots <- list()
397. polarpowerplots[[1]] <-
398. ggplot(data=data_binned[data_binned$band == "delta", ], aes(x = bin, y = Npower, fill=Patient, col=Patient)) +
399. scale_fill_brewer(palette = "Set2", direction = -1) + scale_color_brewer(palette = "Set2", direction = -1) +
400. geom_vline(xintercept = seq(0, 24, by = 3), colour = "grey90") +
401. geom_hline(yintercept = seq(1, 8, by = 1), colour = "grey90") +
402. geom_ribbon(aes(ymin = (9-as.numeric(Patient)), ymax=Npower*2+(9-as.numeric(Patient))), alpha=0.8, colour = NA) +
403. theme_article() +
404. theme(
405. panel.border = element_blank(),
406. legend.text = element_blank(),
407. axis.ticks = element_blank(),
408. axis.text.x = element_blank(),
409. axis.text.y = element_blank(),
410. axis.title.x = element_blank(),
411. axis.title.y = element_blank()) +
412. scale_x_continuous(breaks=seq(0, 21, by = 3),
413. labels = c("0" = "00:00", "3" = "", "6" = "06:00", "9" = "", "12" = "12:00", "15" = "", "18" = "18:00", "21" = "")) +
414. coord_polar(theta = "x", start = 0, clip="off") +
415. ylim(-2, 10) +
416. facet_wrap(~title1)
418. polardelta1power <-
419. ggplot(data=data_binned[data_binned$band == "Delta1", ], aes(x = bin, y = Npower, fill=Patient, col=Patient)) +
420. scale_fill_brewer(palette = "Set2", direction = -1) + scale_color_brewer(palette = "Set2", direction = -1) +
421. geom_vline(xintercept = seq(0, 24, by = 3), colour = "grey90") +
422. geom_hline(yintercept = seq(1, 8, by = 1), colour = "grey90") +
423. geom_ribbon(aes(ymin = (9-as.numeric(Patient)), ymax=Npower*1.5+(9-as.numeric(Patient))), alpha=1, colour = NA) +
424. theme_article() +
425. theme(
426. panel.border = element_blank(),
427. legend.text = element_blank(),
428. axis.ticks = element_blank(),
429. axis.text.x = element_blank(),
430. axis.text.y = element_blank(),
431. axis.title.x = element_blank(),
432. axis.title.y = element_blank()) +
433. scale_x_continuous(breaks=seq(0, 21, by = 3),
434. labels = c("0" = "00:00", "3" = "", "6" = "06:00", "9" = "", "12" = "12:00", "15" = "", "18" = "18:00", "21" = "")) +
435. coord_polar(theta = "x", start = 0, clip="off") +
436. ylim(-2, 10)
438. polardelta2power <-
439. ggplot(data=data_binned[data_binned$band == "Delta2", ], aes(x = bin, y = Npower, fill=Patient, col=Patient)) +
440. scale_fill_brewer(palette = "Set2", direction = -1) + scale_color_brewer(palette = "Set2", direction = -1) +
441. geom_vline(xintercept = seq(0, 24, by = 3), colour = "grey90") +
442. geom_hline(yintercept = seq(0, 8, by = 1), colour = "grey90") +
443. geom_ribbon(aes(ymin = (9-as.numeric(Patient)), ymax=Npower*1.5+(9-as.numeric(Patient))), alpha=1, colour = NA) +
444. theme_article() +
445. theme(
446. panel.border = element_blank(),
447. legend.text = element_blank(),
448. axis.ticks = element_blank(),
449. axis.text.x = element_blank(),
450. axis.text.y = element_blank(),
451. axis.title.x = element_blank(),
452. axis.title.y = element_blank()) +
453. scale_x_continuous(breaks=seq(0, 21, by = 3),
454. labels = c("0" = "00:00", "3" = "", "6" = "06:00", "9" = "", "12" = "12:00", "15" = "", "18" = "18:00", "21" = "")) +
455. coord_polar(theta = "x", start = 0, clip="off") +
456. ylim(-1, 10)
459. polarpowerplots[[2]] <-
460. ggplot(data=data_binned[data_binned$band == "theta", ], aes(x = bin, y = Npower, fill=Patient, col=Patient)) +
461. scale_fill_brewer(palette = "Set2", direction = -1) + scale_color_brewer(palette = "Set2", direction = -1) +
462. geom_vline(xintercept = seq(0, 24, by = 3), colour = "grey90") +
463. geom_hline(yintercept = seq(0, 8, by = 1), colour = "grey90") +
464. geom_ribbon(aes(ymin = (9-as.numeric(Patient)), ymax=Npower*2+(9-as.numeric(Patient))), alpha=0.8, colour = NA) +
465. theme_article() +
466. theme(
467. panel.border = element_blank(),
468. legend.text = element_blank(),
469. axis.ticks = element_blank(),
470. axis.text.x = element_blank(),
471. axis.text.y = element_blank(),
472. axis.title.x = element_blank(),
473. axis.title.y = element_blank()) +
474. scale_x_continuous(breaks=seq(0, 21, by = 3),
475. labels = c("0" = "00:00", "3" = "", "6" = "06:00", "9" = "", "12" = "12:00", "15" = "", "18" = "18:00", "21" = "")) +
476. coord_polar(theta = "x", start = 0, clip="off") +
477. ylim(0, 10) +
478. facet_wrap(~title2)
480. polarpowerplots[[3]] <-
481. ggplot(data=data_binned[data_binned$band == "alpha", ], aes(x = bin, y = Npower, fill=Patient, col=Patient)) +
482. scale_fill_brewer(palette = "Set2", direction = -1) + scale_color_brewer(palette = "Set2", direction = -1) +
483. geom_vline(xintercept = seq(0, 24, by = 3), colour = "grey90") +
484. geom_hline(yintercept = seq(0, 8, by = 1), colour = "grey90") +
485. geom_ribbon(aes(ymin = (9-as.numeric(Patient)), ymax=Npower*2+(9-as.numeric(Patient))), alpha=0.8, colour = NA) +
486. theme_article() +
487. theme(
488. panel.border = element_blank(),
489. legend.text = element_blank(),
490. axis.ticks = element_blank(),
491. axis.text.x = element_blank(),
492. axis.text.y = element_blank(),
493. axis.title.x = element_blank(),
494. axis.title.y = element_blank()) +
495. scale_x_continuous(breaks=seq(0, 21, by = 3),
496. labels = c("0" = "00:00", "3" = "", "6" = "06:00", "9" = "", "12" = "12:00", "15" = "", "18" = "18:00", "21" = "")) +
497. coord_polar(theta = "x", start = 0, clip="off") +
498. ylim(0, 10) +
499. facet_wrap(~title3)
501. polarpowerplots[[4]] <-
502. ggplot(data=data_binned[data_binned$band == "delta_div_alpha", ], aes(x = bin, y = Npower, fill=Patient, col=Patient)) +
503. scale_fill_brewer(palette = "Set2", direction = -1) + scale_color_brewer(palette = "Set2", direction = -1) +
504. geom_vline(xintercept = seq(0, 24, by = 3), colour = "grey90") +
505. geom_hline(yintercept = seq(0, 8, by = 1), colour = "grey90") +
506. geom_ribbon(aes(ymin = (9-as.numeric(Patient)), ymax=Npower*2+(9-as.numeric(Patient))), alpha=0.8, colour = NA) +
507. theme_article() +
508. theme(
509. panel.border = element_blank(),
510. legend.text = element_blank(),
511. axis.ticks = element_blank(),
512. axis.text.y = element_blank(),
513. axis.text.x = element_blank(),
514. axis.title.x = element_blank(),
515. axis.title.y = element_blank()) +
516. scale_x_continuous(breaks=seq(0, 21, by = 3),
517. labels = c("0" = "00:00", "3" = "", "6" = "06:00", "9" = "", "12" = "12:00", "15" = "", "18" = "18:00", "21" = "")) +
518. coord_polar(theta = "x", start = 0, clip="off") +
519. ylim(0, 10) +
520. facet_wrap(~title4)
522. powerplots <- ggarrange(plotlist=polarpowerplots[c(1,2,3)], widths = c(1,1,1,1), heights = c(1,1,1,1), nrow = 1,
523. labels = c("D","E","F"), vjust = 18, hjust = -1,
524. legend = "right", common.legend = TRUE,
525. font.label = list(size = 14, color = "black", face = "bold"))
527. ggarrange(plotlist=polarpowerplots[c(1,2,3)], widths = c(1,1,1,1), heights = c(1,1,1,1), nrow = 1,
528. labels = c("A","B","C"), vjust = 18, hjust = -1,
529. legend = "right", common.legend = TRUE,
530. font.label = list(size = 14, color = "black", face = "bold")) %>%
531. ggexport(filename = "D:/Dropbox/Apps/Overleaf/Hspike/images/polar_power_band.pdf")
534. ######################
535. # IED rate circadian #
536. ######################
538. data_IED <- read.csv(file="Z:/analyses/stephen.whitmarsh/data/hspike/IED_table_PSG.txt", sep=',', header=TRUE, dec='.', na.strings = " ")
539. data_IED[data_IED$hyplabel == "PRE_SLEEP", ] = "AWAKE"
540. data_IED[data_IED$hyplabel == "POST_SLEEP", ] = "AWAKE"
541. data_IED$hyplabel <- factor(data_IED$hyplabel, ordered = TRUE, levels = c("REM", "AWAKE", "PHASE_1", "PHASE_2", "PHASE_3"))
543. data_IED$Patient <- factor(data_IED$patient, levels = c(8:1)) # same order in plot
544. data_IED$part <- factor(data_IED$part)
545. data_IED$marker <- factor(data_IED$marker)
546. data_IED$hour <- data_IED$minute / (60)
547. data_IED$rad <- data_IED$theta
549. # extract distribution statistics
550. dist_IED <- data.frame()
551. d <- list()
552. stat_IED <- list()
553. for (ipatient in 1:8) {
554. d <- density.circular(data_IED$rad[data_IED$patient == ipatient], bw = 100)
555. temp = list()
556. temp$x = as.numeric(d$x)
557. temp$y = d$y
558. le <- lengths(temp)
559. temp$patient <- rep(ipatient,le[1])
560. tempdf <- as.data.frame(temp)
561. colnames(tempdf) = c('rad','density','Patient')
562. dist_IED <- bind_rows(dist_IED, tempdf)
563. temp = list()
565. # Rayleigh Test of Uniformity: General Unimodal Alternative
566. temp <- rayleigh.test(data_IED$rad[data_IED$patient == ipatient])
567. stat_IED$Patient[ipatient] = ipatient
568. stat_IED$p[ipatient] = temp$p.value
569. stat_IED$Rayleigh_stat[ipatient] = temp$statistic
571. # Rayleigh Test of Uniformity: General Unimodal Alternative
572. stat_IED$median[ipatient] <- as.numeric(median(circular(data_IED$rad[data_IED$patient == ipatient]))) / (pi * 2) * 24
573. if (stat_IED$median[ipatient] < 0) {
574. stat_IED$median[ipatient] = stat_IED$median[ipatient] + 24
575. }
576. }
578. # format data for plotting
579. dist_IED$Khour = dist_IED$rad / (pi * 2) * 24
580. dist_IED$Patient = factor(dist_IED$Patient, levels = c(8:1)) # reversed order in plot
581. stat_IED <- as.data.frame(stat_IED)
582. stat_IED$Patient <- factor(stat_IED$Patient, levels = c(8:1))
584. # add some offset in degrees so it shows up from behind the rest
585. stat_IED$medianplot <- stat_IED$median
586. stat_IED$medianplot[1] = stat_IED$median[1] - 0.25
587. stat_IED$medianplot[6] = stat_IED$median[6] + 0.2
588. stat_IED$medianplot[stat_IED$p >= 0.05] = NA # all are significant though
590. IEDpolarplot <-
591. ggplot(data=data_IED, aes(x=hour, y = as.numeric(Patient), fill=Patient)) +
592. geom_vline(xintercept = seq(0, 21, by = 3), colour = "grey90") +
593. geom_hline(yintercept = seq(1, 8, by = 1), colour = "grey90") +
594. geom_ribbon(data=dist_IED, alpha = 1, colour = NA, aes(x = Khour,
595. ymin = (9-as.numeric(Patient)),
596. ymax = density*3 + (9-as.numeric(Patient)),
597. col = Patient), show.legend = TRUE) +
599. geom_segment(data=stat_IED, aes(x=medianplot, y=9.5, xend=medianplot, yend=10, col=Patient),
600. arrow = arrow(length = unit(0.25, "cm"), type="closed"), size = 0.5, show.legend = FALSE) +
602. coord_polar(theta = "x", start = 0, direction = 1, clip = 'off') +
603. scale_x_continuous(breaks = seq(0, 21, by = 3),
604. labels = c("0" = "00:00", "3" = "", "6" = "06:00", "9" = "", "12" = "12:00", "15" = "", "18" = "18:00", "21" = "")) +
605. scale_fill_brewer(palette = "Set2", direction=1) + scale_color_brewer(palette = "Set2", direction=1) +
606. theme_article() +
607. theme(panel.border = element_blank(),
608. #legend.key = element_blank(),
609. axis.ticks = element_blank(),
610. axis.text.y = element_blank(),
611. #axis.text.x = element_blank(),
612. panel.grid = element_blank(),
613. axis.title.x = element_blank(),
614. #legend.text = element_blank(),
615. axis.title.y = element_blank()) +
616. ylim(-2, 10)
618. # LaTeX table
619. library(stringr)
620. stat_IED$time = paste(str_pad( floor(stat_IED$median), 2, pad = "0"), ":", str_pad(floor((stat_IED$median- floor(stat_IED$median)) * 60), 2, pad = "0"), sep = "")
621. stat_IED <- stat_IED[, c("Patient", "time", "Rayleigh_stat", "p")]
622. stat_IED[,4] = ifelse(stat_IED[,4] > .05, paste(round(stat_IED[,4],digits=2),sep=""), ifelse(stat_IED[,4] < .0001, "<.0001\\textsuperscript{***}", ifelse(stat_IED[,4] < .001,"<.001\\textsuperscript{**}", ifelse(stat_IED[,4] < .01, "<.01\\textsuperscript{*}", "<.05"))))
624. kbl(stat_IED, "latex", booktabs = T, linesep = "", label = 'circstat_IED',
625. col.names = c("Patient","Median angle (HH:mm)","Rayleigh", "\\textit{p}"),
626. escape = FALSE, digits = 2,
627. caption = "Circular statistics of circadian IED rate") %>%
628. footnote(general_title = "",
629. footnote_as_chunk = TRUE,
630. threeparttable = TRUE,
631. escape = FALSE,
632. general = c("$. = p<0.05$, $* = p<0.01$, $** = p<0.001$, $*** = p<0.0001$"))%>%
633. kable_styling(latex_options = c("HOLD_position")) %>%
634. save_kable("D:/Dropbox/Apps/Overleaf/Neuronal mechanisms underlying the modulation of interictal epileptic activity during sleep/tables/circstat_IED.tex")
637. ######################
638. # Seizures circadian #
639. ######################
641. # load data: Patients x Units x time window
642. data_seizures <- read.csv(file="Z:/analyses/stephen.whitmarsh/data/hspike/seizuredata_table.txt", sep=',', header=TRUE, dec='.', na.strings = " ")
644. # prepare data
645. data_seizures$Patient <- factor(data_seizures$patient, levels = c(1:8))
646. data_seizures$hour <- data_seizures$minute / 60
647. data_seizures$rad <- data_seizures$minute / 60 / 24 * pi * 2
649. # extract distribution statistics
650. dist_seizures <- data.frame()
651. d <- list()
652. stat_seizures <- list()
653. for (ipatient in 1:8) {
654. d <- density.circular(data_seizures$rad[data_seizures$patient == ipatient], bw = 50)
655. temp = list()
656. temp$x = as.numeric(d$x)
657. temp$y = d$y
658. le <- lengths(temp)
659. temp$patient <- rep(ipatient,le[1])
660. tempdf <- as.data.frame(temp)
661. colnames(tempdf) = c('rad','density','Patient')
662. dist_seizures <- bind_rows(dist_seizures, tempdf)
664. # Rayleigh Test of Uniformity: General Unimodal Alternative
665. temp = list()
666. temp <- rayleigh.test(data_seizures$rad[data_seizures$patient == ipatient])
667. stat_seizures$Patient[ipatient] = ipatient
668. stat_seizures$p[ipatient] = temp$p.value
669. stat_seizures$Rayleigh_stat[ipatient] = temp$statistic
671. # Rayleigh Test of Uniformity: General Unimodal Alternative
672. stat_seizures$median[ipatient] <- as.numeric(median(circular(data_seizures$rad[data_seizures$patient == ipatient]))) / (pi * 2) * 24
673. if (stat_seizures$median[ipatient] < 0) {
674. stat_seizures$median[ipatient] = stat_seizures$median[ipatient] + 24
675. }
676. }
678. # format data for plotting
679. dist_seizures$Khour = dist_seizures$rad / (pi * 2) * 24
680. dist_seizures$Patient = factor(dist_seizures$Patient, levels = c(8:1))
681. stat_seizures <- as.data.frame(stat_seizures)
682. stat_seizures$Patient = factor(stat_seizures$Patient, levels = c(8:1))
683. stat_seizures$significant = stat_seizures$p < 0.05
684. stat_seizures$median_sel = stat_seizures$median
685. stat_seizures$median_sel[stat_seizures$p >= 0.05] = NA # all are significant though
687. # tiny adjustment to make axes line out properly and arrows not overlap
688. data_seizures$hourplot <- data_seizures$hour
689. data_seizures$hourplot[which(data_seizures$hour==max(data_seizures$hour))] = 24
691. # add jitter function for points
692. jitter <- position_jitter(width = 0, height = 0.4)
694. # plot
695. Seizurepolarplot <-
696. ggplot(data=data_seizures, aes(x=hourplot, y = as.numeric(Patient), fill=Patient)) +
697. geom_vline(xintercept = seq(0, 21, by = 3), colour = "grey90") +
698. geom_hline(yintercept = seq(1, 8, by = 1), colour = "grey90") +
699. geom_ribbon(data=dist_seizures, alpha = 1, colour = NA, aes(x=Khour,
700. ymin = (9-as.numeric(Patient)),
701. ymax = density*1.3 + (9-as.numeric(Patient)),
702. col = Patient), show.legend = FALSE) +
704. geom_segment(data=stat_seizures, aes(x=median_sel, y=9.5, xend=median_sel, yend=10, col=Patient),
705. arrow = arrow(length = unit(0.25, "cm"), type="closed"), size = 1, show.legend = FALSE) +
707. geom_point(colour="black", pch=21, size=1, position = jitter) +
708. coord_polar(theta = "x", start = 0, direction = 1, clip = 'off') +
709. scale_x_continuous(breaks = seq(0, 21, by = 3),
710. labels = c("0" = "00:00", "3" = "", "6" = "06:00", "9" = "", "12" = "12:00", "15" = "", "18" = "18:00", "21" = "")) +
711. scale_fill_brewer(palette = "Set2", direction = 1) + scale_color_brewer(palette = "Set2", direction = 1) +
712. theme_article() +
713. theme(panel.border = element_blank(),
714. #legend.key = element_blank(),
715. axis.ticks = element_blank(),
716. axis.text.y = element_blank(),
717. axis.text.x = element_blank(),
718. panel.grid = element_blank(),
719. axis.title.x = element_blank(),
720. #legend.text = element_blank(),
721. axis.title.y = element_blank()) +
722. ylim(-2, 10)
724. # save combined to pdf
725. # ggarrange(IEDpolarplot, Seizurepolarplot,
726. # labels = c("A","B"),
727. # vjust = 15, hjust = -1,
728. # legend = "right",
729. # common.legend = TRUE,
730. # font.label = list(size = 14, color = "black", face = "bold")) %>%
731. # ggexport(filename = "D:/Dropbox/Apps/Overleaf/Hspike/images/polar_density_seizures.pdf")
734. # save combined to pdf, with power
735. # ggarrange(IEDpolarplot, Seizurepolarplot, polardelta1power, polardelta2power,
736. # labels = c("A","B","C", "D"),
737. # vjust = 3, hjust = -1,
738. # legend = "right",
739. # common.legend = TRUE,
740. # font.label = list(size = 14, color = "black", face = "bold")) %>%
741. # ggexport(filename = "D:/Dropbox/Apps/Overleaf/Hspike/images/polar.pdf")
743. # save combined to pdf
744. ggarrange(IEDpolarplot, Seizurepolarplot,
745. labels = c("A","B"),
746. vjust = 15, hjust = -1,
747. ncol = 2, nrow = 1,
748. legend = "right",
749. common.legend = TRUE,
750. font.label = list(size = 14, color = "black", face = "bold")) %>%
751. ggexport(filename = "D:/Dropbox/Apps/Overleaf/Neuronal mechanisms underlying the modulation of interictal epileptic activity during sleep/polar.pdf")
753. # LaTeX table
754. library(stringr)
755. stat_seizures$time = paste(str_pad( floor(stat_seizures$median), 2, pad = "0"), ":", str_pad(floor((stat_seizures$median- floor(stat_seizures$median)) * 60), 2, pad = "0"), sep = "")
756. stat_seizures <- stat_seizures[, c("Patient", "time", "Rayleigh_stat", "p")]
757. stat_seizures[,4] = ifelse(stat_seizures[,4] > .05, paste(round(stat_seizures[,4],digits=2),sep=""), ifelse(stat_seizures[,4] < .0001, "<.0001\\textsuperscript{***}", ifelse(stat_seizures[,4] < .001,"<.001\\textsuperscript{**}", ifelse(stat_seizures[,4] < .01, "<.01\\textsuperscript{*}", "<.05"))))
759. kbl(stat_seizures, "latex", booktabs = T, linesep = "", label = 'circstat_seizures',
760. col.names = c("Patient","Median angle (HH:mm)","Rayleigh", "\\textit{p}"),
761. escape = FALSE, digits = 2,
762. caption = "Circular statistics of circadian seizure occurance")%>%
763. footnote(general_title = "",
764. footnote_as_chunk = TRUE,
765. threeparttable = TRUE,
766. escape = FALSE,
767. general = c("$. = p<0.05$, $* = p<0.01$, $** = p<0.001$, $*** = p<0.0001$"))%>%
768. kable_styling(latex_options = c("HOLD_position"))%>%
769. save_kable("D:/Dropbox/Apps/Overleaf/Neuronal mechanisms underlying the modulation of interictal epileptic activity during sleep/tables/circstat_seizures.tex")
771. # combined table
772. temp <- stat_seizures
773. colnames(temp) <- c("Patient","time2","Rayleigh_stat2","p2")
774. combined = merge(stat_IED,temp)
776. kbl(combined, "latex", booktabs = T, linesep = "", label = 'circstats',
777. col.names = c("Patient","Time","Rayleigh", "\\textit{p}","Time","Rayleigh", "\\textit{p}"),
778. escape = FALSE, digits = 2,
779. caption = "Circular statistics of circadian epileptic activity")%>%
780. add_header_above(c(" ", "Interictal activity" = 3, "Seizures" = 3)) %>%
781. kable_styling(latex_options = c("HOLD_position"))%>%
782. footnote(general_title = "",
783. footnote_as_chunk = TRUE,
784. threeparttable = TRUE,
785. escape = FALSE,
786. general = c("$.=p<0.05$, $*=p<0.01$, $**=p<0.001$, $***=p<0.0001$"))%>%
787. save_kable("D:/Dropbox/Apps/Overleaf/Neuronal mechanisms underlying the modulation of interictal epileptic activity during sleep/tables/circstats.tex")
790. #############################
791. # LFP power per sleep stage #
792. #############################
794. # prepare data
795. data_pow <- read.csv(file="Z:/analyses/stephen.whitmarsh/data/hspike/power_table_long.txt", sep=',', header=TRUE, dec='.', na.strings = " ")
796. data_pow <- data_pow[!data_pow$part > 3, ] # hypnogram is only scored on first three nights
797. data_pow <- data_pow[!data_pow$hyplabel == "NO_SCORE", ]
798. data_pow$hyplabel[data_pow$hyplabel == "PHASE_1"] = "S1"
799. data_pow$hyplabel[data_pow$hyplabel == "PHASE_2"] = "S2"
800. data_pow$hyplabel[data_pow$hyplabel == "PHASE_3"] = "S3"
801. data_pow$hyplabel[data_pow$hyplabel == "AWAKE"] = "WASO"
802. data_pow$hyplabel[data_pow$hyplabel == "PRESLEEP"] = "Pre"
803. data_pow$hyplabel[data_pow$hyplabel == "POSTSLEEP"] = "Post"
804. data_pow$hyplabel <- factor(data_pow$hyplabel, levels = c("Pre", "Post", "REM", "WASO", "S1", "S2", "S3"))
805. data_pow$band <- factor(data_pow$band, ordered = TRUE, levels = c("Delta1", "Delta2"))
806. data_pow$patient <- factor(data_pow$patient, levels = c(8:1))
807. data_pow$part <- factor(data_pow$part)
809. # relative to pre-sleep
810. temp <- setNames(aggregate(data_pow$power, by = list(data_pow$patient, data_pow$part, data_pow$band, data_pow$hyplabel), mean), c("patient", "part", "band", "hyplabel", "Pre"))
811. temp <- temp[temp$hyplabel=="Pre", ]
812. temp <- subset(temp, select = -c(hyplabel))
813. data_pow <- merge(data_pow, temp)
814. data_pow$Zpower <- (data_pow$power-data_pow$Pre) / (data_pow$power+data_pow$Pre)
815. data_pow$Zpower <- (data_pow$power/data_pow$Pre)
816. data_pow_sel <- data_pow[!data_pow$hyplabel=="Pre", ]
818. # averages for plotting
819. data_pow_avg <- setNames(aggregate(data_pow_sel$Zpower, by = list(data_pow_sel$patient, data_pow_sel$part, data_pow_sel$band, data_pow_sel$hyplabel), mean), c("patient", "part", "band", "hyplabel", "power_avg"))
820. data_pow_avg$patient <- factor(data_pow_avg$patient, levels = c(8:1))
821. data_pow_avg$part <- factor(data_pow_avg$part)
823. # plot
824. data_pow_sel$title1 = "SWA (0.1-2.5 Hz) power"
825. data_pow_sel$title2 = "Delta (2.5-4.0 Hz) power"
827. plot_delta1 <-
828. ggplot(data=data_pow_sel[data_pow_sel$band == "Delta1",], aes(y = hyplabel, x = Zpower)) +
829. geom_boxplot(outlier.shape = NA) +
830. scale_fill_brewer(palette = "Set2") + scale_color_brewer(palette = "Set2") +
831. geom_point(data = data_pow_avg[data_pow_avg$band == "Delta1",], aes(group = interaction(patient, part), x = power_avg, y = hyplabel, col = patient),
832. position=position_dodge(width=0.5)) +
833. guides(colour = "none") +
834. coord_cartesian(xlim = c(0, 20)) +
835. scale_x_continuous(breaks=c(0, 20)) +
836. scale_y_discrete(expand=expansion(mult=c(0.1,0.2))) +
837. theme_article() +
838. ylab(NULL) + xlab(NULL) +
839. facet_wrap(~title1)
841. plot_delta2 <-
842. ggplot(data=data_pow_sel[data_pow_sel$band == "Delta2",], aes(y = hyplabel, x = Zpower)) +
843. geom_boxplot(outlier.shape = NA) +
844. scale_fill_brewer(palette = "Set2") + scale_color_brewer(palette = "Set2") +
845. geom_point(data = data_pow_avg[data_pow_avg$band == "Delta2",], aes(group = interaction(patient, part), x = power_avg, y = hyplabel, col = patient),
846. position=position_dodge(width=0.5)) +
847. guides(colour = "none") +
848. coord_cartesian(xlim = c(0, 10)) +
849. scale_x_continuous(breaks=c(0, 10)) +
850. scale_y_discrete(expand=expansion(mult=c(0.1,0.2))) +
851. theme_article() +
852. ylab(NULL) + xlab(NULL) +
853. facet_wrap(~title2)
856. ##########################
857. # IED rate & sleep stage #
858. ##########################
860. # prepare data
861. data_IED <- read.csv(file="Z:/analyses/stephen.whitmarsh/data/hspike/IED_table_PSG.txt", sep=',', header=TRUE, dec='.', na.strings = " ")
862. data_IED <- data_IED[!data_IED$hyplabel == "NO_SCORE", ]
863. data_IED$hyplabel[data_IED$hyplabel == "PHASE_1"] = "S1"
864. data_IED$hyplabel[data_IED$hyplabel == "PHASE_2"] = "S2"
865. data_IED$hyplabel[data_IED$hyplabel == "PHASE_3"] = "S3"
866. data_IED$hyplabel[data_IED$hyplabel == "AWAKE"] = "WASO"
867. data_IED$hyplabel[data_IED$hyplabel == "PRESLEEP"] = "Pre"
868. data_IED$hyplabel[data_IED$hyplabel == "POSTSLEEP"] = "Post"
869. data_IED$hyplabel <- factor(data_IED$hyplabel, levels = c("Pre", "Post", "REM", "WASO", "S1", "S2", "S3"))
870. data_IED$Patient <- factor(data_IED$patient, levels = c(8:1))
871. data_IED$part <- factor(data_IED$part)
872. data_IED$marker <- factor(data_IED$marker)
873. data_IED$hour <- data_IED$minute / (60)
875. # normalize by time spend in sleep stages
876. data_duration <- read.csv(file="Z:/analyses/stephen.whitmarsh/data/hspike/hypnogram_duration.txt", sep=',', header=TRUE, dec='.', na.strings = " ")
877. data_duration <- data_duration %>% rename(
878. 'S1' = 'PHASE_1',
879. 'S2' = 'PHASE_2',
880. 'S3' = 'PHASE_3',
881. 'Pre' = 'PRESLEEP',
882. 'Post' = 'POSTSLEEP',
883. 'Total' = 'TOTAL')
885. data_duration <- melt(data_duration, id = c('patient','part'))
886. colnames(data_duration) = c('Patient','part','hyplabel','duration')
888. # count IEDs per sleepstage
889. data_IEDrate <- na.omit(data_IED %>% dplyr::count(Patient, part, hyplabel))
891. # normalize by time spend in sleep stages
892. data_IEDrate <- merge(data_IEDrate, data_duration)
894. data_IEDrate$IEDrate = data_IEDrate$n / data_IEDrate$duration / 60 # original rate is in Hz, now in minute
896. # normalize rate by Pre rate
897. i <- data_IEDrate[data_IEDrate$hyplabel == 'Pre',]
898. i <- i[, c('Patient','part','IEDrate')]
899. colnames(i) = c('Patient','part','Prerate')
900. data_IEDrate <- merge(data_IEDrate, i)
901. data_IEDrate$IEDrateNorm <- data_IEDrate$IEDrate - data_IEDrate$Prerate
902. data_IEDrate_sel <- data_IEDrate[!data_IEDrate$hyplabel=="Pre", ]
904. # plot
905. data_IEDrate_sel$title1 = "IED count"
906. data_IEDrate_sel$title2 = "IED rate (count/minute)"
907. data_IEDrate_sel$title3 = "IED rate - IED rate Pre-sleep (count/min)"
909. plot_IED_count <- ggplot(data=data_IEDrate_sel, aes(y=hyplabel, x=n)) +
910. scale_fill_brewer(palette = "Set2") + scale_color_brewer(palette = "Set2") +
911. geom_boxplot(outlier.shape = NA) +
912. geom_point(aes(group = interaction(Patient, part), col = Patient), position=position_dodge(width=0.5)) +
913. # guides(colour = "none") +
914. ylab(NULL) + xlab(NULL) +
915. theme_article() +
916. coord_cartesian(xlim = c(1, 3000)) +
917. scale_x_continuous(breaks=c(0, 3000)) +
918. scale_y_discrete(expand=expansion(mult=c(0.1,0.2))) +
919. theme(legend.position="right") +
920. facet_wrap(~title1)
922. plot_IED_rate <- ggplot(data=data_IEDrate_sel, aes(y=hyplabel, x=IEDrate)) +
923. scale_fill_brewer(palette = "Set2") + scale_color_brewer(palette = "Set2") +
924. geom_boxplot(outlier.shape = NA) +
925. geom_point(aes(group = interaction(Patient, part), col = Patient), position=position_dodge(width=0.5)) +
926. # guides(colour = "none") +
927. ylab(NULL) + xlab(NULL) +
928. theme_article() +
929. coord_cartesian(xlim = c(0, 20)) +
930. scale_x_continuous(breaks=c(0, 20)) +
931. scale_y_discrete(expand=expansion(mult=c(0.1,0.2))) +
932. theme(legend.position="right") +
933. facet_wrap(~title2)
935. plot_IED_norm <- ggplot(data=data_IEDrate_sel, aes(y=hyplabel, x=IEDrateNorm)) +
936. scale_fill_brewer(palette = "Set2") + scale_color_brewer(palette = "Set2") +
937. geom_boxplot(outlier.shape = NA) +
938. geom_point(aes(group = interaction(Patient, part), col = Patient), position=position_dodge(width=0.5)) +
939. ylab(NULL) + xlab(NULL) +
940. theme(axis.text.y = element_blank()) +
941. theme_article() +
942. coord_cartesian(xlim = c(0, 20)) +
943. scale_y_discrete(expand=expansion(mult=c(0.1,0.2))) +
944. scale_x_continuous(breaks=c(0, 20)) +
945. # theme(legend.position ="bottom") +
946. # guides(colour = guide_legend(ncol = 1)) +
947. # guides(fill=guide_legend(title="Patient")) +
948. theme(legend.position="right") +
949. labs(color='Patient') +
950. facet_wrap(~title3)
953. #####################################
954. # IED rate explained by sleep stage #
955. #####################################
957. # to create p-values
958. detach(package:lmerTest)
959. library(lmerTest)
960. library(lme4)
962. # determine reference level
963. data_IEDrate$hyplabel <- factor(data_IEDrate$hyplabel, levels = c("Pre","S3", "S2", "S1", "Wake", "REM", "Post"))
964. data_IEDrate$hyplabel = relevel(data_IEDrate$hyplabel, ref="Pre")
966. lIEDrate <- lmer(IEDrate ~ hyplabel + (1 | part) + (1 | Patient), data_IEDrate, control = lmerControl(optimizer ='Nelder_Mead'))
967. summary(lIEDrate)
968. plot_model(lIEDrate)
970. # Coefficients
971. temp = summary(lIEDrate)
972. coefs <- as.data.frame(temp$coefficients)
973. coefs[,5] = ifelse(coefs[,5] > .05, paste(round(coefs[,5],digits=2),sep=""), ifelse(coefs[,5] < .0001, "<.0001\\textsuperscript{***}", ifelse(coefs[,5] < .001,"<.001\\textsuperscript{**}", ifelse(coefs[,5] < .01, "<.01\\textsuperscript{*}", "<.05\\textsuperscript{.}"))))
974. rownames(coefs) <- c("\\textit{Intercept}", "S3", "S2", "S1", "Wake", "REM", "Post")
976. # Post-hoc tests
977. temp = emmeans(lIEDrate, list(pairwise ~ hyplabel), adjust = "tukey")
978. phIEDrate <- as.data.frame(temp$`pairwise differences of hyplabel`)
979. phIEDrate <- phIEDrate[, -4] # remove df since they are at inf
980. phIEDrate[,5] = ifelse(phIEDrate[,5] > .05, paste(round(phIEDrate[,5],digits=2),sep=""), ifelse(phIEDrate[,5] < .0001, "<.0001\\textsuperscript{***}", ifelse(phIEDrate[,5] < .001,"<.001\\textsuperscript{**}", ifelse(phIEDrate[,5] < .01, "<.01\\textsuperscript{*}", "<.05"))))
982. # Concatenate in one LaTeX table
983. coefs <- data.frame(Predictor = row.names(coefs), coefs);
984. rownames(coefs) <- NULL
985. colnames(phIEDrate) <- c("Predictor", "Coef $\\beta$","SE($\\beta$)","z", "\\textit{p}")
986. colnames(coefs) <- c("Predictor", "Coef $\\beta$","SE($\\beta$)", "df", "z","\\textit{p}")
987. stats_IEDrate <- bind_rows(coefs,phIEDrate)
989. options(knitr.kable.NA = '')
991. kbl(stats_IEDrate, "latex", booktabs = T, linesep = "", label = 'stats_IEDrate',
992. escape = FALSE, digits = 2,
993. caption = "Effect of sleep stages on Slow Wave activity (0.1-2.5Hz)")%>%
994. pack_rows("Coefficients", 1, 7) %>%
995. pack_rows("Post-hoc comparisons", 8, 28) %>%
996. kable_styling(latex_options = c("HOLD_position"))%>%
997. footnote(general_title = "",
998. footnote_as_chunk = TRUE,
999. threeparttable = TRUE,
1000. escape = FALSE,
1001. general = c("$. = p<0.05$, $* = p<0.01$, $** = p<0.001$, $*** = p<0.0001$"))%>%
1002. save_kable("D:/Dropbox/Apps/Overleaf/Neuronal mechanisms underlying the modulation of interictal epileptic activity during sleep/stats_IEDrate.tex")
1004. #################################
1005. # IED rate vs. power STATISTICS #
1006. #################################
1008. # prepare data
1009. data_pow_wide <- read.csv(file="Z:/analyses/stephen.whitmarsh/data/hspike/power_table_wide.txt", sep=',', header=TRUE, dec='.', na.strings = " ")
1010. data_pow_wide <- data_pow_wide[!data_pow_wide$hyplabel == "NO_SCORE", ]
1011. data_pow_wide$hyplabel[data_pow_wide$hyplabel == "PHASE_1"] = "S1"
1012. data_pow_wide$hyplabel[data_pow_wide$hyplabel == "PHASE_2"] = "S2"
1013. data_pow_wide$hyplabel[data_pow_wide$hyplabel == "PHASE_3"] = "S3"
1014. data_pow_wide$hyplabel[data_pow_wide$hyplabel == "AWAKE"] = "WASO"
1015. data_pow_wide$hyplabel[data_pow_wide$hyplabel == "PRESLEEP"] = "Pre"
1016. data_pow_wide$hyplabel[data_pow_wide$hyplabel == "POSTSLEEP"] = "Post"
1017. data_pow_wide$hyplabel <- factor(data_pow_wide$hyplabel, levels = c("Pre","S3", "S2", "S1", "WASO", "REM", "Post"))
1018. data_pow_wide$stage <- data_pow_wide$hyplabel
1019. data_pow_wide$Patient <- factor(data_pow_wide$patient, levels = c(8:1))
1020. data_pow_wide$night <- factor(data_pow_wide$part)
1022. # to create p-values
1023. detach(package:lmerTest)
1024. library(lmerTest)
1025. library(lme4)
1028. ###################################
1029. # Delta1 explained by sleep stage #
1030. ###################################
1032. lDelta1 <- lmer(Delta1 ~ stage + (1 | night) + (1 | Patient), data_pow_wide, control = lmerControl(optimizer ='Nelder_Mead'))
1033. summary(lDelta1)
1034. plot_model(lDelta1)
1036. # Coefficients
1037. temp = summary(lDelta1)
1038. coefs <- as.data.frame(temp$coefficients)
1039. coefs[,5] = ifelse(coefs[,5] > .05, paste(round(coefs[,5],digits=2),sep=""), ifelse(coefs[,5] < .0001, "<.0001\\textsuperscript{***}", ifelse(coefs[,5] < .001,"<.001\\textsuperscript{**}", ifelse(coefs[,5] < .01, "<.01\\textsuperscript{*}", "<.05\\textsuperscript{.}"))))
1040. rownames(coefs) <- c("\\textit{Intercept}", "S3", "S2", "S1", "WASO", "REM", "Post")
1042. # Post-hoc tests
1043. temp = emmeans(lDelta1, list(pairwise ~ stage), adjust = "tukey")
1044. phDelta1 <- as.data.frame(temp$`pairwise differences of stage`)
1045. phDelta1 <- phDelta1[, -4] # remove df since they are at inf
1046. phDelta1[,5] = ifelse(phDelta1[,5] > .05, paste(round(phDelta1[,5],digits=2),sep=""), ifelse(phDelta1[,5] < .0001, "<.0001\\textsuperscript{***}", ifelse(phDelta1[,5] < .001,"<.001\\textsuperscript{**}", ifelse(phDelta1[,5] < .01, "<.01\\textsuperscript{*}", "<.05"))))
1048. # Concatenate in one LaTeX table
1049. coefs <- data.frame(Predictor = row.names(coefs), coefs);
1050. rownames(coefs) <- NULL
1051. colnames(phDelta1) <- c("Predictor", "Coef $\\beta$","SE($\\beta$)","z", "\\textit{p}")
1052. colnames(coefs) <- c("Predictor", "Coef $\\beta$","SE($\\beta$)", "df", "z","\\textit{p}")
1053. stats_delta1 <- bind_rows(coefs,phDelta1)
1055. options(knitr.kable.NA = '')
1057. kbl(stats_delta1, "latex", booktabs = T, linesep = "", label = 'stats_delta1',
1058. escape = FALSE, digits = 2,
1059. caption = "Effect of sleep stage on Slow Wave (0.1-2.5Hz) power")%>%
1060. pack_rows("Coefficients", 1, 7) %>%
1061. pack_rows("Post-hoc comparisons", 8, 28) %>%
1062. kable_styling(latex_options = c("HOLD_position"))%>%
1063. footnote(general_title = "",
1064. footnote_as_chunk = TRUE,
1065. threeparttable = TRUE,
1066. escape = FALSE,
1067. general = c("$.=p<0.05$, $*=p<0.01$, $**=p<0.001$, $***=p<0.0001$"))%>%
1068. save_kable("D:/Dropbox/Apps/Overleaf/Neuronal mechanisms underlying the modulation of interictal epileptic activity during sleep/tables/stats_delta1.tex")
1071. ###################################
1072. # Delta2 explained by sleep stage #
1073. ###################################
1075. lDelta2 <- lmer(Delta2 ~ stage + (1 | night) + (1 | Patient), data_pow_wide, control = lmerControl(optimizer ='Nelder_Mead'))
1076. summary(lDelta2)
1077. plot_model(lDelta2)
1079. # Coefficients
1080. temp = summary(lDelta2)
1081. coefs <- as.data.frame(temp$coefficients)
1082. coefs[,5] = ifelse(coefs[,5] > .05, paste(round(coefs[,5],digits=2),sep=""), ifelse(coefs[,5] < .0001, "<.0001\\textsuperscript{***}", ifelse(coefs[,5] < .001,"<.001\\textsuperscript{**}", ifelse(coefs[,5] < .01, "<.01\\textsuperscript{*}", "<.05\\textsuperscript{.}"))))
1083. rownames(coefs) <- c("\\textit{Intercept}", "S3", "S2", "S1", "WASO", "REM", "Post")
1085. # Post-hoc tests
1086. temp = emmeans(lDelta2, list(pairwise ~ stage), adjust = "tukey")
1087. phDelta2 <- as.data.frame(temp$`pairwise differences of stage`)
1088. phDelta2 <- phDelta2[, -4] # remove df since they are at inf
1089. phDelta2[,5] = ifelse(phDelta2[,5] > .05, paste(round(phDelta2[,5],digits=2),sep=""), ifelse(phDelta2[,5] < .0001, "<.0001\\textsuperscript{***}", ifelse(phDelta2[,5] < .001,"<.001\\textsuperscript{**}", ifelse(phDelta2[,5] < .01, "<.01\\textsuperscript{*}", "<.05"))))
1091. # Concatenate in one LaTeX table
1092. coefs <- data.frame(Predictor = row.names(coefs), coefs);
1093. rownames(coefs) <- NULL
1094. colnames(phDelta2) <- c("Predictor", "Coef $\\beta$","SE($\\beta$)","z", "\\textit{p}")
1095. colnames(coefs) <- c("Predictor", "Coef $\\beta$","SE($\\beta$)", "df", "z","\\textit{p}")
1096. stats_delta2 <- bind_rows(coefs,phDelta2)
1098. options(knitr.kable.NA = '')
1100. kbl(stats_delta2, "latex", booktabs = T, linesep = "", label = 'stats_delta2',
1101. escape = FALSE, digits = 2,
1102. caption = "Effect of sleep stage on Delta (2.5-4 Hz) power")%>%
1103. pack_rows("Coefficients", 1, 7) %>%
1104. pack_rows("Post-hoc comparisons", 8, 28) %>%
1105. kable_styling(latex_options = c("HOLD_position"))%>%
1106. footnote(general_title = "",
1107. footnote_as_chunk = TRUE,
1108. threeparttable = TRUE,
1109. escape = FALSE,
1110. general = c("$.=p<0.05$, $*=p<0.01$, $**=p<0.001$, $***=p<0.0001$"))%>%
1111. save_kable("D:/Dropbox/Apps/Overleaf/Neuronal mechanisms underlying the modulation of interictal epileptic activity during sleep/tables/stats_delta2.tex")
1114. ####################################################
1115. # Mixed model with sleep stage explaining IED rate #
1116. ####################################################
1118. # to create p-values
1119. detach(package:lmerTest)
1120. library(lmerTest)
1121. library(lme4)
1123. # determine reference level
1124. data_pow_wide$hyplabel = relevel(data_pow_wide$stage, ref="Pre")
1125. l1 <- lmer(IEDsum ~ stage + (1 | night) + (1 | patient), data_pow_wide)
1127. summary(l1)
1128. plot_model(l1)
1130. # get mathematical description of the model and write to latex
1131. # eq <- equatiomatic::extract_eq(l1)
1132. # fileConn<-file("D:/Dropbox/Apps/Overleaf/Hspike/formula/model1.tex")
1133. # writeLines(c("$$",eq,"$$"), fileConn)
1134. # close(fileConn)
1136. # Coefficients to LaTeX table
1137. temp = summary(l1)
1138. coefs <- as.data.frame(temp$coefficients)
1139. coefs[,5] = ifelse(coefs[,5] > .05, paste(round(coefs[,5],digits=2),sep=""), ifelse(coefs[,5] < .0001, "<.0001\\textsuperscript{***}", ifelse(coefs[,5] < .001,"<.001\\textsuperscript{**}", ifelse(coefs[,5] < .01, "<.01\\textsuperscript{*}", "<.05\\textsuperscript{.}"))))
1140. rownames(coefs) <- c("\\textit{Intercept}", "S3", "S2", "S1", "WASO", "REM", "Post")
1142. # Post-hoc tests to LaTeX table
1143. temp = emmeans(l1, list(pairwise ~ stage), adjust = "tukey")
1144. ph1 <- as.data.frame(temp$`pairwise differences of stage`)
1145. ph1 <- ph1[, -4] # remove df since they are at inf
1146. ph1[,5] = ifelse(ph1[,5] > .05, paste(round(ph1[,5],digits=2),sep=""), ifelse(ph1[,5] < .0001, "<.0001\\textsuperscript{***}", ifelse(ph1[,5] < .001,"<.001\\textsuperscript{**}", ifelse(ph1[,5] < .01, "<.01\\textsuperscript{*}", "<.05"))))
1148. # Concatenate in one LaTeX table
1149. coefs <- data.frame(Predictor = row.names(coefs), coefs);
1150. rownames(coefs) <- NULL
1151. colnames(ph1) <- c("Predictor", "Coef $\\beta$","SE($\\beta$)","z", "\\textit{p}")
1152. colnames(coefs) <- c("Predictor", "Coef $\\beta$","SE($\\beta$)", "df", "z","\\textit{p}")
1153. stats_IEDsum <- bind_rows(coefs,ph1)
1155. kbl(stats_IEDsum, "latex", booktabs = T, linesep = "", label = 'stats_IEDsum',
1156. escape = FALSE, digits = 2,
1157. caption = "Effect of sleep stage on IEDs rate")%>%
1158. kable_styling(latex_options = c("HOLD_position"))%>%
1159. pack_rows("Sleep stages", 1, 7) %>% # latex_gap_space = "2em"
1160. pack_rows("Post-hoc comparisons", 8, 28) %>%
1161. footnote(general_title = "",
1162. footnote_as_chunk = TRUE,
1163. threeparttable = TRUE,
1164. escape = FALSE,
1165. general = c("$.=p<0.05$, $*=p<0.01$, $**=p<0.001$, $***=p<0.0001$"))%>%
1166. save_kable("D:/Dropbox/Apps/Overleaf/Neuronal mechanisms underlying the modulation of interictal epileptic activity during sleep/tables/stats_IEDsum.tex")
1169. #########################################
1170. # correlation between power and IEDrate #
1171. #########################################
1173. # prepare data
1174. data_pow_wide <- read.csv(file="Z:/analyses/stephen.whitmarsh/data/hspike/power_table_wide.txt", sep=',', header=TRUE, dec='.', na.strings = " ")
1175. data_pow_wide <- data_pow_wide[!data_pow_wide$hyplabel == "NO_SCORE", ]
1176. data_pow_wide$hyplabel[data_pow_wide$hyplabel == "PHASE_1"] = "S1"
1177. data_pow_wide$hyplabel[data_pow_wide$hyplabel == "PHASE_2"] = "S2"
1178. data_pow_wide$hyplabel[data_pow_wide$hyplabel == "PHASE_3"] = "S3"
1179. data_pow_wide$hyplabel[data_pow_wide$hyplabel == "AWAKE"] = "WASO"
1180. data_pow_wide$hyplabel[data_pow_wide$hyplabel == "PRESLEEP"] = "Pre"
1181. data_pow_wide$hyplabel[data_pow_wide$hyplabel == "POSTSLEEP"] = "Post"
1182. data_pow_wide$hyplabel <- factor(data_pow_wide$hyplabel, levels = c("Pre","S3", "S2", "S1", "WASE", "REM", "Post"))
1183. data_pow_wide$stage <- data_pow_wide$hyplabel
1184. data_pow_wide$Patient <- factor(data_pow_wide$patient)
1185. data_pow_wide$night <- factor(data_pow_wide$part)
1186. data_pow_wide$Delta1_log <- log(data_pow_wide$Delta1)
1187. data_pow_wide$Delta2_log <- log(data_pow_wide$Delta2)
1188. data_pow_wide$IEDsum_log <- log(data_pow_wide$IEDsum)
1189. data_pow_wide$Patient <- factor(data_pow_wide$patient, levels = c(8:1))
1191. # Delta 1
1192. data_pow_wide2 <- data_pow_wide %>% mutate(Delta_log_bin = cut(Delta1_log, breaks=seq(-1.5,10.5,1)))
1193. levels(data_pow_wide2$Delta_log_bin) <- seq(1:12)-2
1194. data_pow_wide2 <- data_pow_wide2 %>% group_by(Delta_log_bin, IEDsum)
1195. data_pow_wide2 <- data_pow_wide2 %>% summarise(count = n(), Patient)
1196. data_pow_wide2 <- data_pow_wide2[order(data_pow_wide2$Patient), ]
1198. # # Delta 1 for S3
1199. # data_pow_wide2 <- data_pow_wide[data_pow_wide$hyplabel == "S3", ] %>% mutate(Delta_log_bin = cut(Delta1_log, breaks=seq(-1.5,10.5,1)))
1200. # levels(data_pow_wide2$Delta_log_bin) <- seq(1:12)-2
1201. # data_pow_wide2 <- data_pow_wide2 %>% group_by(Delta_log_bin, IEDsum)
1202. # data_pow_wide2 <- data_pow_wide2 %>% summarise(count = n(), Patient)
1203. # data_pow_wide2 <- data_pow_wide2[order(data_pow_wide2$Patient), ]
1205. # ggplot(data=data_pow_wide2[data_pow_wide2$IEDsum > 0, ], aes(x=Delta1_log_bin, y=IEDsum*6)) +
1207. D1 <- ggplot(data=data_pow_wide2, aes(x=Delta_log_bin, y=IEDsum*6, color = Patient)) +
1208. geom_count() +
1209. # geom_count(aes(size = after_stat(prop), group = Delta_log_bin, color = Patient)) +
1210. # scale_size_area(max_size = 3) +
1211. # geom_smooth(data=data_pow_wide, aes(x=Delta1_log, y=IEDsum), method="lm", fullrange = FALSE, color="red") +
1212. scale_fill_brewer(palette = "Set2") + scale_color_brewer(palette = "Set2") +
1213. theme(axis.text.y = element_blank()) +
1214. theme_article() +
1215. theme(
1216. strip.background = element_blank(),
1217. strip.text.x = element_blank()
1218. ) +
1219. xlab("") +
1220. ylab("IEDs per minute") +
1221. # labs(size = "Proportion") +
1222. # coord_cartesian(xlim = c(2,11)) +
1223. # coord_cartesian(ylim = c(2,11)) +
1224. # facet_wrap(~Patient, ncol = 4, scales="free_x")
1225. facet_wrap(~Patient, ncol = 4)
1226. # ggsave("D:/Dropbox/Apps/Overleaf/Hspike/images/IEDrate_Delta1_count.pdf", width = 7, height = 4, units = "in")
1229. # Delta 2
1230. data_pow_wide2 <- data_pow_wide %>% mutate(Delta_log_bin = cut(Delta2_log, breaks=seq(-1.5,10.5,1)))
1231. levels(data_pow_wide2$Delta_log_bin) <- seq(1:12)-2
1232. data_pow_wide2 <- data_pow_wide2 %>% group_by(Delta_log_bin, IEDsum)
1233. data_pow_wide2 <- data_pow_wide2 %>% summarise(count = n(), Patient)
1234. data_pow_wide2 <- data_pow_wide2[order(data_pow_wide2$Patient), ]
1236. # ggplot(data=data_pow_wide2[data_pow_wide2$IEDsum > 0, ], aes(x=Delta1_log_bin, y=IEDsum*6)) +
1237. D2 <- ggplot(data=data_pow_wide2, aes(x=Delta_log_bin, y=IEDsum*6, color = Patient)) +
1238. # geom_count(aes(size = after_stat(prop), group = Delta_log_bin, color = Patient)) +
1239. geom_count() +
1240. # scale_size_area(max_size = 3) +
1241. # geom_smooth(data=data_pow_wide, aes(x=Delta1_log, y=IEDsum), method="lm", fullrange = FALSE, color="red") +
1242. scale_fill_brewer(palette = "Set2") + scale_color_brewer(palette = "Set2") +
1243. theme(axis.text.y = element_blank()) +
1244. theme_article() +
1245. theme(
1246. strip.background = element_blank(),
1247. strip.text.x = element_blank()
1248. ) +
1249. # xlab("Delta2 power (log)")
1250. #xlab = expression(mu "Volts" / "Hertz" ^ 2) +
1251. # xlab = expression("Force spaces with ~" ~ mu ~ pi * sigma ~ pi) +
1252. # ylab( units~are~(mu*g)/L )
1253. # xlab(TeX(r'($\alpha x^\alpha$, where $\alpha \in \{1 \ldots 5\}$)')) +
1254. # xlab(TeX(r'($log(\mu V/Hz^2$))')) +
1255. xlab(TeX(r'($log(\mu V^2/Hz$))')) +
1256. ylab("IEDs per minute") +
1257. labs(size = "Observations") +
1258. # coord_cartesian(xlim = c(2,11)) +
1259. # coord_cartesian(ylim = c(2,11)) +
1260. # facet_wrap(~Patient, ncol = 4, scales="free_x")
1261. facet_wrap(~Patient, ncol = 4)
1262. # ggsave("D:/Dropbox/Apps/Overleaf/Hspike/images/IEDrate_Delta2_count.pdf", width = 7, height = 4, units = "in")
1264. ggarrange(D1, D2,
1265. ncol = 1, nrow = 2,
1266. vjust = 1, hjust = 0,
1267. labels = c("A","B"),
1268. legend = "right",
1269. common.legend = TRUE,
1270. font.label = list(size = 14, color = "black", face = "bold")) %>%
1271. ggexport(filename = "D:/Dropbox/Apps/Overleaf/Neuronal mechanisms underlying the modulation of interictal epileptic activity during sleep/images/IEDrate_Delta_count.pdf")
1275. ###############################################
1276. # Mixed model with power explaining IED rate #
1277. ###############################################
1279. # to create p-values
1280. detach(package:lmerTest)
1281. library(lmerTest)
1282. library(lme4)
1284. # determine reference level
1285. data_pow_wide$hyplabel = relevel(data_pow_wide$stage, ref="Pre")
1286. l1 <- lmer(IEDsum ~ Delta1 + Delta2 + (1 | night) + (1 | patient), data_pow_wide)
1288. summary(l1)
1289. plot_model(l1)
1291. # # get mathematical description of the model and write to latex
1292. # eq <- equatiomatic::extract_eq(l1)
1293. # fileConn<-file("D:/Dropbox/Apps/Overleaf/Hspike/formula/model1.tex")
1294. # writeLines(c("$$",eq,"$$"), fileConn)
1295. # close(fileConn)
1297. # Coefficients to LaTeX table
1298. temp = summary(l1)
1299. coefs <- as.data.frame(temp$coefficients)
1300. coefs[,5] = ifelse(coefs[,5] > .05, paste(round(coefs[,5],digits=2),sep=""), ifelse(coefs[,5] < .0001, "<.0001\\textsuperscript{***}", ifelse(coefs[,5] < .001,"<.001\\textsuperscript{**}", ifelse(coefs[,5] < .01, "<.01\\textsuperscript{*}", "<.05\\textsuperscript{.}"))))
1301. rownames(coefs) <- c("\\textit{Intercept}", "Slow Wave", "Delta")
1303. kbl(coefs, "latex", booktabs = T, linesep = "", label = 'stats_IEDsum_vs_power',
1304. escape = FALSE, digits = 2,
1305. caption = "Effect of Slow Wave activity (0.1-2.5Hz) and Delta power (2.5-4Hz) on rate of IEDs")%>%
1306. kable_styling(latex_options = c("HOLD_position")) %>%
1307. footnote(general_title = "",
1308. footnote_as_chunk = TRUE,
1309. threeparttable = TRUE,
1310. escape = FALSE,
1311. general = c("$.=p<0.05$, $*=p<0.01$, $**=p<0.001$, $***=p<0.0001$"))%>%
1312. save_kable("D:/Dropbox/Apps/Overleaf/Neuronal mechanisms underlying the modulation of interictal epileptic activity during sleep/tables/stats_IEDsum_vs_power.tex")
1317. # ###################################################################
1318. # # Mixed model with both sleep stage and power explaining IED rate #
1319. # ###################################################################
1320. #
1321. # # to create p-values
1322. # detach(package:lmerTest)
1323. # library(lmerTest)
1324. # library(lme4)
1325. #
1326. # # determine reference level
1327. # data_pow_wide$hyplabel = relevel(data_pow_wide$stage, ref="Pre")
1328. # l1 <- lmer(IEDsum ~ stage + Delta1 + Delta2 + (1 | night) + (1 | patient), data_pow_wide)
1329. #
1330. # summary(l1)
1331. # plot_model(l1)
1332. #
1333. # # # get mathematical description of the model and write to latex
1334. # # eq <- equatiomatic::extract_eq(l1)
1335. # # fileConn<-file("D:/Dropbox/Apps/Overleaf/Hspike/formula/model1.tex")
1336. # # writeLines(c("$$",eq,"$$"), fileConn)
1337. # # close(fileConn)
1338. #
1339. # # Coefficients to LaTeX table
1340. # temp = summary(l1)
1341. # coefs <- as.data.frame(temp$coefficients)
1342. # coefs[,5] = ifelse(coefs[,5] > .05, paste(round(coefs[,5],digits=2),sep=""), ifelse(coefs[,5] < .0001, "<.0001\\textsuperscript{***}", ifelse(coefs[,5] < .001,"<.001\\textsuperscript{**}", ifelse(coefs[,5] < .01, "<.01\\textsuperscript{*}", "<.05\\textsuperscript{.}"))))
1343. # rownames(coefs) <- c("\\textit{Intercept}", "S3", "S2", "S1", "Wake", "REM", "Post", "SW", "Delta")
1344. #
1345. # # Post-hoc tests to LaTeX table
1346. # temp = emmeans(l1, list(pairwise ~ stage), adjust = "tukey")
1347. # ph1 <- as.data.frame(temp$`pairwise differences of stage`)
1348. # ph1 <- ph1[, -4] # remove df since they are at inf
1349. # ph1[,5] = ifelse(ph1[,5] > .05, paste(round(ph1[,5],digits=2),sep=""), ifelse(ph1[,5] < .0001, "<.0001\\textsuperscript{***}", ifelse(ph1[,5] < .001,"<.001\\textsuperscript{**}", ifelse(ph1[,5] < .01, "<.01\\textsuperscript{*}", "<.05"))))
1350. #
1351. # # Concatenate in one LaTeX table
1352. # coefs <- data.frame(Predictor = row.names(coefs), coefs);
1353. # rownames(coefs) <- NULL
1354. # colnames(ph1) <- c("Predictor", "Coef $\\beta$","SE($\\beta$)","z", "\\textit{p}")
1355. # colnames(coefs) <- c("Predictor", "Coef $\\beta$","SE($\\beta$)", "df", "z","\\textit{p}")
1356. # stats_IEDsum <- bind_rows(coefs,ph1)
1357. #
1358. # kbl(stats_IEDsum, "latex", booktabs = T, linesep = "", label = 'stats_IEDsum_power',
1359. # escape = FALSE, digits = 2,
1360. # caption = "Effect of sleepstages on IEDrate")%>%
1361. # kable_styling(latex_options = c("HOLD_position"))%>%
1362. # pack_rows("Sleep stages", 1, 7) %>% # latex_gap_space = "2em"
1363. # pack_rows("Power", 8, 9) %>%
1364. # pack_rows("Post-hoc comparisons", 10, 30) %>%
1365. # save_kable("D:/Dropbox/Apps/Overleaf/Hspike/tables/stats_IEDsum.tex")
1368. # Plot models
1370. set_theme(
1371. base = theme_article(),
1372. # panel.bordercol = NA
1373. )
1375. plot_IEDrate_model <- plot_model(
1376. lIEDrate,
1377. title = "",
1378. colors = "bw",
1379. axis.labels = "",
1380. axis.title = "",
1381. show.values = TRUE,
1382. show.p = TRUE,
1383. decimals = 4,
1384. digits = 2,
1385. value.offset = 0.5,
1386. value.size = 2.5) +
1387. ylim(-3, 9) +
1388. font_size(labels.x = 9, labels.y = 9) +
1389. scale_x_discrete(expand=expansion(mult=c(0.1,0.2)),
1390. labels=c("Post","REM","WASO","S1","S2","S3")) +
1391. scale_y_continuous(breaks=c(-2, 0, 8))
1392. plot_IEDrate_model$data$title = "IED rate (count/minute) model"
1393. plot_IEDrate_model <- plot_IEDrate_model + facet_wrap(~title, scales="free_y")
1395. plot_delta1_model <- plot_model(
1396. lDelta1,
1397. title = "",
1398. colors = "bw",
1399. axis.labels = "",
1400. axis.title = "",
1401. show.values = TRUE,
1402. show.p = TRUE,
1403. decimals = 4,
1404. digits = 2,
1405. value.offset = 0.5,
1406. value.size = 2.5) +
1407. ylim(0, 300) +
1408. font_size(labels.x = 9, labels.y = 9) +
1409. scale_x_discrete(expand=expansion(mult=c(0.1,0.2)),
1410. labels=c("Post","REM","WASO","S1","S2","S3")) +
1411. scale_y_continuous(expand=expansion(mult=c(0.1,0.1)), breaks=c(0, 270))
1412. plot_delta1_model$data$title = "SWA (0.1-2.5 Hz) model"
1413. plot_delta1_model <- plot_delta1_model + facet_wrap(~title, scales="free_y")
1415. plot_delta2_model <- plot_model(
1416. lDelta2,
1417. title = "",
1418. colors = "bw",
1419. axis.labels = "",
1420. axis.title = "",
1421. show.values = TRUE,
1422. show.p = TRUE,
1423. decimals = 4,
1424. digits = 2,
1425. value.offset = 0.5,
1426. value.size = 2.5) +
1427. ylim(-1, 41) +
1428. font_size(labels.x = 9, labels.y = 9) +
1429. scale_x_discrete(expand=expansion(mult=c(0.1,0.2)),
1430. labels=c("Post","REM","WASO","S1","S2","S3")) +
1431. scale_y_continuous(breaks=c(0, 40))
1432. plot_delta2_model$data$title = "Delta (2.5-4.0 Hz) model"
1433. plot_delta2_model <- plot_delta2_model + facet_wrap(~title, scales="free_y")
1436. # Boxplots together in Figure
1438. ggarrange(plot_IED_rate, plot_IEDrate_model, plot_delta1, plot_delta1_model, plot_delta2, plot_delta2_model,
1439. ncol = 2, nrow = 3,
1440. vjust = 1.5, hjust = -1,
1441. labels = c("A","B","C","D","E","F","G","H"),
1442. legend = "right",
1443. common.legend = TRUE,
1444. font.label = list(size = 14, color = "black", face = "bold")) %>%
1445. ggexport(filename = "D:/Dropbox/Apps/Overleaf/Neuronal mechanisms underlying the modulation of interictal epileptic activity during sleep/images/IEDrate_delta_boxplots.pdf")
1448. ###############################
1449. # IED amplitude & sleep stage #
1450. ###############################
1452. # prepare data
1453. data_amp <- read.csv(file="Z:/analyses/stephen.whitmarsh/data/hspike/amplitude_table.txt", sep=',', header=TRUE, dec='.', na.strings = " ")
1454. data_amp <- data_amp[!data_amp$hyplabel == "NO_SCORE", ]
1455. data_amp <- data_amp[!data_amp$hyplabel == "NO_SCORE", ]
1456. data_amp$hyplabel[data_amp$hyplabel == "PHASE_1"] = "S1"
1457. data_amp$hyplabel[data_amp$hyplabel == "PHASE_2"] = "S2"
1458. data_amp$hyplabel[data_amp$hyplabel == "PHASE_3"] = "S3"
1459. data_amp$hyplabel[data_amp$hyplabel == "AWAKE"] = "WASO"
1460. data_amp$hyplabel[data_amp$hyplabel == "PRESLEEP"] = "Pre"
1461. data_amp$hyplabel[data_amp$hyplabel == "POSTSLEEP"] = "Post"
1462. data_amp$stage <- factor(data_amp$hyplabel, levels = c("Pre","S3", "S2", "S1", "Wake", "REM", "Post"))
1463. data_amp$Patient <- factor(data_amp$patient, levels = c(8:1))
1464. data_amp$night <- factor(data_amp$part)
1465. data_amp$Template <- factor(data_amp$template)
1467. # # average over patient/part/sleepstage/template for plotting
1468. # posamp <- setNames(aggregate(data_amp$posamp, by = c(list(data_amp$patient, data_amp$night, data_amp$Template, data_amp$hyplabel)), mean), c("Patient", "part", "Template", "hyplabel", "posamp"))
1469. # negamp <- setNames(aggregate(data_amp$negamp, by = c(list(data_amp$patient, data_amp$night, data_amp$Template, data_amp$hyplabel)), mean), c("Patient", "part", "Template", "hyplabel", "negamp"))
1470. # data_amp_mean <- merge(posamp, negamp)
1472. # do not separate into templates
1473. posamp <- setNames(aggregate(data_amp$posamp, by = c(list(data_amp$patient, data_amp$night, data_amp$hyplabel)), mean), c("Patient", "night", "hyplabel", "posamp"))
1474. negamp <- setNames(aggregate(data_amp$negamp, by = c(list(data_amp$patient, data_amp$night, data_amp$hyplabel)), mean), c("Patient", "night", "hyplabel", "negamp"))
1475. data_amp_mean <- merge(posamp, negamp)
1478. # # relative to pre-sleep
1479. # temp_pos <- setNames(aggregate(data_amp$posamp, by = list(data_amp$patient, data_amp$part, data_amp$Template, data_amp$hyplabel), mean), c("Patient", "part", "Template", "hyplabel", "Pre_pos"))
1480. # temp_neg <- setNames(aggregate(data_amp$negamp, by = list(data_amp$patient, data_amp$part, data_amp$Template, data_amp$hyplabel), mean), c("Patient", "part", "Template", "hyplabel", "Pre_neg"))
1481. # temp_pos <- temp_pos[temp_pos$hyplabel=="Pre", ]
1482. # temp_neg <- temp_neg[temp_neg$hyplabel=="Pre", ]
1483. # temp_pos <- subset(temp_pos, select = -c(hyplabel))
1484. # temp_neg <- subset(temp_neg, select = -c(hyplabel))
1485. # data_amp_mean <- merge(data_amp_mean, temp_pos)
1486. # data_amp_mean <- merge(data_amp_mean, temp_neg)
1487. # data_amp_mean$Zposamp <- (data_amp_mean$posamp-data_amp_mean$Pre_pos)
1488. # data_amp_mean$Znegamp <- (data_amp_mean$negamp-data_amp_mean$Pre_neg)
1489. # data_amp_mean$diffamp <- data_amp_mean$Zposamp + data_amp_mean$Znegamp
1491. # # Normalized over all trials
1492. # y1 <- setNames(aggregate(data_amp_mean$posamp, by = c(list(data_amp_mean$Patient, data_amp_mean$Template)), mean), c("Patient", "Template", "Mposamp"))
1493. # y1sd <- setNames(aggregate(data_amp_mean$posamp, by = c(list(data_amp_mean$Patient, data_amp_mean$Template)), sd), c("Patient", "Template", "SDposamp"))
1494. # y2 <- setNames(aggregate(data_amp_mean$negamp, by = c(list(data_amp_mean$Patient, data_amp_mean$Template)), mean), c("Patient", "Template", "Mnegamp"))
1495. # y2sd <- setNames(aggregate(data_amp_mean$negamp, by = c(list(data_amp_mean$Patient, data_amp_mean$Template)), sd), c("Patient", "Template", "SDnegamp"))
1498. # do not separate into templates
1499. y1 <- setNames(aggregate(data_amp_mean$posamp, by = c(list(data_amp_mean$Patient)), mean), c("Patient", "Mposamp"))
1500. y1sd <- setNames(aggregate(data_amp_mean$posamp, by = c(list(data_amp_mean$Patient)), sd), c("Patient", "SDposamp"))
1501. y2 <- setNames(aggregate(data_amp_mean$negamp, by = c(list(data_amp_mean$Patient)), mean), c("Patient", "Mnegamp"))
1502. y2sd <- setNames(aggregate(data_amp_mean$negamp, by = c(list(data_amp_mean$Patient)), sd), c("Patient", "SDnegamp"))
1505. data_amp_mean <- merge(data_amp_mean,y1)
1506. data_amp_mean <- merge(data_amp_mean,y2)
1507. data_amp_mean <- merge(data_amp_mean,y1sd)
1508. data_amp_mean <- merge(data_amp_mean,y2sd)
1510. data_amp_mean$Zposamp <- (data_amp_mean$posamp-data_amp_mean$Mposamp)/data_amp_mean$SDposamp
1511. data_amp_mean$Znegamp <- (data_amp_mean$negamp-data_amp_mean$Mnegamp)/data_amp_mean$SDnegamp
1512. data_amp_mean$diffamp <- data_amp_mean$Zposamp - data_amp_mean$Znegamp
1514. data_amp_mean$Patient <- as.factor(data_amp_mean$Patient)
1515. data_amp_mean$hyplabel <- factor(data_amp_mean$hyplabel, levels = c("Pre","Post", "REM", "WASO", "S1", "S2", "S3"))
1517. data_amp_mean_sel <- data_amp_mean[!data_amp_mean$hyplabel=="Pre", ]
1519. # plot
1520. data_amp_mean_sel$title1 = "Standardized spike amplitude"
1521. data_amp_mean_sel$title2 = "Standardized slow-wave amplitude"
1522. data_amp_mean_sel$title3 = "Standardized spike vs. slow-wave amplitude"
1523. data_amp_mean_sel$Patient <- factor(data_amp_mean_sel$Patient, levels = c(8:1))
1525. plot_amp_pos <- ggplot(data=data_amp_mean_sel, aes(y=hyplabel, x=Zposamp)) +
1526. scale_fill_brewer(palette = "Set2") + scale_color_brewer(palette = "Set2") +
1527. geom_boxplot(outlier.shape = NA) +
1528. geom_point(aes(group = interaction(Patient, night), col = Patient), position=position_dodge(width=0.5), size = 1) +
1529. ylab(NULL) + xlab(NULL) +
1530. theme_article() +
1531. theme(legend.position="bottom") +
1532. coord_cartesian(xlim = c(-2.5, 2.5)) +
1533. facet_wrap(~title1)
1535. plot_amp_neg <- ggplot(data=data_amp_mean_sel, aes(y=hyplabel, x=Znegamp)) +
1536. scale_fill_brewer(palette = "Set2") + scale_color_brewer(palette = "Set2") +
1537. geom_boxplot(outlier.shape = NA) +
1538. geom_point(aes(group = interaction(Patient, night), col = Patient), position=position_dodge(width=0.5),size = 1) +
1539. ylab(NULL) + xlab(NULL) +
1540. theme_article() +
1541. theme(legend.position="bottom") +
1542. coord_cartesian(xlim = c(-3, 3)) +
1543. # coord_cartesian(xlim = c(-300, 100)) +
1544. facet_wrap(~title2)
1546. plot_amp_diff <- ggplot(data=data_amp_mean_sel, aes(y=hyplabel, x=diffamp)) +
1547. scale_fill_brewer(palette = "Set2") + scale_color_brewer(palette = "Set2") +
1548. geom_boxplot(outlier.shape = NA) +
1549. geom_point(aes(group = interaction(Patient, night), col = Patient), position=position_dodge(width=0.5), size = 1) +
1550. ylab(NULL) + xlab(NULL) +
1551. theme_article() +
1552. theme(legend.position="bottom") +
1553. coord_cartesian(xlim = c(-5, 5)) +
1554. scale_x_continuous(breaks=c(-5, 0, 5)) +
1555. facet_wrap(~title3)
1558. # plot pos vs neg in count scatter plot
1559. data_bin <- data_amp %>% group_by(Patient) %>% mutate(pos_bin = cut(posamp, 20))
1560. data_bin <- data_bin %>% group_by(Patient) %>% mutate(neg_bin = cut(negamp, 20))
1562. ggplot(data=data_bin, aes(y=pos_bin, x=neg_bin)) +
1563. geom_count(aes(color = Patient)) +
1564. scale_fill_brewer(palette = "Set2") + scale_color_brewer(palette = "Set2") +
1565. theme(axis.text.y = element_blank()) +
1566. theme_article() +
1567. theme(
1568. strip.background = element_blank(),
1569. strip.text.x = element_blank(),
1570. axis.text.x = element_blank(),
1571. axis.text.y = element_blank(),
1572. axis.ticks = element_blank()
1573. ) +
1574. ylab("Peak amplitude") +
1575. xlab("Slow wave amplitude") +
1576. # labs(size = "Proportion") +
1577. facet_wrap(~Patient, ncol = 4, scales = "free") +
1578. theme(aspect.ratio = 1)
1579. ggsave("D:/Dropbox/Apps/Overleaf/Hspike/images/amp_corr.pdf", width = 7, height = 4, units = "in")
1582. #########################
1583. ## STATISTICS LFP peaks #
1584. #########################
1586. # prepare data
1587. data_amp <- read.csv(file="Z:/analyses/stephen.whitmarsh/data/hspike/amplitude_table.txt", sep=',', header=TRUE, dec='.', na.strings = " ")
1588. data_amp <- data_amp[!data_amp$hyplabel == "NO_SCORE", ]
1589. data_amp <- data_amp[!data_amp$hyplabel == "NO_SCORE", ]
1590. data_amp$hyplabel[data_amp$hyplabel == "PHASE_1"] = "S1"
1591. data_amp$hyplabel[data_amp$hyplabel == "PHASE_2"] = "S2"
1592. data_amp$hyplabel[data_amp$hyplabel == "PHASE_3"] = "S3"
1593. data_amp$hyplabel[data_amp$hyplabel == "AWAKE"] = "WASO"
1594. data_amp$hyplabel[data_amp$hyplabel == "PRESLEEP"] = "Pre"
1595. data_amp$hyplabel[data_amp$hyplabel == "POSTSLEEP"] = "Post"
1596. data_amp$stage <- factor(data_amp$hyplabel, levels = c("Pre","S3", "S2", "S1", "WASO", "REM", "Post"))
1597. data_amp$Patient <- factor(data_amp$patient, levels = c(8:1))
1598. data_amp$night <- factor(data_amp$part)
1599. data_amp$template <- factor(data_amp$template)
1601. # Standardize data
1602. y1 <- setNames(aggregate(data_amp$posamp, by = c(list(data_amp$patient, data_amp$night, data_amp$template)), mean), c("patient", "night", "template", "Mposamp"))
1603. y2 <- setNames(aggregate(data_amp$negamp, by = c(list(data_amp$patient, data_amp$night, data_amp$template)), mean), c("patient", "night", "template", "Mnegamp"))
1604. y1sd <- setNames(aggregate(data_amp$posamp, by = c(list(data_amp$patient, data_amp$night, data_amp$template)), sd), c("patient", "night", "template", "SDposamp"))
1605. y2sd <- setNames(aggregate(data_amp$negamp, by = c(list(data_amp$patient, data_amp$night, data_amp$template)), sd), c("patient", "night", "template", "SDnegamp"))
1606. data_amp <- merge(data_amp,y1)
1607. data_amp <- merge(data_amp,y2)
1608. data_amp <- merge(data_amp,y1sd)
1609. data_amp <- merge(data_amp,y2sd)
1610. data_amp$Zposamp <- (data_amp$posamp-data_amp$Mposamp)/data_amp$SDposamp
1611. data_amp$Znegamp <- (data_amp$negamp-data_amp$Mnegamp)/data_amp$SDnegamp
1612. data_amp$Zdiffamp <- data_amp$Zposamp - data_amp$Znegamp
1614. # fit model
1615. data_amp$stage = relevel(data_amp$stage, ref="Pre")
1617. # to create p-values
1618. detach(package:lmerTest)
1619. library(lmerTest)
1620. library(lme4)
1622. lpos <- lmer(posamp ~ stage + (1 | patient) + (1 | night), data_amp, control = lmerControl(optimizer ='Nelder_Mead'))
1623. lneg <- lmer(negamp ~ stage + (1 | night) + (1 | template) + (1 | patient), data_amp, control = lmerControl(optimizer ='Nelder_Mead'))
1624. ldiff <- lmer(Zdiffamp ~ stage + (1 | night) + (1 | template) + (1 | patient), data_amp, control = lmerControl(optimizer ='Nelder_Mead'))
1626. summary(lpos)
1627. plot_model(lpos)
1628. summary(lneg)
1629. plot_model(lneg)
1630. summary(ldiff)
1631. plot_model(ldiff)
1633. # Post-hoc tests
1634. temp = emmeans(lpos, list(pairwise ~ stage), adjust = "tukey")
1635. phpos <- as.data.frame(temp$`pairwise differences of stage`)
1636. colnames(phpos) <- c("Comparison","Estimate","SE","df","Z ratio","p")
1637. phpos$df <- NA
1639. temp = emmeans(lneg, list(pairwise ~ stage), adjust = "tukey")
1640. phneg <- as.data.frame(temp$`pairwise differences of stage`)
1641. colnames(phneg) <- c("Comparison","EstimateNeg","SENeg","dfNeg","Z ratioNeg","pNeg")
1642. phneg$dfNeg <- NA
1644. temp = emmeans(ldiff, list(pairwise ~ stage), adjust = "tukey")
1645. phdiff <- as.data.frame(temp$`pairwise differences of stage`)
1646. colnames(phdiff) <- c("Comparison","EstimateDiff","SEDiff","dfDiff","Z ratioDiff","pDiff")
1647. phdiff$dfDiff <- NA
1649. # Mathematical description of the model and write to latex
1650. # eq <- equatiomatic::extract_eq(lpos)
1651. # fileConn<-file("D:/Dropbox/Apps/Overleaf/Hspike/formula/model_posamp.tex")
1652. # writeLines(c("$$",eq,"$$"), fileConn)
1653. # close(fileConn)
1656. #############################################
1657. # Coefficients to LaTeX table (and reorder) #
1658. #############################################
1660. # Model coefficients for LaTeX table
1661. temp = summary(lpos)
1662. spos <- temp$coefficients
1663. spos <- data.frame(Predictor = row.names(spos), spos);
1664. rownames(spos) <- NULL
1665. colnames(spos) <- c("Predictor","Estimate","SD","df","t","p")
1667. temp = summary(lneg)
1668. sneg <- temp$coefficients
1669. sneg <- data.frame(Predictor = row.names(sneg), sneg)
1670. colnames(sneg) <- c("Predictor","EstimateNeg","SDNeg","dfNeg","tNeg","pNeg")
1672. temp = summary(lneg)
1673. sdiff <- temp$coefficients
1674. sdiff <- data.frame(Predictor = row.names(sdiff), sdiff)
1675. colnames(sdiff) <- c("Predictor","EstimatePos","SDDiff","dfDiff","tDiff","pDiff")
1677. sneg$id <- 1:nrow(sneg)
1678. coef <- merge(spos, sneg)
1679. coef <- merge(coef, sdiff)
1680. coef <- coef[order(coef$id), ]
1681. coef <- coef[, c(-1, -12)]
1682. rownames(coef) <- NULL
1684. coef[, 3] = round(coef[,3],digits=0)
1685. coef[, 8] = round(coef[,8],digits=0)
1686. coef[,13] = round(coef[,8],digits=0)
1688. coef$Predictor <- c("\\textit{Intercept}","S3", "S2", "S1", "WASO", "Post", "REM")
1689. coef <- coef[, c(16,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15)]
1691. # Post-hoc comparisons
1692. phneg$id <- 1:nrow(phneg)
1693. ph <- merge(phpos,phneg)
1694. ph <- merge(ph,phdiff)
1695. ph <- ph[order(ph$id), ]
1696. ph <- ph[, -12]
1697. rownames(ph) <- NULL
1699. # Concatenate in one LaTeX table
1700. colnames(ph) <- c("Predictor", "Coef $\\beta$","SE($\\beta$)","df", "z", "\\textit{p}",
1701. "Coef $\\beta$","SE($\\beta$)","df","z", "\\textit{p}",
1702. "Coef $\\beta$","SE($\\beta$)","df","z", "\\textit{p}")
1703. colnames(coef) <- c("Predictor", "Coef $\\beta$","SE($\\beta$)","df", "z", "\\textit{p}",
1704. "Coef $\\beta$","SE($\\beta$)","df","z", "\\textit{p}",
1705. "Coef $\\beta$","SE($\\beta$)","df","z", "\\textit{p}")
1707. stats_amp <- bind_rows(coef,ph)
1708. colnames(stats_amp) <- c("", "Coef $\\beta$","SE($\\beta$)","df", "z", "\\textit{p}",
1709. "Coef $\\beta$","SE($\\beta$)","df","z", "\\textit{p}",
1710. "Coef $\\beta$","SE($\\beta$)","df","z", "\\textit{p}")
1712. stats_amp[,6] = ifelse(stats_amp[,6] > .05, paste(round(stats_amp[,6],digits=2),sep=""),
1713. ifelse(stats_amp[,6] < .0001, "<.0001\\textsuperscript{***}",
1714. ifelse(stats_amp[,6] < .001,"<.001\\textsuperscript{**}",
1715. ifelse(stats_amp[,6] < .01, "<.01\\textsuperscript{*}", "<.05\\textsuperscript{.}"))))
1716. stats_amp[,11] = ifelse(stats_amp[,11] > .05, paste(round(stats_amp[,11],digits=2),sep=""),
1717. ifelse(stats_amp[,11] < .0001, "<.0001\\textsuperscript{***}",
1718. ifelse(stats_amp[,11] < .001,"<.001\\textsuperscript{**}",
1719. ifelse(stats_amp[,11] < .01, "<.01\\textsuperscript{*}", "<.05\\textsuperscript{.}"))))
1720. stats_amp[,16] = ifelse(stats_amp[,16] > .05, paste(round(stats_amp[,16],digits=2),sep=""),
1721. ifelse(stats_amp[,16] < .0001, "<.0001\\textsuperscript{***}",
1722. ifelse(stats_amp[,16] < .001,"<.001\\textsuperscript{**}",
1723. ifelse(stats_amp[,16] < .01, "<.01\\textsuperscript{*}", "<.05\\textsuperscript{.}"))))
1725. options(knitr.kable.NA = '')
1728. kbl(stats_amp, "latex", booktabs = T, linesep = "", label = 'stats_amp',
1729. escape = FALSE, digits = 2,
1730. caption = "Effect of sleep stage on ERP peak amplitude")%>%
1731. kable_styling(latex_options = c("HOLD_position"))%>%
1732. pack_rows("Sleep stages", 1, 7) %>% # latex_gap_space = "2em"
1733. pack_rows("Post-hoc comparisons", 8, 28) %>%
1734. add_header_above(c(" ", "Spike amplitude" = 5, "Slow wave amplitude" = 5, "Difference amplitude" = 5)) %>%
1735. kable_styling(latex_options = c("scale_down"))%>%
1736. footnote(general_title = "",
1737. footnote_as_chunk = TRUE,
1738. threeparttable = TRUE,
1739. escape = FALSE,
1740. general = c("$.=p<0.05$, $*=p<0.01$, $**=p<0.001$, $***=p<0.0001$"))%>%
1741. save_kable("D:/Dropbox/Apps/Overleaf/Neuronal mechanisms underlying the modulation of interictal epileptic activity during sleep/tables/stats_amp.tex")
1744. #################
1745. ## Plot models ##
1746. #################
1748. set_theme(
1749. base = theme_article(),
1750. # panel.bordercol = NA
1751. )
1753. plot_amp_pos_model <- plot_model(
1754. lpos,
1755. title = "",
1756. colors = "bw",
1757. axis.labels = "",
1758. axis.title = "",
1759. show.values = TRUE,
1760. show.p = TRUE,
1761. decimals = 4,
1762. digits = 2,
1763. value.offset = 0.5,
1764. value.size = 2.5) +
1765. font_size(labels.x = 9, labels.y = 9) +
1766. scale_x_discrete(expand=expansion(mult=c(0.1,0.2)),
1767. labels=c("Post","REM","WASO","S1","S2","S3")) +
1768. scale_y_continuous(breaks=c(-20, 0, 50, 100, 150))
1769. plot_amp_pos_model$data$title = "Model fixed effects"
1770. plot_amp_pos_model <- plot_amp_pos_model + facet_wrap(~title, scales="free_y")
1772. plot_amp_neg_model <- plot_model(
1773. lneg,
1774. title = "",
1775. colors = "bw",
1776. axis.labels = "",
1777. axis.title = "",
1778. show.values = TRUE,
1779. show.p = TRUE,
1780. decimals = 4,
1781. digits = 2,
1782. value.offset = 0.5,
1783. value.size = 2.5) +
1784. font_size(labels.x = 9, labels.y = 9) +
1785. scale_x_discrete(expand=expansion(mult=c(0.1,0.2)),
1786. labels=c("Post","REM","WASO","S1","S2","S3")) +
1787. scale_y_continuous(breaks=c(0, 50, 100, 150))
1788. plot_amp_neg_model$data$title = "Model fixed effects"
1789. plot_amp_neg_model <- plot_amp_neg_model + facet_wrap(~title, scales="free_y")
1791. plot_amp_diff_model <- plot_model(
1792. ldiff,
1793. title = "",
1794. colors = "bw",
1795. axis.labels = "",
1796. axis.title = "",
1797. show.values = TRUE,
1798. show.p = TRUE,
1799. decimals = 4,
1800. digits = 2,
1801. value.offset = 0.5,
1802. value.size = 2.5) +
1803. font_size(labels.x = 9, labels.y = 9) +
1804. scale_x_discrete(expand=expansion(mult=c(0.1,0.2)),
1805. labels=c("Post","REM","WASO","S1","S2","S3")) +
1806. scale_y_continuous(breaks=c(-0.4, -0.2, 0, 0.2, 0.4))
1807. plot_amp_diff_model$data$title = "Model fixed effects"
1808. plot_amp_diff_model <- plot_amp_diff_model + facet_wrap(~title, scales="free_y")
1811. ###############################
1812. # Boxplots together in Figure #
1813. ###############################
1815. ggarrange(plot_amp_pos, plot_amp_pos_model, plot_amp_neg, plot_amp_neg_model, plot_amp_diff, plot_amp_diff_model,
1816. ncol = 2, nrow = 3,
1817. vjust = 1.5, hjust = -1,
1818. labels = c("A","B","C","D","E","F","G","H"),
1819. legend = "right",
1820. common.legend = TRUE,
1821. font.label = list(size = 14, color = "black", face = "bold")) %>%
1822. ggexport(filename = "D:/Dropbox/Apps/Overleaf/Neuronal mechanisms underlying the modulation of interictal epileptic activity during sleep/images/amp_boxplots.pdf")
1825. ######################################
1826. # Plot normalized PSTH & sleep stage #
1827. ######################################
1829. # prepare data
1830. data_psth <- read.csv(file="Z:/analyses/stephen.whitmarsh/data/hspike/psth_table.txt", sep=',', header=TRUE, dec='.', na.strings = " ")
1831. data_psth$hyplabel[data_psth$hyplabel == "Wake"] = "WASO"
1833. data_psth$hyplabel <- factor(data_psth$hyplabel, levels = c("Pre", "Post", "REM", "WASO", "S1", "S2", "S3"))
1834. data_psth$Patient <- factor(data_psth$Patient, levels = c(8:1))
1835. data_psth$part <- factor(data_psth$part)
1836. data_psth$template <- factor(data_psth$template)
1837. data_psth$unit <- factor(data_psth$unit)
1838. data_psth$Type <- factor(data_psth$SUA)
1840. # data_sel <- data_psth[data_psth$hyplabel=="S1", ]
1841. data_sel <- setNames(aggregate(data_psth$responsive, by = c(list(data_psth$Patient, data_psth$part, data_psth$unit, data_psth$Type)), mean), c("Patient", "part", "unit", "Type", "responsive"))
1842. data_sel$responsive = as.integer(data_sel$responsive > 0)
1844. t = data_sel %>% count(Patient, part, responsive, Type)
1845. t2 <- pivot_wider(t, names_from = "Type", names_prefix = "SUA", values_from = "n")
1846. t3 <- pivot_wider(t2, names_from = "responsive", names_prefix = "responsive", values_from = c("SUA1", "SUA0"))
1847. # t4 <- pivot_wider(t, names_from = c("responsive", "Type"), values_from = c("n"))
1849. #
1850. # ############# USE THISFOR TABLE OF SUA/MUA reporting ###############
1851. #
1852. #
1853. # test <- setNames(aggregate(data_psth$responsive, by = c(list(data_psth$Patient, data_psth$unit, data_psth$template, data_psth$Type)), mean), c("Patient", "unit", "template", "Type", "responsive"))
1854. #
1855. # # select responsive units
1856. # data_psth <- data_psth[c(data_psth$responsive == 1), ]
1857. #
1858. # # remove unresponsive patient - already done in MATLAB
1859. # # data_psth <- data_psth[-c(data_psth$Patient == 7), ]
1860. #
1861. # # select SUA
1862. # # data_psth <- data_psth[c(data_psth$SUA == 1), ]
1864. posrate <- setNames(aggregate(data_psth$posrate, by = c(list(data_psth$Patient, data_psth$unit, data_psth$template, data_psth$hyplabel, data_psth$Type)), mean), c("Patient", "unit", "template", "hyplabel", "Type", "posrate"))
1865. negrate <- setNames(aggregate(data_psth$negrate, by = c(list(data_psth$Patient, data_psth$unit, data_psth$template, data_psth$hyplabel, data_psth$Type)), mean), c("Patient", "unit", "template", "hyplabel", "Type", "negrate"))
1866. data_psth_mean <- merge(posrate, negrate)
1868. y1 <- setNames(aggregate(data_psth_mean$posrate, by = c(list(data_psth_mean$Patient, data_psth_mean$unit, data_psth_mean$template, data_psth_mean$Type)), mean), c("Patient", "unit", "template", "Type", "Mposrate"))
1869. y1sd <- setNames(aggregate(data_psth_mean$posrate, by = c(list(data_psth_mean$Patient, data_psth_mean$unit, data_psth_mean$template, data_psth_mean$Type)), sd), c("Patient", "unit", "template", "Type", "SDposrate"))
1870. y2 <- setNames(aggregate(data_psth_mean$negrate, by = c(list(data_psth_mean$Patient, data_psth_mean$unit, data_psth_mean$template, data_psth_mean$Type)), mean), c("Patient", "unit", "template", "Type", "Mnegrate"))
1871. y2sd <- setNames(aggregate(data_psth_mean$negrate, by = c(list(data_psth_mean$Patient, data_psth_mean$unit, data_psth_mean$template, data_psth_mean$Type)), sd), c("Patient", "unit", "template", "Type", "SDnegrate"))
1873. data_psth_mean <- merge(data_psth_mean,y1)
1874. data_psth_mean <- merge(data_psth_mean,y2)
1875. data_psth_mean <- merge(data_psth_mean,y1sd)
1876. data_psth_mean <- merge(data_psth_mean,y2sd)
1878. data_psth_mean$Zposrate <- (data_psth_mean$posrate-data_psth_mean$Mposrate)/data_psth_mean$SDposrate
1879. data_psth_mean$Znegrate <- (data_psth_mean$negrate-data_psth_mean$Mnegrate)/data_psth_mean$SDnegrate
1880. data_psth_mean$diffrate <- data_psth_mean$Zposrate - data_psth_mean$Znegrate
1882. # plot
1883. data_psth_mean$title1 = "Standardized firingrate during spike"
1884. data_psth_mean$title2 = "Standardized firingrate during slow wave"
1885. data_psth_mean$title3 = "Relative difference"
1887. data_psth_mean_sel <- data_psth_mean[!data_psth_mean$hyplabel=="Pre", ]
1890. plot_cnt_pos <- ggplot(data=data_psth_mean_sel, aes(y=hyplabel, x=Zposrate)) +
1891. scale_fill_brewer(palette = "Set2") + scale_color_brewer(palette = "Set2") +
1892. geom_boxplot(outlier.shape = NA) +
1893. geom_point(aes(group = interaction(Patient, unit), col = Patient, shape = Type), position=position_dodge(width=0.5), alpha = 0.2, size = 0.8) +
1894. geom_point(data=data_psth_mean_sel[data_psth_mean_sel$Type == 1, ],
1895. aes(group = interaction(Patient, unit), col = Patient, shape = Type), position=position_dodge(width=0.5), size = 0.8) +
1896. scale_shape_discrete(label = c("MUA", "SUA"))+
1897. ylab(NULL) + xlab(NULL) +
1898. theme_article() +
1899. theme(legend.position="bottom") +
1900. coord_cartesian(xlim = c(-2.5, 2.5)) +
1901. scale_x_continuous(breaks = c(-2.5, 0, 2.5)) +
1902. scale_y_discrete(expand=expansion(mult=c(0.1,0.2))) +
1903. facet_wrap(~title1)
1905. plot_cnt_neg <- ggplot(data=data_psth_mean_sel, aes(y=hyplabel, x=Znegrate)) +
1906. scale_fill_brewer(palette = "Set2") + scale_color_brewer(palette = "Set2") +
1907. geom_boxplot(outlier.shape = NA) +
1908. geom_point(aes(group = interaction(Patient, unit), col = Patient, shape = Type), position=position_dodge(width=0.5), alpha = 0.2, size = 0.8) +
1909. geom_point(data=data_psth_mean_sel[data_psth_mean_sel$Type == 1, ],
1910. aes(group = interaction(Patient, unit), col = Patient, shape = Type), position=position_dodge(width=0.5), size = 0.8) +
1911. scale_shape_discrete(label = c("MUA", "SUA"))+
1912. ylab(NULL) + xlab(NULL) +
1913. theme_article() +
1914. theme(legend.position="bottom") +
1915. coord_cartesian(xlim = c(-2.5, 2.5)) +
1916. scale_x_continuous(breaks = c(-2.5, 0, 2.5)) +
1917. scale_y_discrete(expand=expansion(mult=c(0.1,0.2))) +
1918. facet_wrap(~title2)
1920. plot_cnt_diff <- ggplot(data=data_psth_mean_sel, aes(y=hyplabel, x=diffrate)) +
1921. scale_fill_brewer(palette = "Set2") + scale_color_brewer(palette = "Set2") +
1922. geom_boxplot(outlier.shape = NA) +
1923. geom_point(aes(group = interaction(Patient, unit), col = Patient, shape = Type), position=position_dodge(width=0.5), alpha = 0.2, size = 0.8) +
1924. geom_point(data=data_psth_mean_sel[data_psth_mean_sel$Type == 1, ],
1925. aes(group = interaction(Patient, unit), col = Patient, shape = Type), position=position_dodge(width=0.5), size = 0.8) +
1926. scale_shape_discrete(label = c("MUA", "SUA"))+
1927. ylab(NULL) + xlab(NULL) +
1928. theme_article() +
1929. theme(legend.position="bottom") +
1930. scale_y_discrete(expand=expansion(mult=c(0.1,0.2))) +
1931. coord_cartesian(xlim = c(-4, 4)) +
1932. scale_x_continuous(breaks = c(-4, 0, 4)) +
1933. facet_wrap(~title3)
1936. ###################################
1937. ## STATISTICS Positive peaks PSTH #
1938. ###################################
1940. # to create p-values
1941. detach(package:lmerTest)
1942. library(lmerTest)
1943. library(lme4)
1945. # prepare data
1946. data_psth <- read.csv(file="Z:/analyses/stephen.whitmarsh/data/hspike/psth_table.txt", sep=',', header=TRUE, dec='.', na.strings = " ")
1947. data_psth$hyplabel[data_psth$hyplabel == "Wake"] = "WASO"
1948. data_psth$hyplabel <- factor(data_psth$hyplabel, levels = c("Pre", "Post", "REM", "WASO", "S1", "S2", "S3"))
1949. data_psth$Patient <- factor(data_psth$Patient, levels = c(8:1))
1950. data_psth$part <- factor(data_psth$part)
1952. data_psth$template <- factor(as.integer(data_psth$template) + as.integer(data_psth$part) * 100 + as.integer(data_psth$Patient) * 1000) # beacuse template nr. resets for each night/patient
1953. data_psth$unit <- factor(data_psth$unit)
1954. data_psth$responsive <- factor(data_psth$responsive)
1955. data_psth$diffrate <- data_psth$posrate - data_psth$negrate
1957. # select responsive units
1958. data_psth <- data_psth[c(data_psth$responsive == 1), ]
1960. # remove unresponsive patient
1961. # data_psth <- data_psth[-c(data_psth$Patient == 7), ]
1963. # remove NaN
1964. data_psth <- data_psth[!is.nan(data_psth$posrate), ]
1966. # Rename for plotting
1967. data_psth$night <- factor(data_psth$part)
1969. # Reorder for table
1970. data_psth$stage <- factor(data_psth$hyplabel, levels = c("Pre","S3", "S2", "S1", "WASO", "REM", "Post"))
1972. # fit model
1973. data_psth$stage = relevel(data_psth$stage, ref="Pre")
1975. lpos_SUA <- lmer(posrate ~ stage + (1 | Patient) + (1 | template) + (1 | unit), data_psth[data_psth$SUA == 1, ], control = lmerControl(optimizer ='Nelder_Mead'))
1976. lneg_SUA <- lmer(negrate ~ stage + (1 | Patient) + (1 | template) + (1 | unit), data_psth[data_psth$SUA == 1, ], control = lmerControl(optimizer ='Nelder_Mead'))
1977. ldiff_SUA <- lmer(diffrate ~ stage + (1 | Patient) + (1 | template) + (1 | unit), data_psth[data_psth$SUA == 1, ], control = lmerControl(optimizer ='Nelder_Mead'))
1979. lpos_MUA <- lmer(posrate ~ stage + (1 | Patient) + (1 | template) + (1 | unit), data_psth[data_psth$SUA == 0, ], control = lmerControl(optimizer ='Nelder_Mead'))
1980. lneg_MUA <- lmer(negrate ~ stage + (1 | Patient) + (1 | template) + (1 | unit), data_psth[data_psth$SUA == 0, ], control = lmerControl(optimizer ='Nelder_Mead'))
1981. ldiff_MUA <- lmer(diffrate ~ stage + (1 | Patient) + (1 | template) + (1 | unit), data_psth[data_psth$SUA == 0, ], control = lmerControl(optimizer ='Nelder_Mead'))
1983. lpos <- lmer(posrate ~ stage + (1 | Patient) + (1 | template) + (1 | SUA), data_psth, control = lmerControl(optimizer ='Nelder_Mead'))
1984. lneg <- lmer(negrate ~ stage + (1 | Patient) + (1 | template) + (1 | SUA), data_psth, control = lmerControl(optimizer ='Nelder_Mead'))
1985. ldiff <- lmer(diffrate ~ stage + (1 | Patient) + (1 | template) + (1 | SUA), data_psth, control = lmerControl(optimizer ='Nelder_Mead'))
1987. summary(lpos_SUA)
1988. plot_model(lpos_SUA)
1989. summary(lneg_SUA)
1990. plot_model(lneg_SUA)
1991. summary(ldiff_SUA)
1992. plot_model(ldiff_SUA)
1994. summary(lpos_MUA)
1995. plot_model(lpos_MUA)
1996. summary(lneg_MUA)
1997. plot_model(lneg_MUA)
1998. summary(ldiff_MUA)
1999. plot_model(ldiff_MUA)
2001. summary(lpos)
2002. plot_model(lpos)
2003. summary(lneg)
2004. plot_model(lneg)
2005. summary(ldiff)
2006. plot_model(ldiff)
2008. # Post-hoc tests
2009. temp = emmeans(lpos, list(pairwise ~ stage), adjust = "tukey")
2010. phpos <- as.data.frame(temp$`pairwise differences of stage`)
2011. colnames(phpos) <- c("Comparison","Estimate","SE","df","Z ratio","p")
2012. phpos$df <- NA
2014. temp = emmeans(lneg, list(pairwise ~ stage), adjust = "tukey")
2015. phneg <- as.data.frame(temp$`pairwise differences of stage`)
2016. colnames(phneg) <- c("Comparison","EstimateNeg","SENeg","dfNeg","Z ratioNeg","pNeg")
2018. temp = emmeans(lneg, list(pairwise ~ stage), adjust = "tukey")
2019. phneg <- as.data.frame(temp$`pairwise differences of stage`)
2020. colnames(phneg) <- c("Comparison","Estimate","SE","df","Z ratio","p")
2021. phneg$df <- NA
2023. temp = emmeans(ldiff, list(pairwise ~ stage), adjust = "tukey")
2024. phdiff <- as.data.frame(temp$`pairwise differences of stage`)
2025. colnames(phdiff) <- c("Comparison","Estimate","SE","df","Z ratio","p")
2026. phdiff$df <- NA
2028. #########################################
2029. # Model coefficients to table for LaTeX #
2030. #########################################
2032. # Model coefficients for table
2033. temp = summary(lpos)
2034. spos <- temp$coefficients
2035. spos <- data.frame(Predictor = row.names(spos), spos);
2036. rownames(spos) <- NULL
2037. colnames(spos) <- c("Predictor","Estimate","SD","df","t","p")
2039. temp = summary(lneg)
2040. sneg <- temp$coefficients
2041. sneg <- data.frame(Predictor = row.names(sneg), sneg);
2042. rownames(sneg) <- NULL
2043. colnames(sneg) <- c("Predictor","Estimate","SD","df","t","p")
2045. temp = summary(ldiff)
2046. sdiff <- temp$coefficients
2047. sdiff <- data.frame(Predictor = row.names(sdiff), sdiff);
2048. rownames(sdiff) <- NULL
2049. colnames(sdiff) <- c("Predictor","Estimate","SD","df","t","p")
2051. # to LaTeX table (and reorder)
2052. sneg$id <- 1:nrow(sneg)
2053. coef <- merge(sneg, spos, by="Predictor")
2054. coef <- merge(coef, sdiff, by="Predictor")
2055. coef <- coef[order(coef$id), ]
2056. coef <- coef[, c(-1, -7)]
2057. rownames(coef) <- NULL
2059. coef[,3] = round(coef[,3],digits=0)
2060. coef[,8] = round(coef[,8],digits=0)
2061. coef[,13] = round(coef[,13],digits=0)
2063. coef$Predictor <- c("\\textit{Intercept}","S3", "S2", "S1", "WASO", "Post", "REM")
2064. coef <- coef[, c(16,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15)]
2066. # posthoc coefficients
2067. phneg$id <- 1:nrow(phneg)
2069. ph <- merge(phneg,phpos, by = "Comparison")
2070. ph <- merge(ph,phdiff, by = "Comparison")
2071. ph <- ph[order(ph$id), ]
2072. ph <- ph[, -7]
2073. rownames(ph) <- NULL
2075. # Concatenate in one LaTeX table
2076. colnames(ph) <- c("Predictor",
2077. "Coef $\\beta$","SE($\\beta$)","df", "z", "\\textit{p}",
2078. "Coef $\\beta$","SE($\\beta$)","df","z", "\\textit{p}",
2079. "Coef $\\beta$","SE($\\beta$)","df","z", "\\textit{p}")
2080. colnames(coef) <- c("Predictor",
2081. "Coef $\\beta$","SE($\\beta$)","df", "z", "\\textit{p}",
2082. "Coef $\\beta$","SE($\\beta$)","df","z", "\\textit{p}",
2083. "Coef $\\beta$","SE($\\beta$)","df","z", "\\textit{p}")
2085. stats_cnt <- bind_rows(coef,ph)
2086. colnames(stats_cnt) <- c("", "Coef $\\beta$","SE($\\beta$)","df", "z",
2087. "\\textit{p}", "Coef $\\beta$","SE($\\beta$)","df","z",
2088. "\\textit{p}", "Coef $\\beta$","SE($\\beta$)","df","z",
2089. "\\textit{p}")
2091. stats_cnt[,6] = ifelse(stats_cnt[,6] > .05, paste(round(stats_cnt[,6],digits=2),sep=""),
2092. ifelse(stats_cnt[,6] < .0001, "<.0001\\textsuperscript{***}",
2093. ifelse(stats_cnt[,6] < .001,"<.001\\textsuperscript{**}",
2094. ifelse(stats_cnt[,6] < .01, "<.01\\textsuperscript{*}", "<.05\\textsuperscript{.}"))))
2095. stats_cnt[,11] = ifelse(stats_cnt[,11] > .05, paste(round(stats_cnt[,11],digits=2),sep=""),
2096. ifelse(stats_cnt[,11] < .0001, "<.0001\\textsuperscript{***}",
2097. ifelse(stats_cnt[,11] < .001,"<.001\\textsuperscript{**}",
2098. ifelse(stats_cnt[,11] < .01, "<.01\\textsuperscript{*}", "<.05\\textsuperscript{.}"))))
2099. stats_cnt[,16] = ifelse(stats_cnt[,16] > .05, paste(round(stats_cnt[,16],digits=2),sep=""),
2100. ifelse(stats_cnt[,16] < .0001, "<.0001\\textsuperscript{***}",
2101. ifelse(stats_cnt[,16] < .001,"<.001\\textsuperscript{**}",
2102. ifelse(stats_cnt[,16] < .01, "<.01\\textsuperscript{*}", "<.05\\textsuperscript{.}"))))
2104. options(knitr.kable.NA = '')
2106. kbl(stats_cnt, "latex", booktabs = T, linesep = "", label = 'unit_stats',
2107. escape = FALSE, digits = 2,
2108. caption = "Effect of sleep stage on firing rates during IED")%>%
2109. kable_styling(latex_options = c("HOLD_position"))%>%
2110. pack_rows("Sleep stages", 1, 7) %>% # latex_gap_space = "2em"
2111. pack_rows("Post-hoc comparisons", 8, 28) %>%
2112. add_header_above(c(" ", "Spike" = 5, "Wave" = 5, "Ratio" = 5)) %>%
2113. kable_styling(latex_options = c("scale_down"))%>%
2114. footnote(general_title = "",
2115. footnote_as_chunk = TRUE,
2116. threeparttable = TRUE,
2117. escape = FALSE,
2118. general = c("$.=p<0.05$, $*=p<0.01$, $**=p<0.001$, $***=p<0.0001$"))%>%
2119. save_kable("D:/Dropbox/Apps/Overleaf/Neuronal mechanisms underlying the modulation of interictal epileptic activity during sleep/tables/unit_stats.tex")
2122. #################
2123. ## Plot models ##
2124. #################
2126. set_theme(
2127. base = theme_article(),
2128. # panel.bordercol = NA
2129. )
2131. plot_cnt_pos_model <- plot_model(
2132. lpos,
2133. title = "",
2134. colors = "bw",
2135. axis.labels = "",
2136. axis.title = "",
2137. show.values = TRUE,
2138. show.p = TRUE,
2139. decimals = 4,
2140. digits = 2,
2141. value.offset = 0.5,
2142. value.size = 2.5) +
2143. # ylim(-5, 6) +
2144. font_size(labels.y = 9, labels.x = 9) +
2145. scale_x_discrete(expand=expansion(mult=c(0.1,0.2)),
2146. labels=c("Post","REM","WASO","S1","S2","S3")) +
2147. scale_y_continuous(limits = c(-6, 8), breaks = c(-6, 0, 8))
2148. plot_cnt_pos_model$data$title = "Model of spike count during LFP spike (vs. Pre-sleep)"
2149. plot_cnt_pos_model$data$title = "Model fixed effects"
2150. plot_cnt_pos_model <- plot_cnt_pos_model + facet_wrap(~title, scales="free_y")
2152. plot_cnt_neg_model <- plot_model(
2153. lneg,
2154. title = "",
2155. colors = "bw",
2156. axis.labels = "",
2157. axis.title = "",
2158. show.values = TRUE,
2159. show.p = TRUE,
2160. decimals = 4,
2161. digits = 2,
2162. value.offset = 0.5,
2163. value.size = 2.5) +
2164. ylim(-3, 2) +
2165. font_size(labels.y = 9, labels.x = 9) +
2166. scale_x_discrete(expand=expansion(mult=c(0.1,0.2)),
2167. labels=c("Post","REM","WASO","S1","S2","S3")) +
2168. scale_y_continuous(limits = c(-3, 2), breaks = c(-3, 0, 2))
2169. plot_cnt_neg_model$data$title = "Model of spike count during slow-wave (vs. Pre-sleep)"
2170. plot_cnt_neg_model$data$title = "Model fixed effects"
2171. plot_cnt_neg_model <- plot_cnt_neg_model + facet_wrap(~title, scales="free_y")
2173. plot_cnt_diff_model <- plot_model(
2174. ldiff,
2175. title = "",
2176. colors = "bw",
2177. axis.labels = "",
2178. axis.title = "",
2179. show.values = TRUE,
2180. show.p = TRUE,
2181. decimals = 4,
2182. digits = 2,
2183. value.offset = 0.5,
2184. value.size = 2.5) +
2185. ylim(-5, 7) +
2186. font_size(labels.y = 9, labels.x = 9) +
2187. scale_x_discrete(expand=expansion(mult=c(0.1,0.2)),
2188. labels=c("Post","REM","WASO","S1","S2","S3")) +
2189. scale_y_continuous(limits = c(-5, 7), breaks = c(-5, 0, 7))
2190. plot_cnt_diff_model$data$title = "Model of spike count difference between LFP spike & slow wave (vs. Pre-sleep)"
2191. plot_cnt_diff_model$data$title = "Model fixed effects"
2192. plot_cnt_diff_model <- plot_cnt_diff_model + facet_wrap(~title, scales="free_y")
2195. ggarrange(plot_cnt_pos, plot_cnt_pos_model, plot_cnt_neg, plot_cnt_neg_model, plot_cnt_diff, plot_cnt_diff_model,
2196. ncol = 2, nrow = 3,
2197. vjust = 1.5, hjust = -1,
2198. labels = c("A","B","C","D","E","F","G","H"),
2199. legend = "right",
2200. common.legend = TRUE,
2201. font.label = list(size = 14, color = "black", face = "bold")) %>%
2202. ggexport(filename = "D:/Dropbox/Apps/Overleaf/Neuronal mechanisms underlying the modulation of interictal epileptic activity during sleep/images/unit_timelocked_boxplots.pdf")
2204. #########################################
2205. # Plot correlation between LFP and PSTH #
2206. #########################################
2208. # prepare data
2209. data_rho <- read.csv(file="//lexport/iss02.charpier/analyses/stephen.whitmarsh/data/hspike/rho_table.txt", sep=',', header=TRUE, dec='.', na.strings = " ")
2210. data_rho$Patient <- factor(data_rho$Patient, levels = c(8:1))
2211. data_rho$part <- factor(data_rho$part)
2212. data_rho$template <- factor(data_rho$template)
2213. data_rho$unit <- factor(data_rho$unit)
2215. # select responsive units
2216. # data_rho <- data_rho[c(data_rho$responsive == 1), ]
2218. plot_rho <- ggplot(data=data_rho, aes(y=Patient, x=corr_rho)) +
2219. scale_fill_brewer(palette = "Set2") + scale_color_brewer(palette = "Set2") +
2220. geom_boxplot(outlier.shape = NA) +
2221. geom_point(aes(group = interaction(Patient, unit), col = Patient), position=position_dodge(width=0.3), alpha = 0.5) +
2222. # guides(colour = "none") +
2223. ylab(NULL) + xlab(NULL) +
2224. theme_article() +
2225. theme(legend.position="bottom") +
2226. coord_cartesian(xlim = c(-0.6, 0.6))
2228. final_plot = ggarrange(plot_rho,
2229. nrow = 1, ncol = 1,
2230. vjust = 1.5, hjust = -1,
2231. legend = "bottom",
2232. font.label = list(size = 14, color = "black", face = "bold")) %>%
2233. ggexport(filename = "D:/Dropbox/Apps/Overleaf/Hspike/images/rho_boxplots.pdf")
2237. ##########
2238. # Combine:
2240. data_psth
2241. data_amp
2243. ###########################
2244. # Unit baseline behaviour #
2245. ###########################
2247. data_window <- read.csv(file="Z:/analyses/stephen.whitmarsh/data/hspike/window_spike_table.txt", sep=',', header=TRUE, dec='.', na.strings = " ")
2248. data_window$hyplabel[data_window$hyplabel == "Wake"] = "WASO"
2249. data_window$hyplabel <- factor(data_window$hyplabel, levels = c("Pre", "Post", "REM", "WASO", "S1", "S2", "S3"))
2250. data_window$Patient <- factor(data_window$patient, levels = c(8:1))
2251. data_window$Night <- factor(data_window$part)
2252. data_window$Unit <- factor(data_window$unit)
2253. data_window$responsive <- factor(data_window$responsive)
2254. data_window$Type <- factor(data_window$Type)
2255. levels(data_window$Type)[levels(data_window$Type)=="good"] <- "SUA"
2256. levels(data_window$Type)[levels(data_window$Type)=="mua"] <- "MUA"
2257. data_window$Type <- factor(ordered(data_window$Type, levels = c("MUA", "SUA")))
2259. data_sel <- na.omit(data_window) # removing rows with missing data
2260. data_sel <- data_sel[data_sel$BAD_cnt == 0, ] # removing rows with artefacts
2261. data_sel <- data_sel[data_sel$IEDsum == 0, ] # removing windows with IEDs
2263. data_CV2 <- setNames(aggregate(data_sel$CV2_trial, by = c(list(data_sel$Patient, data_sel$Unit, data_sel$Type, data_sel$hyplabel)), mean), c("Patient", "unit", "Type", "hyplabel", "CV2"))
2264. data_CV2_burst <- setNames(aggregate(data_sel$CV2_intraburst_trial,
2265. by = c(list(data_sel$Patient, data_sel$Unit, data_sel$Type, data_sel$hyplabel)), mean), c("Patient", "unit", "Type", "hyplabel", "CV2_burst"))
2266. data_FR <- setNames(aggregate(data_sel$trialfreq, by = c(list(data_sel$Patient, data_sel$Unit, data_sel$Type, data_sel$hyplabel)), mean), c("Patient", "unit", "Type", "hyplabel", "FR"))
2267. data_FRcor <- setNames(aggregate(data_sel$trialfreq_corrected, by = c(list(data_sel$Patient, data_sel$Unit, data_sel$Type, data_sel$hyplabel)), mean), c("Patient", "unit", "Type", "hyplabel", "FRcor"))
2268. data_burst <- setNames(aggregate(data_sel$burst_trialsum, by = c(list(data_sel$Patient, data_sel$Unit, data_sel$Type, data_sel$hyplabel)), mean), c("Patient", "unit", "Type", "hyplabel", "Bursts"))
2269. data_amp <- setNames(aggregate(data_sel$amplitude, by = c(list(data_sel$Patient, data_sel$Unit, data_sel$Type, data_sel$hyplabel)), mean), c("Patient", "unit", "Type", "hyplabel", "Amplitude"))
2270. data_rpv <- setNames(aggregate(data_sel$RPV, by = c(list(data_sel$Patient, data_sel$Unit, data_sel$Type, data_sel$hyplabel)), mean), c("Patient", "unit", "Type", "hyplabel", "RPV"))
2271. data_mean <- merge(data_CV2, data_FR)
2272. data_mean <- merge(data_mean, data_FRcor)
2273. data_mean <- merge(data_mean, data_burst)
2274. data_mean <- merge(data_mean, data_amp)
2275. data_mean <- merge(data_mean, data_rpv)
2276. data_mean <- merge(data_mean, data_CV2_burst)
2278. # relative to pre-sleep
2279. temp <- data_CV2
2280. names(temp)[names(temp) == "CV2"] <- "Pre_CV2"
2281. temp <- temp[temp$hyplabel=="Pre", ]
2282. temp <- subset(temp, select = -c(hyplabel))
2283. data_mean <- merge(data_mean, temp)
2284. data_mean$CV2_rel <- (data_mean$CV2-data_mean$Pre_CV2) / (data_mean$CV2+data_mean$Pre_CV2)
2286. temp <- data_CV2_burst
2287. names(temp)[names(temp) == "CV2_burst"] <- "Pre_CV2_burst"
2288. temp <- temp[temp$hyplabel=="Pre", ]
2289. temp <- subset(temp, select = -c(hyplabel))
2290. data_mean <- merge(data_mean, temp)
2291. data_mean$CV2_burst_rel <- (data_mean$CV2_burst-data_mean$Pre_CV2_burst) / (data_mean$CV2_burst+data_mean$Pre_CV2_burst)
2293. temp <- data_FR
2294. names(temp)[names(temp) == "FR"] <- "Pre_FR"
2295. temp <- temp[temp$hyplabel=="Pre", ]
2296. temp <- subset(temp, select = -c(hyplabel))
2297. data_mean <- merge(data_mean, temp)
2298. data_mean$FR_rel <- (data_mean$FR-data_mean$Pre_FR) / (data_mean$FR+data_mean$Pre_FR)
2300. temp <- data_FRcor
2301. names(temp)[names(temp) == "FRcor"] <- "Pre_FRcor"
2302. temp <- temp[temp$hyplabel=="Pre", ]
2303. temp <- subset(temp, select = -c(hyplabel))
2304. data_mean <- merge(data_mean, temp)
2305. data_mean$FRcor_rel <- (data_mean$FRcor-data_mean$Pre_FRcor) / (data_mean$FRcor+data_mean$Pre_FRcor)
2307. temp <- data_amp
2308. names(temp)[names(temp) == "Amplitude"] <- "Pre_Amplitude"
2309. temp <- temp[temp$hyplabel=="Pre", ]
2310. temp <- subset(temp, select = -c(hyplabel))
2311. data_mean <- merge(data_mean, temp)
2312. data_mean$Amplitude_rel <- (data_mean$Amplitude-data_mean$Pre_Amplitude) / (data_mean$Amplitude+data_mean$Pre_Amplitude)
2314. temp <- data_burst
2315. names(temp)[names(temp) == "Bursts"] <- "Pre_Bursts"
2316. temp <- temp[temp$hyplabel=="Pre", ]
2317. temp <- subset(temp, select = -c(hyplabel))
2318. data_mean <- merge(data_mean, temp)
2319. data_mean$Bursts_rel <- (data_mean$Bursts-data_mean$Pre_Bursts) / (data_mean$Bursts+data_mean$Pre_Bursts)
2321. # Normalize per unit
2322. # FRmean <- setNames(aggregate(data_mean_sel$FR_corrected, by = c(list(data_mean_sel$unit)), mean), c("unit", "FRmean"))
2323. # FRsd <- setNames(aggregate(data_mean_sel$FR_corrected, by = c(list(data_mean_sel$unit)), sd), c("unit", "FRsd"))
2324. # data_mean_sel <- merge(data_mean_sel,FRmean)
2325. # data_mean_sel <- merge(data_mean_sel,FRsd)
2326. # data_mean_sel$FR_standardized <- (data_mean_sel$FR_corrected-data_mean_sel$FRmean)/data_mean_sel$FRsd
2327. #
2328. # Ampmean <- setNames(aggregate(data_mean_sel$Amplitude, by = c(list(data_mean_sel$unit)), mean), c("unit", "Ampmean"))
2329. # Ampsd <- setNames(aggregate(data_mean_sel$Amplitude, by = c(list(data_mean_sel$unit)), sd), c("unit", "Ampsd"))
2330. # data_mean_sel <- merge(data_mean_sel,Ampmean)
2331. # data_mean_sel <- merge(data_mean_sel,Ampsd)
2332. # data_mean_sel$Amp_standardized <- (data_mean_sel$Amplitude-data_mean_sel$Ampmean)/data_mean_sel$Ampsd
2335. # plot
2336. data_mean_sel <- data_mean[!data_mean$hyplabel == "Pre", ] # Remove for plotting
2337. data_mean_sel$title_CV2 = "CV2 (inter-spike intervals)"
2338. data_mean_sel$title_CV2_burst = "CV2 (inter-burst intervals)"
2339. data_mean_sel$title_FR = "Firing rate"
2340. data_mean_sel$title_FRcor = "Firing rate corrected"
2341. data_mean_sel$title_SUA_FRcor = "Standardized firing rate (SUA)"
2342. data_mean_sel$title_MUA_FRcor = "Standardized firing rate (MUA)"
2343. data_mean_sel$title_Amp = "Standardized amplitude"
2344. data_mean_sel$title_Burst = "Standardized burst rate"
2346. FR_plot <-
2347. ggplot(data=data_mean_sel, aes(y=hyplabel, x=FR_rel )) +
2348. scale_fill_brewer(palette = "Set2") + scale_color_brewer(palette = "Set2") +
2349. geom_boxplot(outlier.shape = NA) +
2350. geom_point(size = 0.8,
2351. aes(group = interaction(Patient, unit), col = Patient, shape = Type, fill = Patient),
2352. position = position_dodge(width = 0.5)) +
2353. scale_shape_manual(name = "Unit type", labels = c("MUA", "SUA"), values = c(15, 17)) +
2354. scale_alpha(guide = 'none') +
2355. ylab(NULL) + xlab(NULL) +
2356. theme_article() +
2357. coord_cartesian(xlim = c(-1, 1)) +
2358. scale_x_continuous(breaks = c(-1, 0, 1)) +
2359. scale_y_discrete(expand=expansion(mult=c(0.1,0.2))) +
2360. facet_wrap(~title_FR)
2362. FRcor_SUA_plot <-
2363. ggplot(data=data_mean_sel[data_mean_sel$Type=="SUA", ], aes(y=hyplabel, x=FRcor_rel )) +
2364. scale_fill_brewer(palette = "Set2") + scale_color_brewer(palette = "Set2") +
2365. geom_boxplot(outlier.shape = NA) +
2366. geom_point(size = 0.8,
2367. aes(group = interaction(Patient, unit), col = Patient, fill = Patient),
2368. position = position_dodge(width = 0.5)) +
2369. ylab(NULL) + xlab(NULL) +
2370. theme_article() +
2371. coord_cartesian(xlim = c(-1, 1)) +
2372. scale_x_continuous(breaks = c(-1, 0, 1)) +
2373. scale_y_discrete(expand=expansion(mult=c(0.1,0.2))) +
2374. facet_wrap(~title_SUA_FRcor)
2376. FRcor_MUA_plot <-
2377. ggplot(data=data_mean_sel[data_mean_sel$Type=="MUA", ], aes(y=hyplabel, x=FRcor_rel )) +
2378. scale_fill_brewer(palette = "Set2") + scale_color_brewer(palette = "Set2") +
2379. geom_boxplot(outlier.shape = NA) +
2380. geom_point(size = 0.8,
2381. aes(group = interaction(Patient, unit), col = Patient, fill = Patient),
2382. position = position_dodge(width = 0.5)) +
2383. ylab(NULL) + xlab(NULL) +
2384. theme_article() +
2385. coord_cartesian(xlim = c(-1, 1)) +
2386. scale_x_continuous(breaks = c(-1, 0, 1)) +
2387. scale_y_discrete(expand=expansion(mult=c(0.1,0.2))) +
2388. facet_wrap(~title_MUA_FRcor)
2390. Amp_plot <-
2391. ggplot(data=data_mean_sel[data_mean_sel$Type=="SUA", ], aes(y=hyplabel, x=Amplitude_rel)) +
2392. scale_fill_brewer(palette = "Set2") + scale_color_brewer(palette = "Set2") +
2393. geom_boxplot(outlier.shape = NA) +
2394. geom_point(size = 0.8,
2395. aes(group = interaction(Patient, unit), col = Patient, fill = Patient),
2396. position = position_dodge(width = 0.5)) +
2397. ylab(NULL) + xlab(NULL) +
2398. theme_article() +
2399. coord_cartesian(xlim = c(-0.2, 0.2)) +
2400. scale_x_continuous(breaks = c(-0.2, 0, 0.2)) +
2401. scale_y_discrete(expand=expansion(mult=c(0.1,0.2))) +
2402. facet_wrap(~title_Amp)
2404. CV2_plot <- # original data, not relative, since CV2 is already standardized
2405. ggplot(data=data_mean_sel[data_mean_sel$Type == "SUA", ], aes(y=hyplabel, x=CV2)) +
2406. scale_fill_brewer(palette = "Set2") + scale_color_brewer(palette = "Set2") +
2407. geom_boxplot(outlier.shape = NA) +
2408. geom_point(size = 0.8, # shape = 17,
2409. aes(group = interaction(Patient, unit), col = Patient, fill = Patient),
2410. position = position_dodge(width = 0.5)) +
2411. ylab(NULL) + xlab(NULL) +
2412. theme_article() +
2413. # coord_cartesian(xlim = c(-0.1, 0.1)) +
2414. coord_cartesian(xlim = c(0.6, 1.3)) +
2415. # scale_x_continuous(breaks = c(-0.1, 0, 0.1)) +
2416. scale_x_continuous(breaks = c(0.6, 1, 1.3)) +
2417. scale_y_discrete(expand=expansion(mult=c(0.1,0.2))) +
2418. facet_wrap(~title_CV2)
2420. CV2_burst_plot <- # original data, not relative, since CV2 is already standardized
2421. ggplot(data=data_mean_sel[data_mean_sel$Type == "SUA", ], aes(y=hyplabel, x=CV2_burst)) +
2422. scale_fill_brewer(palette = "Set2") + scale_color_brewer(palette = "Set2") +
2423. geom_boxplot(outlier.shape = NA) +
2424. geom_point(size = 0.8, # shape = 17,
2425. aes(group = interaction(Patient, unit), col = Patient, fill = Patient),
2426. position = position_dodge(width = 0.5)) +
2427. ylab(NULL) + xlab(NULL) +
2428. theme_article() +
2429. # coord_cartesian(xlim = c(-0.6, 0.6)) +
2430. # scale_x_continuous(breaks = c(-0.6, 0, 0.6)) +
2431. coord_cartesian(xlim = c(0, 1)) +
2432. scale_x_continuous(breaks = c(0, 0.5, 1)) +
2433. scale_y_discrete(expand=expansion(mult=c(0.1,0.2))) +
2434. facet_wrap(~title_CV2_burst)
2436. Burst_plot <- ggplot(data=data_mean_sel[data_mean_sel$Type == "SUA", ], aes(y=hyplabel, x=Bursts_rel)) +
2437. scale_fill_brewer(palette = "Set2") + scale_color_brewer(palette = "Set2") +
2438. geom_boxplot(outlier.shape = NA) +
2439. geom_point(size = 0.8, # shape = 17,
2440. aes(group = interaction(Patient, unit), col = Patient, fill = Patient),
2441. position = position_dodge(width = 0.5)) +
2442. # scale_shape_manual(name = "Unit type", labels = c("MUA", "SUA"), values = c(15, 17)) +
2443. # scale_shape_discrete(label = c("MUA", "SUA")) +
2444. ylab(NULL) + xlab(NULL) +
2445. theme_article() +
2446. coord_cartesian(xlim = c(-0.8, 0.8)) +
2447. scale_x_continuous(breaks = c(-0.8, 0, 0.8)) +
2448. scale_y_discrete(expand=expansion(mult=c(0.1,0.2))) +
2449. facet_wrap(~title_Burst)
2452. ################
2453. ## STATISTICS ##
2454. ################
2456. # to create p-values
2457. detach(package:lmerTest)
2458. library(lmerTest)
2459. library(lme4)
2461. # load and prepare data
2462. data_window <- read.csv(file="Z:/analyses/stephen.whitmarsh/data/hspike/window_spike_table.txt", sep=',', header=TRUE, dec='.', na.strings = " ")
2463. data_window$hyplabel[data_window$hyplabel == "Wake"] = "WASO"
2464. data_window$hyplabel <- factor(data_window$hyplabel, levels = c("Pre", "Post", "REM", "WASO", "S1", "S2", "S3"))
2465. data_window$Patient <- factor(data_window$patient, levels = c(8:1))
2466. data_window$Night <- factor(data_window$part)
2467. data_window$Unit <- factor(data_window$unit + data_window$part * 100) # add becaue template nr repeats over nights
2468. data_window$responsive <- factor(data_window$responsive)
2469. data_window$Type <- factor(data_window$Type)
2470. data_window$FR <- data_window$trialfreq
2471. data_window$FRcor <- data_window$trialfreq_corrected
2472. levels(data_window$Type)[levels(data_window$Type)=="good"] <- "SUA"
2473. levels(data_window$Type)[levels(data_window$Type)=="mua"] <- "MUA"
2474. data_window$Type <- ordered(data_window$Type, levels = c("MUA", "SUA"))
2476. data_sel <- na.omit(data_window) # remove rows with missing data
2477. data_sel <- data_sel[data_sel$BAD_cnt == 0, ] # remove rows with artefacts
2478. data_sel <- data_sel[data_sel$IEDsum == 0, ] # remove windows with IEDs
2480. # Reorder for table
2481. data_sel$stage <- factor(data_sel$hyplabel, levels = c("Pre","S3", "S2", "S1", "WASO", "REM", "Post"))
2483. # fit model
2484. data_sel$stage = relevel(data_sel$stage, ref="Pre")
2486. lFR <- lmer(FR ~ stage + (1 | Patient) + (1 | Unit), data_sel[data_sel$Type=="MUA", ])
2487. lFRcor_SUA <- lmer(FRcor ~ stage + (1 | Patient) + (1 | Unit), data_sel[data_sel$Type=="SUA", ])
2488. lFRcor_MUA <- lmer(FRcor ~ stage + (1 | Patient) + (1 | Unit), data_sel[data_sel$Type=="MUA", ])
2489. lAmp <- lmer(amplitude ~ stage + (1 | Patient) + (1 | Unit), data_sel)
2490. lCV2 <- lmer(CV2_trial ~ stage + (1 | Patient) + (1 | Unit), data_sel)
2491. lCV2_burst <- lmer(CV2_intraburst_trial ~ stage + (1 | Patient) + (1 | Unit), data_sel)
2492. lBurst <- lmer(burst_trialsum ~ stage + (1 | Patient) + (1 | Unit), data_sel)
2494. summary(lFR)
2495. summary(lFRcor_SUA)
2496. summary(lFRcor_MUA)
2497. summary(lAmp)
2498. summary(lCV2)
2499. summary(lCV2_burst)
2500. summary(lBurst)
2502. plot_model(lFR)
2503. plot_model(lFRcor_SUA)
2504. plot_model(lFRcor_MUA)
2505. plot_model(lAmp)
2506. plot_model(lCV2)
2507. plot_model(lCV2_burst)
2508. plot_model(lBurst)
2510. # Post-hoc tests
2511. temp = emmeans(lFR, list(pairwise ~ stage), adjust = "tukey")
2512. phFR <- as.data.frame(temp$`pairwise differences of stage`)
2513. colnames(phFR) <- c("Comparison","Estimate","SE","df","Z ratio","p")
2514. phFR$df <- NA
2516. temp = emmeans(lFRcor_SUA, list(pairwise ~ stage), adjust = "tukey")
2517. phFR_SUA <- as.data.frame(temp$`pairwise differences of stage`)
2518. colnames(phFR_SUA) <- c("Comparison","Estimate","SE","df","Z ratio","p")
2519. phFR_SUA$df <- NA
2521. temp = emmeans(lFRcor_MUA, list(pairwise ~ stage), adjust = "tukey")
2522. phFR_MUA <- as.data.frame(temp$`pairwise differences of stage`)
2523. colnames(phFR_MUA) <- c("Comparison","Estimate","SE","df","Z ratio","p")
2524. phFR_MUA$df <- NA
2526. temp = emmeans(lFRcor, list(pairwise ~ stage), adjust = "tukey")
2527. phFRcor <- as.data.frame(temp$`pairwise differences of stage`)
2528. colnames(phFRcor) <- c("Comparison","Estimate","SE","df","Z ratio","p")
2529. phFRcor$df <- NA
2531. temp = emmeans(lAmp, list(pairwise ~ stage), adjust = "tukey")
2532. phAmp <- as.data.frame(temp$`pairwise differences of stage`)
2533. colnames(phAmp) <- c("Comparison","Estimate","SE","df","Z ratio","p")
2534. phAmp$df <- NA
2536. temp = emmeans(lCV2, list(pairwise ~ stage), adjust = "tukey")
2537. phCV2 <- as.data.frame(temp$`pairwise differences of stage`)
2538. colnames(phCV2) <- c("Comparison","Estimate","SE","df","Z ratio","p")
2539. phCV2$df <- NA
2541. temp = emmeans(lCV2_burst, list(pairwise ~ stage), adjust = "tukey")
2542. phCV2_burst <- as.data.frame(temp$`pairwise differences of stage`)
2543. colnames(phCV2_burst) <- c("Comparison","Estimate","SE","df","Z ratio","p")
2544. phCV2_burst$df <- NA
2546. temp = emmeans(lBurst, list(pairwise ~ stage), adjust = "tukey")
2547. phBurst <- as.data.frame(temp$`pairwise differences of stage`)
2548. colnames(phBurst) <- c("Comparison","Estimate","SE","df","Z ratio","p")
2549. phBurst$df <- NA
2551. #########################################
2552. # Model coefficients to table for LaTeX #
2553. #########################################
2555. # Model coefficients for table
2556. temp = summary(lFR)
2557. sFR <- temp$coefficients
2558. sFR <- data.frame(Predictor = row.names(sFR), sFR);
2559. rownames(sFR) <- NULL
2560. colnames(sFR) <- c("Predictor","Estimate","SD","df","t","p")
2562. temp = summary(lFRcor_SUA)
2563. sFRcor_SUA <- temp$coefficients
2564. sFRcor_SUA <- data.frame(Predictor = row.names(sFRcor_SUA), sFRcor_SUA);
2565. rownames(sFRcor_SUA) <- NULL
2566. colnames(sFRcor_SUA) <- c("Predictor","Estimate","SD","df","t","p")
2568. temp = summary(lFRcor_MUA)
2569. sFRcor_MUA <- temp$coefficients
2570. sFRcor_MUA <- data.frame(Predictor = row.names(sFRcor_MUA), sFRcor_MUA);
2571. rownames(sFRcor_MUA) <- NULL
2572. colnames(sFRcor_MUA) <- c("Predictor","Estimate","SD","df","t","p")
2574. temp = summary(lAmp)
2575. sAmp <- temp$coefficients
2576. sAmp <- data.frame(Predictor = row.names(sAmp), sAmp);
2577. rownames(sAmp) <- NULL
2578. colnames(sAmp) <- c("Predictor","Estimate","SD","df","t","p")
2580. temp = summary(lCV2)
2581. sCV2 <- temp$coefficients
2582. sCV2 <- data.frame(Predictor = row.names(sCV2), sCV2);
2583. rownames(sCV2) <- NULL
2584. colnames(sCV2) <- c("Predictor","Estimate","SD","df","t","p")
2586. temp = summary(lCV2_burst)
2587. sCV2_burst <- temp$coefficients
2588. sCV2_burst <- data.frame(Predictor = row.names(sCV2_burst), sCV2_burst);
2589. rownames(sCV2_burst) <- NULL
2590. colnames(sCV2_burst) <- c("Predictor","Estimate","SD","df","t","p")
2592. temp = summary(lBurst)
2593. sBurst <- temp$coefficients
2594. sBurst <- data.frame(Predictor = row.names(sBurst), sBurst);
2595. rownames(sBurst) <- NULL
2596. colnames(sBurst) <- c("Predictor","Estimate","SD","df","t","p")
2598. # Model coefficients to LaTeX table (and reorder)
2599. sFRcor_SUA$id <- 1:nrow(sFRcor_SUA)
2600. coef <- merge(sFRcor_SUA, sFRcor_MUA, by="Predictor")
2601. coef <- merge(coef, sBurst, by="Predictor")
2602. coef <- merge(coef, sAmp, by="Predictor")
2603. coef <- merge(coef, sCV2, by="Predictor")
2604. coef <- merge(coef, sCV2_burst, by="Predictor")
2606. coef <- coef[order(coef$id), ]
2607. coef <- coef[, c(-1, -7)]
2608. rownames(coef) <- NULL
2610. coef[,3] = round(coef[,3],digits=0)
2611. coef[,8] = round(coef[,8],digits=0)
2612. coef[,13] = round(coef[,13],digits=0)
2613. coef[,18] = round(coef[,18],digits=0)
2614. coef[,23] = round(coef[,23],digits=0)
2616. coef$Predictor <- c("\\textit{Intercept}","S3", "S2", "S1", "WASO", "REM", "Post")
2617. coef <- coef[, c(31,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30)]
2619. # posthoc coefficients
2620. phFR_SUA$id <- 1:nrow(phFR_SUA)
2622. ph <- merge(phFR_SUA, phFR_MUA, by="Comparison")
2623. ph <- merge(ph,phBurst, by = "Comparison")
2624. ph <- merge(ph,phAmp, by = "Comparison")
2625. ph <- merge(ph,phCV2, by = "Comparison")
2626. ph <- merge(ph,phCV2_burst, by = "Comparison")
2627. ph <- ph[order(ph$id), ]
2628. ph <- ph[, -7]
2629. rownames(ph) <- NULL
2631. # Concatenate in one LaTeX table
2632. colnames(ph) <- c("Predictor", "Coef $\\beta$","SE($\\beta$)","df", "z",
2633. "\\textit{p}", "Coef $\\beta$","SE($\\beta$)","df","z",
2634. "\\textit{p}", "Coef $\\beta$","SE($\\beta$)","df","z",
2635. "\\textit{p}", "Coef $\\beta$","SE($\\beta$)","df","z",
2636. "\\textit{p}", "Coef $\\beta$","SE($\\beta$)","df","z",
2637. "\\textit{p}", "Coef $\\beta$","SE($\\beta$)","df","z",
2638. "\\textit{p}")
2639. colnames(coef) <- c("Predictor", "Coef $\\beta$","SE($\\beta$)","df", "z",
2640. "\\textit{p}", "Coef $\\beta$","SE($\\beta$)","df","z",
2641. "\\textit{p}", "Coef $\\beta$","SE($\\beta$)","df","z",
2642. "\\textit{p}", "Coef $\\beta$","SE($\\beta$)","df","z",
2643. "\\textit{p}", "Coef $\\beta$","SE($\\beta$)","df","z",
2644. "\\textit{p}", "Coef $\\beta$","SE($\\beta$)","df","z",
2645. "\\textit{p}")
2647. stats_window <- bind_rows(coef,ph)
2648. colnames(stats_window) <- c("", "Coef $\\beta$","SE($\\beta$)","df", "z",
2649. "\\textit{p}", "Coef $\\beta$","SE($\\beta$)","df","z",
2650. "\\textit{p}", "Coef $\\beta$","SE($\\beta$)","df","z",
2651. "\\textit{p}", "Coef $\\beta$","SE($\\beta$)","df","z",
2652. "\\textit{p}", "Coef $\\beta$","SE($\\beta$)","df","z",
2653. "\\textit{p}", "Coef $\\beta$","SE($\\beta$)","df","z",
2654. "\\textit{p}")
2655. stats_window[,6] = ifelse(stats_window[,6] > .05, paste(round(stats_window[,6], digits=2),sep=""),
2656. ifelse(stats_window[,6] < .0001, "<.0001\\textsuperscript{***}",
2657. ifelse(stats_window[,6] < .001, "<.001\\textsuperscript{**}",
2658. ifelse(stats_window[,6] < .01, "<.01\\textsuperscript{*}", "<.05\\textsuperscript{.}"))))
2659. stats_window[,11] = ifelse(stats_window[,11] > .05, paste(round(stats_window[,11],digits=2),sep=""),
2660. ifelse(stats_window[,11] < .0001, "<.0001\\textsuperscript{***}",
2661. ifelse(stats_window[,11] < .001,"<.001\\textsuperscript{**}",
2662. ifelse(stats_window[,11] < .01, "<.01\\textsuperscript{*}", "<.05\\textsuperscript{.}"))))
2663. stats_window[,16] = ifelse(stats_window[,16] > .05, paste(round(stats_window[,16],digits=2),sep=""),
2664. ifelse(stats_window[,16] < .0001, "<.0001\\textsuperscript{***}",
2665. ifelse(stats_window[,16] < .001,"<.001\\textsuperscript{**}",
2666. ifelse(stats_window[,16] < .01, "<.01\\textsuperscript{*}", "<.05\\textsuperscript{.}"))))
2667. stats_window[,21] = ifelse(stats_window[,21] > .05, paste(round(stats_window[,21],digits=2),sep=""),
2668. ifelse(stats_window[,21] < .0001, "<.0001\\textsuperscript{***}",
2669. ifelse(stats_window[,21] < .001,"<.001\\textsuperscript{**}",
2670. ifelse(stats_window[,21] < .01, "<.01\\textsuperscript{*}", "<.05\\textsuperscript{.}"))))
2671. stats_window[,26] = ifelse(stats_window[,26] > .05, paste(round(stats_window[,26],digits=2),sep=""),
2672. ifelse(stats_window[,26] < .0001, "<.0001\\textsuperscript{***}",
2673. ifelse(stats_window[,26] < .001,"<.001\\textsuperscript{**}",
2674. ifelse(stats_window[,26] < .01, "<.01\\textsuperscript{*}", "<.05\\textsuperscript{.}"))))
2675. stats_window[,31] = ifelse(stats_window[,31] > .05, paste(round(stats_window[,31],digits=2),sep=""),
2676. ifelse(stats_window[,31] < .0001, "<.0001\\textsuperscript{***}",
2677. ifelse(stats_window[,31] < .001,"<.001\\textsuperscript{**}",
2678. ifelse(stats_window[,31] < .01, "<.01\\textsuperscript{*}", "<.05\\textsuperscript{.}"))))
2679. options(knitr.kable.NA = '')
2681. kbl(stats_window, "latex", booktabs = T, linesep = "", label = 'unit_window_stats',
2682. escape = FALSE, digits = 2,
2683. caption = "Effect of sleep stage on spontanious neuronal behaviour")%>%
2684. kable_styling(latex_options = c("HOLD_position"))%>%
2685. pack_rows("Sleep stages", 1, 7) %>% # latex_gap_space = "2em"
2686. pack_rows("Post-hoc comparisons", 8, 28) %>%
2687. add_header_above(c(" ", "Firing Rate (SUA)" = 5, "Firing Rate (MUA)" = 5, "Bursts" = 5, "Amplitude" = 5, "CV2" = 5, "CV2 (bursts)" = 5)) %>%
2688. kable_styling(latex_options = c("scale_down"))%>%
2689. footnote(general_title = "",
2690. footnote_as_chunk = TRUE,
2691. threeparttable = TRUE,
2692. escape = FALSE,
2693. general = c("$.=p<0.05$, $*=p<0.01$, $**=p<0.001$, $***=p<0.0001$"))%>%
2694. save_kable("D:/Dropbox/Apps/Overleaf/Neuronal mechanisms underlying the modulation of interictal epileptic activity during sleep/tables/unit_window_stats.tex")
2697. set_theme(
2698. base = theme_article(),
2699. # panel.bordercol = NA
2700. axis.textsize = 0.9
2701. )
2703. FRcor_SUA_plot_model <- plot_model(
2704. lFRcor_SUA,
2705. title = "",
2706. colors = "bw",
2707. axis.title = "",
2708. show.values = TRUE,
2709. show.p = TRUE,
2710. decimals = 4,
2711. digits = 2,
2712. value.offset = 0.5,
2713. value.size = 2.5) +
2714. font_size(labels.y = 9, labels.x = 9) +
2715. scale_x_discrete(expand=expansion(mult=c(0.1,0.2)),
2716. labels=c("Post","REM","WASO","S1","S2","S3")) +
2717. scale_y_continuous(limits = c(-4.2, -1.8), breaks = c(-4, -2))
2719. FRcor_MUA_plot_model <- plot_model(
2720. lFRcor_MUA,
2721. title = "",
2722. colors = "bw",
2723. axis.title = "",
2724. show.values = TRUE,
2725. show.p = TRUE,
2726. decimals = 4,
2727. digits = 2,
2728. value.offset = 0.5,
2729. value.size = 2.5) +
2730. font_size(labels.y = 9, labels.x = 9) +
2731. scale_x_discrete(expand=expansion(mult=c(0.1,0.2)),
2732. labels=c("Post","REM","WASO","S1","S2","S3")) +
2733. scale_y_continuous(limits = c(-7, 1), breaks = c(-7, 0, 1))
2735. CV2_plot_model <- plot_model(
2736. lCV2,
2737. title = "",
2738. colors = "bw",
2739. axis.title = "",
2740. show.values = TRUE,
2741. show.p = TRUE,
2742. decimals = 4,
2743. digits = 2,
2744. value.offset = 0.5,
2745. value.size = 2.5,
2746. axis.labels = c("REM","Post","WASO","S1","S2","S3")) +
2747. font_size(labels.y = 9, labels.x = 9) +
2748. scale_x_discrete(expand=expansion(mult=c(0.1, 0.2)),
2749. labels=c("Post","REM","WASO","S1","S2","S3")) +
2750. scale_y_continuous(limits = c(0, 0.12), breaks = c(0, 0.1))
2753. CV2_burst_plot_model <- plot_model(
2754. lCV2_burst,
2755. title = "",
2756. colors = "bw",
2757. axis.title = "",
2758. show.values = TRUE,
2759. show.p = TRUE,
2760. decimals = 4,
2761. digits = 2,
2762. value.offset = 0.5,
2763. value.size = 2.5,
2764. axis.labels = c("REM","Post","WASO","S1","S2","S3")) +
2765. font_size(labels.y = 9, labels.x = 9) +
2766. ylim(-0.05, 0.05) +
2767. scale_x_discrete(expand=expansion(mult=c(0.1,0.2)),
2768. labels=c("Post","REM","WASO","S1","S2","S3")) +
2769. scale_y_continuous(limits = c(-0.04, 0.01), breaks = c(-0.04, 0, 0.01))
2771. Amp_plot_model <- plot_model(
2772. lAmp,
2773. title = "",
2774. colors = "bw",
2775. axis.title = "",
2776. show.values = TRUE,
2777. show.p = TRUE,
2778. decimals = 4,
2779. digits = 2,
2780. value.offset = 0.5,
2781. value.size = 2.5,
2782. axis.labels = c("REM","Post","WASO","S1","S2","S3")) +
2783. font_size(labels.y = 9, labels.x = 9) +
2784. scale_x_discrete(expand=expansion(mult=c(0.1,0.2)),
2785. labels=c("Post","REM","WASO","S1","S2","S3")) +
2786. scale_y_continuous(limits = c(-1, 2), breaks = c(-1, 0, 2))
2788. Burst_plot_model <- plot_model(
2789. lBurst,
2790. title = "",
2791. colors = "bw",
2792. axis.title = "",
2793. show.values = TRUE,
2794. show.p = TRUE,
2795. decimals = 4,
2796. digits = 2,
2797. value.offset = 0.5,
2798. value.size = 2.5,
2799. axis.labels = c("REM","Post","WASO","S1","S2","S3")) +
2800. font_size(labels.y = 9, labels.x = 9) +
2801. scale_x_discrete(expand=expansion(mult=c(0.1,0.2)),
2802. labels=c("Post","REM","WASO","S1","S2","S3")) +
2803. scale_y_continuous(limits = c(-20, 0), breaks = c(-20, 0))
2805. Burst_plot_model$data$title = "Model fixed effects"
2806. Burst_plot_model <- Burst_plot_model + facet_wrap(~title, scales="free_y")
2807. # FR_plot_model$data$title = "Model fixed effects"
2808. # FR_plot_model <- FR_plot_model + facet_wrap(~title, scales="free_y")
2809. FRcor_SUA_plot_model$data$title = "Model fixed effects"
2810. FRcor_SUA_plot_model <- FRcor_SUA_plot_model + facet_wrap(~title, scales="free_y")
2811. FRcor_MUA_plot_model$data$title = "Model fixed effects"
2812. FRcor_MUA_plot_model <- FRcor_MUA_plot_model + facet_wrap(~title, scales="free_y")
2813. Amp_plot_model$data$title = "Model fixed effects"
2814. Amp_plot_model <- Amp_plot_model + facet_wrap(~title, scales="free_y")
2815. CV2_plot_model$data$title = "Model fixed effects"
2816. CV2_plot_model <- CV2_plot_model + facet_wrap(~title, scales="free_y")
2817. CV2_burst_plot_model$data$title = "Model fixed effects"
2818. CV2_burst_plot_model <- CV2_burst_plot_model + facet_wrap(~title, scales="free_y")
2820. ggarrange(FRcor_SUA_plot, FRcor_SUA_plot_model, FRcor_MUA_plot, FRcor_MUA_plot_model, Burst_plot, Burst_plot_model,
2821. ncol = 2, nrow = 3,
2822. vjust = 1.5, hjust = -1,
2823. labels = c("A","B","C","D","E","F","G","H"),
2824. legend = "right",
2825. common.legend = TRUE,
2826. font.label = list(size = 14, color = "black", face = "bold")) %>%
2827. ggexport(filename = "D:/Dropbox/Apps/Overleaf/Neuronal mechanisms underlying the modulation of interictal epileptic activity during sleep/images/window_FR_boxplots.pdf")
2830. ggarrange(Amp_plot, Amp_plot_model, CV2_plot, CV2_plot_model, CV2_burst_plot, CV2_burst_plot_model,
2831. ncol = 2, nrow = 3,
2832. vjust = 1.5, hjust = -1,
2833. labels = c("A","B","C","D","E","F","G","H"),
2834. legend = "right",
2835. common.legend = TRUE,
2836. font.label = list(size = 14, color = "black", face = "bold")) %>%
2837. ggexport(filename = "D:/Dropbox/Apps/Overleaf/Neuronal mechanisms underlying the modulation of interictal epileptic activity during sleep/images/window_CV2_boxplots.pdf")
2840. ##########################################################################
2841. ### Combine spike and power data
2842. ##########################################################################
2846. # load and prepare data
2847. data_window <- read.csv(file="Z:/analyses/stephen.whitmarsh/data/hspike/window_spike_table.txt", sep=',', header=TRUE, dec='.', na.strings = " ")
2848. data_window$hyplabel <- factor(data_window$hyplabel, levels = c("NO_SCORE", "Pre", "Post", "REM", "Wake", "S1", "S2", "S3"))
2849. data_window$Patient <- factor(data_window$patient, levels = c(8:1))
2850. data_window$Night <- factor(data_window$part)
2851. data_window$Unit <- factor(data_window$unit + data_window$part * 100) # add becaue template nr repeats over nights
2852. data_window$responsive <- factor(data_window$responsive)
2853. data_window$Type <- factor(data_window$Type)
2854. data_window$FR <- data_window$trialfreq
2855. data_window$FRcor <- data_window$trialfreq_corrected
2856. levels(data_window$Type)[levels(data_window$Type)=="good"] <- "SUA"
2857. levels(data_window$Type)[levels(data_window$Type)=="mua"] <- "MUA"
2858. data_window$Type <- ordered(data_window$Type, levels = c("MUA", "SUA"))
2861. # prepare data
2862. data_power <- read.csv(file="Z:/analyses/stephen.whitmarsh/data/hspike/power_table_long.txt", sep=',', header=TRUE, dec='.', na.strings = " ")
2863. data_power$hyplabel <- factor(data_power$hyplabel, ordered = TRUE, levels = c("NO_SCORE", "REM", "AWAKE", "PHASE_1", "PHASE_2", "PHASE_3"))
2864. # data_power$band <- factor(data_power$band, ordered = TRUE, levels = c("delta", "theta", "alpha", "beta", "delta_div_alpha"))
2865. data_power$band <- factor(data_power$band, ordered = TRUE, levels = c("Delta1", "Delta2"))
2866. data_power$Patient <- factor(data_power$patient, levels = c(8:1))
2867. data_power$part <- factor(data_power$part)
2868. data_power$power_log <- log(data_power$power)
2870. # data_sel <- na.omit(data_window) # remove rows with missing data
2871. # data_sel <- data_sel[data_sel$BAD_cnt == 0, ] # remove rows with artefacts
2872. # data_sel <- data_sel[data_sel$IEDsum == 0, ] # remove windows with IEDs
2874. d <- merge(data_window, data_power, by = "starttime")
2875. d <- d[d$band == "Delta2", ]
2877. d2 <- d %>% mutate(power_bin = cut(power, breaks=seq(-1.5,10.5,1)))
2878. levels(d2$power_bin) <- seq(1:12)-2
2880. d$IEDsum.x[105:108] - d$IEDsum.y[105:108]
2881. d$IEDsum.x[665:672] - d$IEDsum.y[665:672]
2882. d$IEDsum.x
2883. d$IEDsum.y
2885. test <- d$IEDsum.x - d$IEDsum.y
2886. hist(test[test != 0])
2888. h1 <- hist(d$IEDsum.y[d$IEDsum.y != 0])
2889. h2 <- hist(d$IEDsum.x[d$IEDsum.x != 0])
2890. h3 <- hist(d$IEDsum.x - d$IEDsum.y, 10)
2892. d[104:109, ]
2894. # bin data
2896. # prepare datastructure
2897. d2 <- d2 %>% group_by(power_bin, IEDsum)
2898. d2 <- d2 %>% summarise(count = n(), Patient)
2899. d2 <- d2[order(data_pow_wide2$Patient), ]
2901. # ggplot(data=data_pow_wide2[data_pow_wide2$IEDsum > 0, ], aes(x=Delta1_log_bin, y=IEDsum*6)) +
2902. D2 <- ggplot(data=data_pow_wide2, aes(x=Delta_log_bin, y=IEDsum*6, color = Patient)) +
2903. # geom_count(aes(size = after_stat(prop), group = Delta_log_bin, color = Patient)) +
2904. geom_count() +
2905. # scale_size_area(max_size = 3) +
2906. # geom_smooth(data=data_pow_wide, aes(x=Delta1_log, y=IEDsum), method="lm", fullrange = FALSE, color="red") +
2907. scale_fill_brewer(palette = "Set2") + scale_color_brewer(palette = "Set2") +
2908. theme(axis.text.y = element_blank()) +
2909. theme_article() +
2910. theme(
2911. strip.background = element_blank(),
2912. strip.text.x = element_blank()
2913. ) +
2914. # xlab("Delta2 power (log)")
2915. #xlab = expression(mu "Volts" / "Hertz" ^ 2) +
2916. # xlab = expression("Force spaces with ~" ~ mu ~ pi * sigma ~ pi) +
2917. # ylab( units~are~(mu*g)/L )
2918. # xlab(TeX(r'($\alpha x^\alpha$, where $\alpha \in \{1 \ldots 5\}$)')) +
2919. # xlab(TeX(r'($log(\mu V/Hz^2$))')) +
2920. xlab(TeX(r'($log(\mu V^2/Hz$))')) +
2921. ylab("IEDs per minute") +
2922. labs(size = "Observations") +
2923. # coord_cartesian(xlim = c(2,11)) +
2924. # coord_cartesian(ylim = c(2,11)) +
2925. # facet_wrap(~Patient, ncol = 4, scales="free_x")
2926. facet_wrap(~Patient, ncol = 4)
2927. # ggsave("D:/Dropbox/Apps/Overleaf/Hspike/images/IEDrate_Delta2_count.pdf", width = 7, height = 4, units = "in")
2929. ggarrange(D1, D2,
2930. ncol = 1, nrow = 2,
2931. vjust = 1, hjust = 0,
2932. labels = c("A","B"),
2933. legend = "right",
2934. common.legend = TRUE,
2935. font.label = list(size = 14, color = "black", face = "bold")) %>%
2936. ggexport(filename = "D:/Dropbox/Apps/Overleaf/Hspike/images/IEDrate_Delta_count.pdf")
2952. #
2953. # ## POLAR ##
2954. #
2955. # library("ggplot2")
2956. # library("dplyr")
2957. #
2958. # ## plot power values PER HOUR
2959. #
2960. # data$bin <- as.integer(cut(data$minute, seq(0, 24*60, by = 60)))
2961. # data_binned <- setNames(aggregate(data$alpha, c(list(data$patient), list(data$bin)), mean), c("patient", "bin", "alpha"))
2962. # data_binned <- as.data.frame(data_binned %>% group_by(patient) %>% mutate(Nalpha = alpha/max(alpha)))
2963. #
2964. # ggplot(data=data_binned, aes(x = bin, y = alpha, fill=patient, col=patient)) +
2965. # scale_fill_brewer(palette = "Set2") + scale_color_brewer(palette = "Set2") +
2966. # geom_bar(stat="identity", width = 1) +
2967. # geom_vline(xintercept = seq(0.5, 21.5, by = 3), colour = "grey90") +
2968. # theme_bw() +
2969. # theme(panel.border = element_blank(),
2970. # legend.key = element_blank(),
2971. # axis.ticks = element_blank(),
2972. # axis.text.y = element_blank(),
2973. # panel.grid = element_blank(),
2974. # ) +
2975. # scale_x_continuous(breaks=seq(0.5, 21.5, by = 3), labels = c("0" = "00:00", "3" = "03:00", "6" = "06:00", "9" = "09:00", "12" = "12:00", "15" = "15:00", "18" = "18:00", "21" = "21:00")) +
2976. # coord_polar(theta = "x", start = 0)
2977. #
2978. #
2979. # ## plot power values (use bindivision for hour divisions)
2980. # bindivision = 1
2981. # data$bin <- as.integer(cut(data$minute, seq(0, 24*60, by = 60/bindivision)))
2982. # y1 <- setNames(aggregate(data$alpha, c(list(data$patient), list(data$bin)), median), c("patient", "bin", "alpha"))
2983. # y2 <- setNames(aggregate(data$theta, c(list(data$patient), list(data$bin)), median), c("patient", "bin", "theta"))
2984. # y3 <- setNames(aggregate(data$beta, c(list(data$patient), list(data$bin)), median), c("patient", "bin", "beta"))
2985. # y4 <- setNames(aggregate(data$delta, c(list(data$patient), list(data$bin)), median), c("patient", "bin", "delta"))
2986. # y1sd <- setNames(aggregate(data$alpha, c(list(data$patient), list(data$bin)), sd), c("patient", "bin", "SDalpha"))
2987. # y2sd <- setNames(aggregate(data$theta, c(list(data$patient), list(data$bin)), sd), c("patient", "bin", "SDtheta"))
2988. # y3sd <- setNames(aggregate(data$beta, c(list(data$patient), list(data$bin)), sd), c("patient", "bin", "SDbeta"))
2989. # y4sd <- setNames(aggregate(data$delta, c(list(data$patient), list(data$bin)), sd), c("patient", "bin", "SDdelta"))
2990. # data_binned <- merge(y1,y2)
2991. # data_binned <- merge(data_binned,y3)
2992. # data_binned <- merge(data_binned,y4)
2993. # data_binned <- merge(data_binned,y1sd)
2994. # data_binned <- merge(data_binned,y2sd)
2995. # data_binned <- merge(data_binned,y3sd)
2996. # data_binned <- merge(data_binned,y4sd)
2997. # data_binned <- as.data.frame(data_binned %>% group_by(patient) %>% mutate(Nalpha = (alpha-min(alpha)) / max(alpha-min(alpha))))
2998. # data_binned <- as.data.frame(data_binned %>% group_by(patient) %>% mutate(Ntheta = (theta-min(theta)) / max(theta-min(theta))))
2999. # data_binned <- as.data.frame(data_binned %>% group_by(patient) %>% mutate(Ndelta = (delta-min(delta)) / max(delta-min(delta))))
3000. # data_binned <- as.data.frame(data_binned %>% group_by(patient) %>% mutate(Nbeta = (beta-min(beta)) / max(beta-min(beta))))
3001. # data_binned_melt <-melt(data_binned[c("Ndelta","Ntheta","Nalpha","Nbeta","patient","bin")], id=c("patient","bin"))
3002. # data_binned_melt$variable <- factor(data_binned_melt$variable, labels = c("\u03B4 power (1-4Hz)","\u03B8 power (5-7Hz)","\u03B1 power (8-14Hz)","\u03B2 power (15-25Hz)"))
3003. # data_binned_melt$patient <- factor(data_binned_melt$patient, levels = c(8:1))
3004. #
3005. # pdf(file="//lexport/iss02.charpier/analyses/stephen.whitmarsh/images/hspike/R/polar_power_combined.pdf")
3006. #
3007. # ggplot(data=data_binned_melt, aes(x = bin, y = value, fill=patient, col=patient)) +
3008. # ggtitle("Normalized circadian LFP power") +
3009. # scale_fill_brewer(palette = "Set2") + scale_color_brewer(palette = "Set2") +
3010. # geom_vline(xintercept = seq(0.5, 21.5 * bindivision, by = 3 * bindivision), colour = "grey90") +
3011. # geom_hline(yintercept = seq(1, 8, by = 1), colour = "grey90") +
3012. # geom_bar(stat="identity", width = 1) +
3013. # theme_bw() +
3014. # theme(panel.border = element_blank(),
3015. # legend.key = element_blank(),
3016. # axis.ticks = element_blank(),
3017. # axis.text.y = element_blank(),
3018. # panel.grid = element_blank(),
3019. # axis.title.x = element_blank(),
3020. # axis.title.y = element_blank(),
3021. # strip.background = element_blank(),
3022. # ) +
3023. # scale_x_continuous(breaks=seq(0.5, 21.5*bindivision, by = 3 * bindivision),
3024. # labels = c("0" = "00:00", "3" = "3:00", "6" = "6:00", "9" = "9:00",
3025. # "12" = "12:00", "15" = "15:00", "18" = "18:00", "21" = "21:00")) +
3026. # coord_polar(theta = "x", start = 0) +
3027. # ylim(-1, 7) +
3028. # facet_wrap(~variable)
3029. #
3030. # dev.off()
3031. #
3032. # library(circular)
3033. #
3034. # data$radian = data$minute / (24*60) * pi * 2
3035. #
3036. # stats_circ <- median.circular(data$radian)
3037. #
3038. #
3039. # # load data: Patients x Units x time window
3040. # data_IED <- read.csv(file="//lexport/iss02.charpier/analyses/stephen.whitmarsh/data/hspike/hypdata_table.txt", sep=',', header=TRUE, dec='.', na.strings = " ")
3041. #
3042. # # prepare data
3043. # data_IED$patient <- factor(data_IED$patient, levels = c(8:1))
3044. # data_IED$hour <- data_IED$minute / (60)
3045. # data_IED$part <- factor(data_IED$part,)
3046. #
3047. # # load data: Patients x Units x time window
3048. # data_seizures <- read.csv(file="//lexport/iss02.charpier/analyses/stephen.whitmarsh/data/hspike/seizuredata_table.txt", sep=',', header=TRUE, dec='.', na.strings = " ")
3049. #
3050. # # prepare data
3051. # data_seizures$patient <- factor(data_seizures$patient, levels = c(8:1))
3052. # data_seizures$hour <- data_seizures$minute / (60)
3053. #
3054. # # graphics.off()
3055. #
3056. # # display.brewer.all(colorblindFriendly = TRUE)
3057. #
3058. # pdf(file="//lexport/iss02.charpier/analyses/stephen.whitmarsh/images/hspike/R/polar_density_combined.pdf")
3059. #
3060. # ggplot(data=data_IED, aes(x=hour, fill=patient, col=patient)) +
3061. # ggtitle("Circadian interictal density & seizure occurance") +
3062. # scale_fill_brewer(palette = "Set2") + scale_color_brewer(palette = "Set2") +
3063. # geom_vline(xintercept = seq(0, 21, by = 3), colour = "grey90") +
3064. # geom_hline(yintercept = seq(-0.25+0.025, -0.05, by = 0.025), colour = "grey90") +
3065. # geom_histogram(position = 'stack', aes(y = stat(density)), binwidth=0.1, center = 0.05) + # make sure binwidth is multiple of 24/60
3066. # coord_polar(theta = "x") +
3067. # scale_x_continuous(breaks=seq(0, 21, by = 3), labels = c("0" = "00:00", "3" = "3:00", "6" = "6:00", "9" = "9:00", "12" = "12:00", "15" = "15:00", "18" = "18:00", "21" = "21:00")) +
3068. # coord_polar(theta = "x", start = 0, direction = 1) +
3069. # theme_bw() +
3070. # theme(panel.border = element_blank(),
3071. # legend.key = element_blank(),
3072. # axis.ticks = element_blank(),
3073. # axis.text.y = element_blank(),
3074. # panel.grid = element_blank(),
3075. # ) +
3076. # geom_point(size=1, data = data_seizures, aes(x = hour, y = -0.05 - as.numeric(patient)*0.025+0.025, fill = patient, col = patient)) +
3077. # ylim(-0.3, 0.85)
3078. #
3079. # dev.off()
3080. #
3081. #
3082. #
3083. #
3084. #
3085. #
3086. #
3087. #
3088. # ggplot(data=data_IED, aes(x=hour, fill=patient, col=patient)) +
3089. # ggtitle("Circadian interictal density & seizure occurance") +
3090. # scale_fill_brewer(palette = "Set2") + scale_color_brewer(palette = "Set2") +
3091. # geom_vline(xintercept = seq(0, 21, by = 3), colour = "grey90") +
3092. # #geom_segment(aes(x = 0, xend = 0, y = 1, yend = 0.1)) +
3093. # geom_hline(yintercept = seq(-0.25+0.025, -0.05, by = 0.025), colour = "grey90") +
3094. # geom_histogram(position = 'stack', aes(y = stat(density)), binwidth=0.1, center = 0.05) + # make sure binwidth is multiple of 24/60
3095. # coord_polar(theta = "x") +
3096. # scale_x_continuous(breaks=seq(0, 21, by = 3), labels = c("0" = "00:00", "3" = "03:00", "6" = "06:00", "9" = "09:00", "12" = "12:00", "15" = "15:00", "18" = "18:00", "21" = "21:00")) +
3097. # coord_polar(theta = "x", start = 0, direction = 1) +
3098. # theme_bw() +
3099. # theme(panel.border = element_blank(),
3100. # legend.key = element_blank(),
3101. # axis.ticks = element_blank(),
3102. # axis.text.y = element_blank(),
3103. # panel.grid = element_blank(),
3104. # ) +
3105. # geom_point(size=0.5, data = data_seizures, aes(x = hour, y = -0.05 - as.numeric(patient)*0.025+0.025, fill = patient, col = patient)) +
3106. # ylim(-0.3, 0.85)
3107. #
3108. #
3109. #
3110. #
3111. #
3112. #
3113. #
3114. #
3115. #
3116. #
3117. #
3118. #
3119. #
3120. # library(circular)
3121. # library(units)
3122. #
3123. # plot.circular(data_seizures$radian, pch = 16, cex = 1, stack = TRUE,
3124. # axes = TRUE, start.sep=0, sep = 0.025, shrink = 1,
3125. # bins = 23, ticks = TRUE, tcl = 0.025, tcl.text = 0.125,
3126. # col = NULL, tol = 0.04, uin = NULL,
3127. # xlim = c(-1, 1), ylim = c(-1, 1), digits = 2, units = "rad",
3128. # template = NULL, zero = 0.1, rotation = NULL,
3129. # main = NULL, sub=NULL, xlab = "", ylab = "")
3130. #
3131. #
3132. # plot(sin, -pi, 2*pi) # see ?plot.function
3133. #
3134. # pi_rad <- as_units(pi, "radians")
3135. #
3136. #
3137. #
3138. # pdf(file="//lexport/iss02.charpier/analyses/stephen.whitmarsh/images/hspike/R/polar_density_per patient.pdf")
3139. #
3140. # ggplot(data=data, aes(x=hour, fill=patient, col=patient)) +
3141. # ggtitle("Circadian interictal density per patient") +
3142. # scale_fill_brewer(palette = "Set2") + scale_color_brewer(palette = "Set2") +
3143. # geom_histogram(position = 'stack', aes(y = stat(density)), binwidth=0.2, center = 0.1) + # make sure binwidth is multiple of 24/60
3144. # coord_polar(theta = "x") +
3145. # geom_vline(xintercept = seq(0, 21, by = 3), colour = "grey90") +
3146. # scale_x_continuous(breaks=seq(0, 21, by = 3), labels = c("0" = "00:00", "3" = "03:00", "6" = "06:00", "9" = "09:00", "12" = "12:00", "15" = "15:00", "18" = "18:00", "21" = "21:00")) +
3147. # coord_polar(theta = "x", start = 0, direction = 1) +
3148. # theme_bw() +
3149. #
3150. # theme(panel.border = element_blank(),
3151. # legend.key = element_blank(),
3152. # axis.ticks = element_blank(),
3153. # axis.text.y = element_blank(),
3154. # panel.grid = element_blank(),
3155. # strip.background = element_blank(),
3156. # ) + facet_wrap(~patient)
3157. # dev.off()
3158. #
3159. #
3160. #
3161. #
3162. #
3163. #
3164. # # get ylim and xlim from plot
3165. # #layer_scales(p)$y$range$range
3166. # #layer_scales(p)$x$range$range
3167. #
3168. #
3169. #
3170. #
3171. #
3172. # # data$SU <- data$percRPV < 0.25
3173. # data$SU <- factor(ifelse(data$percRPV > 0.25, "MUA", "SUA"))
3174. # data$SU <- c("MUA","SUA","MUA","MUA","MUA","MUA","SUA","MUA","MUA","MUA","MUA","SUA","MUA","MUA","SUA","MUA","MUA","MUA","SUA","MUA","SUA","SUA","MUA","SUA","MUA")
3175. # # data$PI <- factor(data$PI)
3176. # data$PI <- factor(ifelse(data$template_tp < 550, "Int", "Pyr"))
3177. # # data$PI <- revalue(data$PI, c("Int"="Interneuron", "Pyr"="Pyramidal cell"))
3178. #
3179. #
3180. #
3181. #
3182. # data <- read.csv(file="//lexport/iss02.charpier/analyses/stephen.whitmarsh/data/aurelie/statAH.csv", sep=';', header=TRUE, dec=',', na.strings = " ")
3183. # data$P <- data$ï..Patient
3184. #
3185. # # get some p-values
3186. # detach(package:lmerTest)
3187. # library(lmerTest)
3188. # l1 <- lmer(EEG ~ NSE +S100 + (1 | P) + (1 | Time), data); summary(l1)
3189. # l1 <- lmer(EEG ~ NSE + S100 + (1 | P), data); summary(l1)
3190. #
3191. # l1 <- lm(EEG ~ NSE + S100, data); summary(l1)

analysis.R at commit 46bf8dc, under GPL-3.0 · at the source