30 April 2021 GC, R code to make output violin plot
#install.packages("RColorBrewer")
library("readxl")
library('data.table')
library("ggplot2")
library("SparseM")
library("quantreg")
library("plyr")
library("dplyr")
Attaching package: ‘SparseM’ The following object is masked from ‘package:base’: backsolve Attaching package: ‘dplyr’ The following objects are masked from ‘package:plyr’: arrange, count, desc, failwith, id, mutate, rename, summarise, summarize The following objects are masked from ‘package:data.table’: between, first, last The following objects are masked from ‘package:stats’: filter, lag The following objects are masked from ‘package:base’: intersect, setdiff, setequal, union
read_diff_SL <- function(path_in, sht) {
dt <- read_excel(path_in, sheet = sht)
dt <- dt[ -c(1, 4, 5, 6, 7, 8, 9, 21, 22)]
# rename columns
clm = colnames(dt)
clm= clm[-c(1, 2)]
newco <- c()
for (fb in clm){
f <- strsplit(fb, "_")[[1]][1]
newco <- c(newco, f)
}
newco <- c('Study_ID', 'Device', newco)
print(newco)
colnames(dt) <- newco
return(dt)
}
path1 <- '/media/guido/LACIE/Cingle_Guido/Master/Headband/Diff_SL.xlsx'
sht1 <- 'Diff_SL_65dB_270deg'
dsl <- read_diff_SL(path1, sht1)
New names: * `` -> ...1
[1] "Study_ID" "Device" "500" "630" "800" "1000" [7] "1250" "1600" "2000" "2500" "3150" "4000" [13] "5000"
bh5 <- dsl[dsl$Device == 'BAHA5P',]
bh5 <- subset(bh5, select = -c(Device))
sbj <- seq(from = 1, to = nrow(bh5))
bh5$BAHA5P_group <- sbj
bh5_long <- melt(setDT(bh5), id.vars = c('Study_ID', 'BAHA5P_group'),
value.name = 'Diff_SL', variable.name = 'Frequency_band')
bh5_long
Study_ID | BAHA5P_group | Frequency_band | Diff_SL |
---|---|---|---|
<dbl> | <int> | <fct> | <dbl> |
21 | 1 | 500 | -34 |
38 | 2 | 500 | -41 |
48 | 3 | 500 | -57 |
52 | 4 | 500 | -34 |
54 | 5 | 500 | -41 |
55 | 6 | 500 | -48 |
56 | 7 | 500 | -46 |
57 | 8 | 500 | -56 |
58 | 9 | 500 | -23 |
59 | 10 | 500 | -55 |
60 | 11 | 500 | -47 |
61 | 12 | 500 | -43 |
62 | 13 | 500 | -53 |
63 | 14 | 500 | -51 |
64 | 15 | 500 | -36 |
65 | 16 | 500 | -47 |
66 | 17 | 500 | -43 |
67 | 18 | 500 | -48 |
68 | 19 | 500 | -49 |
69 | 20 | 500 | -53 |
70 | 21 | 500 | -49 |
71 | 22 | 500 | -49 |
72 | 23 | 500 | -53 |
73 | 24 | 500 | -49 |
74 | 25 | 500 | -52 |
75 | 26 | 500 | -60 |
76 | 27 | 500 | -58 |
77 | 28 | 500 | -40 |
78 | 29 | 500 | -43 |
79 | 30 | 500 | -49 |
⋮ | ⋮ | ⋮ | ⋮ |
56 | 7 | 5000 | -14 |
57 | 8 | 5000 | -38 |
58 | 9 | 5000 | -32 |
59 | 10 | 5000 | -51 |
60 | 11 | 5000 | -38 |
61 | 12 | 5000 | -36 |
62 | 13 | 5000 | -17 |
63 | 14 | 5000 | -41 |
64 | 15 | 5000 | -38 |
65 | 16 | 5000 | -25 |
66 | 17 | 5000 | -23 |
67 | 18 | 5000 | -45 |
68 | 19 | 5000 | -42 |
69 | 20 | 5000 | -23 |
70 | 21 | 5000 | -30 |
71 | 22 | 5000 | -44 |
72 | 23 | 5000 | -49 |
73 | 24 | 5000 | -46 |
74 | 25 | 5000 | -31 |
75 | 26 | 5000 | -1 |
76 | 27 | 5000 | -31 |
77 | 28 | 5000 | -28 |
78 | 29 | 5000 | -33 |
79 | 30 | 5000 | -42 |
80 | 31 | 5000 | -20 |
81 | 32 | 5000 | -38 |
82 | 33 | 5000 | -48 |
83 | 34 | 5000 | -50 |
84 | 35 | 5000 | -23 |
85 | 36 | 5000 | -44 |
q <- t(do.call("rbind", tapply(bh5_long$Diff_SL, bh5_long$Frequency_band, quantile, c(0.1, 0.5, 0.9))))
q <- as.data.frame(q)
p <- c("P10", "P50", "P90")
q$Percentile <- p
q
500 | 630 | 800 | 1000 | 1250 | 1600 | 2000 | 2500 | 3150 | 4000 | 5000 | Percentile | |
---|---|---|---|---|---|---|---|---|---|---|---|---|
<dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <chr> | |
10% | -58 | -47.5 | -38.0 | -38.0 | -34.5 | -30 | -31.5 | -28.0 | -27.5 | -32.5 | -48.5 | P10 |
50% | -49 | -39.5 | -29.5 | -28.5 | -26.0 | -24 | -22.0 | -20.5 | -19.0 | -21.0 | -38.0 | P50 |
90% | -38 | -29.0 | -24.0 | -21.0 | -19.0 | -18 | -12.5 | -12.5 | -7.0 | -3.5 | -20.0 | P90 |
P10 = q %>% slice(1)
P50 = q %>% slice(2)
P90 = q %>% slice(3)
P10_long <- melt(setDT(P10), id.vars = c('Percentile'),
value.name = 'Diff_SL', variable.name = 'Frequency_band')
P50_long <- melt(setDT(P50), id.vars = c('Percentile'),
value.name = 'Diff_SL', variable.name = 'Frequency_band')
P90_long <- melt(setDT(P90), id.vars = c('Percentile'),
value.name = 'Diff_SL', variable.name = 'Frequency_band')
# Make plots wider
options(repr.plot.width=8, repr.plot.height=6)
ttl = "BAHA5P-group (N=36): diff. SL BC vs. AC-path, ISTS 65 dB, at best-ear-side"
p <- ggplot(bh5_long, aes(x=Frequency_band, y=Diff_SL, group=Study_ID, colour=BAHA5P_group)) +
geom_line() + scale_colour_gradientn(colours=rainbow(5)) + theme_bw() +
geom_hline(yintercept=0, linetype="dotted", color = "red", size = 1.5)
p <- p + geom_line(P10_long, mapping = aes(x=Frequency_band, y=Diff_SL, group = Percentile), colour="black", size=1.5, linetype='dashed')
p <- p + geom_line(P50_long, mapping = aes(x=Frequency_band, y=Diff_SL, group = Percentile), colour="black", size=1.5, linetype='solid')
p <- p + geom_line(P90_long, mapping = aes(x=Frequency_band, y=Diff_SL, group = Percentile), colour="black", size=1.5, linetype='twodash')
p <- p + ggtitle("") + xlab("CF 1/3 octave band [Hz]") + ylab("Diff. SL BC vs. AC [dB]")
p <- p + theme(plot.margin = margin(1, 1, 1, 1, "cm"),
plot.title = element_text(colour="black", size=16, face="bold", hjust=0.5, vjust=1),
axis.title.x = element_text(colour="black", size=15, face="plain", hjust=0.5, margin = margin(t = .8, r = 0, b = 0, l = 0, 'cm')),
axis.title.y = element_text(colour="black", size=15, face="plain", hjust=0.5, vjust=0.5, margin = margin(t = 0, r = .8, b = 0, l = 0, 'cm')),
axis.text.x = element_text(colour = "black", size = 11, vjust=0),
axis.text.y = element_text(colour = "black", size = 11, hjust=0),
legend.position = "none")
#legend.position = c(0.8, 0.19)
#legend.background = element_rect(size=0.3, linetype="solid", colour ="black"),
#legend.title = element_text(colour="black", size=14, face="plain", margin=margin(0.25, 0.25, 0, 1.1, 'cm')),
#legend.text = element_text(colour="black", size=12, margin=margin(0.25, 0.25, 0.25, 0.25, 'cm')),
#legend.margin = margin(0.25, 0.25, 0.25, 0.25, 'cm'),
p
tiff("/home/guido/R/cingle/figures/diff_SL_270deg_BH5.tiff", units="in", width=8, height=6, res=300)
p
dev.off()