R/abr_final.R
In aballas22/autoabr:
Defines functions
LIF_V
LIF_III
LIF_I
get_stim_side
get_ms
get_intensity
get_age
fix_jewett_single
fix_jewett
find_amplitude
dB_to_V
create_area_analysis_df
create_analysis_df
change_jewetts
analyze
abr_test
abr_test <- function(patient){
#patient = list.files(pattern=paste("ID-",pat_id,sep=""))
patient <- sort(patient)
df <- data.frame("Files" = patient,stringsAsFactors = FALSE)
samples <- seq(0,14.98, by=0.03211991)
assign("samples",samples,envir = .GlobalEnv)
intensities <- c()
for(file in df$Files) {
intensity <- get_intensity(file)
intensities <- as.numeric(c(intensities, intensity))
}
df<-cbind(df,Intensity = intensities)
# Loop that find diff volumes
v = 1
volumes <- c()
for (i in 1:length(intensities)){
if (i==length(intensities)){
volumes[v] = intensities[i]
}
else {
if (intensities[i] == intensities[i+1]) {
v = v
}
else {
volumes[v] = intensities[i]
v = v + 1
}
}
}
valleys_list <- list()
peaks_list <- list()
avg_list <- list()
jewett_list <- list()
assign("volumes",volumes,envir = .GlobalEnv)
for(v in 1:length(volumes)) {
start <- head(grep(volumes[v],df$Intensity),n=1)
end <- tail(grep(volumes[v],df$Intensity),n=1)
amplitudes<-data.frame()
count = 1
sum = 0
for (i in start:end) {
arxeio<-df$Files[i]
amp<-find_amplitude(arxeio)
amplitudes<-rbind(amplitudes,amp)
sum = sum + amplitudes[count,]
count = count + 1
}
avg<-as.numeric(sum/(count-1))
peaks<-find_peaks(avg[1:220], m=4,error = 5)
valleys<-find_peaks(-avg[1:220], m=4)
#akrotata<-sort(c(peaks,valleys))
jewett<-as.roman(find_jewett(peaks))
#assign(paste("avg", v, sep = ""), avg, envir = .GlobalEnv)
#assign(paste("peaks", v, sep = ""), peaks, envir = .GlobalEnv)
#assign(paste("jewett", v, sep = ""), jewett, envir = .GlobalEnv)
#rm(amplitudes,sum,akrotata,amp,avg,count,peaks,valleys)
valleys_list[[v]]<-valleys
peaks_list[[v]] <- peaks
avg_list[[v]] <- avg
jewett_list[[v]] <- jewett
}
assign("valleys_list",valleys_list,envir = .GlobalEnv)
assign("peaks_list",peaks_list,envir = .GlobalEnv)
assign("avg_list",avg_list,envir = .GlobalEnv)
assign("jewett_list",jewett_list,envir = .GlobalEnv)
## Loop that fixes Jewett Labels
if (length(volumes) == 1) {
Jewett <- fix_jewett_single(peaks_list[[1]],jewett_list[[1]],avg_list[[1]])
}else if (length(volumes) == 2) {
Jewett <- fix_jewett(peaks_list[[1]],jewett_list[[1]],peaks_list[[2]],jewett_list[[2]])
}else if (length(volumes) == 3) {
Jewett <- list()
Jewett1 <- fix_jewett(peaks_list[[1]],jewett_list[[1]],peaks_list[[2]],jewett_list[[2]])
Jewett2 <- fix_jewett_single(peaks_list[[3]],jewett_list[[3]])
Jewett[[1]]<-Jewett1[[1]]
Jewett[[2]]<-Jewett1[[2]]
Jewett[[3]]<-Jewett2[[1]]
Jewett[[4]]<-Jewett1[[3]]
Jewett[[5]]<-Jewett1[[4]]
Jewett[[6]]<-Jewett2[[2]]
}else {
Jewett <- list()
peaks_temp <- list()
jewett_temp <-list()
for (v in 1:length(volumes)){
if (v == length(volumes)){
Jewett<-c(peaks_temp,jewett_temp)
}else {
j_temp <- fix_jewett(peaks_list[[v]],jewett_list[[v]],peaks_list[[v+1]],jewett_list[[v+1]])
peaks_temp[[v]] <- j_temp[[1]]
peaks_temp[[v+1]] <- j_temp[[2]]
jewett_temp[[v]] <- j_temp[[3]]
jewett_temp[[v+1]] <- j_temp[[4]]
}
}
}
if (length(volumes)<4) {
par(mfrow=c(length(volumes),1),xpd=TRUE)
}else{
par(mfrow=c(length(volumes)/2,2),xpd=TRUE)
}
assign("Jewett",Jewett,envir = .GlobalEnv)
samples <- seq(0,14.98, by=0.03211991)
## Creates ABR Plots
for (v in 1:length(volumes)){
l <- length(Jewett)
a <- avg_list[[v]]
p <- Jewett[[v]]
lab <- Jewett[[((l/2)+v)]]
plot(samples[1:225],a[1:225],type='l',yaxt='n',main = paste(volumes[v], 'nHL'),
col='red',panel.first = grid(),
xlab="Time (ms)", ylab="Amplitude (μV)")
points(samples[p],a[p])
plot_without = recordPlot()
assign("plot_without",plot_without,envir = .GlobalEnv)
text(samples[p],a[p],labels = lab, pos = 3)
}
plot = recordPlot()
assign("plot",plot,envir = .GlobalEnv)
return(plot)
}
analyze <- function(Jewett, avg_list, name) {
samples <- seq(0,14.98, by=0.03211991)
avg<-avg_list[[1]]
peaks <- Jewett[[1]]
name <- name
amplitudes <-c()
for(i in 1:length(peaks)) {
if (!is.na(peaks[i])) {
amp <- avg[peaks[i]]
}else {
amp <- NA
}
amplitudes <- c(amplitudes,amp)
}
latencies <- c()
for(i in 1:length(peaks)) {
if (!is.na(peaks[i])) {
lat <- samples[peaks[i]]
}else {
lat <- NA
}
latencies <- c(latencies,lat)
}
inter_latencies <-c()
for (i in 1:(length(latencies)-1)) {
for (j in (i+1):length(latencies)){
if (!is.na(latencies[i])) {
inter_lat <- latencies[j] - latencies[i]
}
else {
inter_lat <- NA
}
inter_latencies <- c(inter_latencies,inter_lat)
}
}
latency_analysis <- c(name,amplitudes,latencies,inter_latencies)
return(latency_analysis)
}
change_jewetts <- function(x) {
if (length(grep(x,c('Yes' , 'yes' , 'y' , 'Y' , 'YES')))!= 0) {
points <- dlg_input(message = "Enter num of labels you would like to mark (digit): ")
msgBox <- tkmessageBox(title = "Information",
message = paste("After you press OK, click max ", points$res," labels",sep=""), icon = "info", type = "ok")
replayPlot(plot_without)
new_peaks <- identify(samples[1:225],avg_list[[1]][1:225],n=as.numeric(points$res))
jewetts_new<-as.roman(find_jewett(new_peaks))
Jewett[[1]]<-new_peaks
Jewett[[2]]<-jewetts_new
# analyze latencies of changes labels
pat_id <- strsplit(patient[1]," ")[[1]][3]
name <- paste(pat_id,"_",volumes[1],"_dB",sep="")
latency_analysis <- analyze(Jewett,avg_list,name)
analysis_data[nrow(analysis_data) + 1,] <- latency_analysis
for (v in 1:length(volumes)){
l <- length(Jewett)
a <- avg_list[[v]]
p <- Jewett[[v]]
lab <- Jewett[[((l/2)+v)]]
plot(samples[1:225],a[1:225],type='l',yaxt='n',main = paste(volumes[v], 'nHL'),
col='red',panel.first = grid(),
xlab="Time (ms)", ylab="Amplitude (μV)")
points(samples[p],a[p])
text(samples[p],a[p],labels = lab, pos = 3)
}
final_plot = recordPlot()
}else {
# analyze latencies of unchanged labels
pat_id <- strsplit(patient[1]," ")[[1]][3]
name <- paste(pat_id,"_",volumes[1],"_dB",sep="")
latency_analysis <- analyze(Jewett,avg_list,name)
analysis_data[nrow(analysis_data) + 1,] <- latency_analysis
final_plot <- plot
}
}
create_analysis_df <- function(){
names <- c('Patient','I','II','III','IV','V','latency_I','latency_II','latency_III','latency_IV','latency_V','inter_I_II','inter_I_III','inter_I_IV','inter_I_V')
names <- c(names,'inter_II_III','inter_II_IV','inter_II_V','inter_III_IV','inter_III_V','inter_IV_V')
analysis_data<-data.frame(matrix(ncol=length(names)),nrow=0)
colnames(analysis_data)<-names
analysis_data <- analysis_data[-1,] # deletes 1st row
analysis_data <- analysis_data[-22] # deletes 22nd column
return(analysis_data)
}
create_area_analysis_df <- function() {
names<-c('Patient','Curve I','Curve II','Curve III','Curve IV','Curve V')
area_analysis <- data.frame(matrix(ncol=length(names)),nrow=0)
colnames(area_analysis)<-names
area_analysis <- area_analysis[-1,]
area_analysis <- area_analysis[-7]
return(area_analysis)
}
dB_to_V <- function(value) {
value <- 10^(value/20)
return(value)
}
fill_peaks <- function (Jewett) {
for (i in 1:5) {
if (is.na(match(as.roman(i),Jewett[[2]]))) {
x<-append(Jewett[[1]],NA,after=(i-1))
y<-append(Jewett[[2]],as.roman(i),after=(i-1))
}else {
x<- Jewett[[1]]
y<- Jewett[[2]]
}
Jewett[[1]] <- x
Jewett[[2]] <- as.roman(y)
}
return(Jewett)
}
find_amplitude <- function(arxeio){
file = xmlParse(arxeio)
xml_data <- xmlToList(file)
stimuli_side <- xml_data[["Waveform"]][[".attrs"]][["StimuliSide"]]
samples_num <- as.numeric(xml_data[["Waveform"]][["NumberOfSamples"]])
samples <- seq(0,14.98, by=0.03211991)
IPSI_A_Raw <- as.numeric(xml_data[["Waveform"]][["Response"]][["IPSI_A_Raw"]])
IPSI_B_Raw <- as.numeric(xml_data[["Waveform"]][["Response"]][["IPSI_B_Raw"]])
IPSI_Avg <- (IPSI_A_Raw + IPSI_B_Raw)/2
###### Metatropi Raw data se μV
gain <- strsplit(xml_data[["Waveform"]][["Gain"]]," ")[[1]]
gain <- as.numeric(gain[1])
gain <- dB_to_V(gain)
data_res <- (1.6/gain)/32768
amp <- (IPSI_Avg * data_res)*10^6
return(amp)
}
find_jewett <- function (x){
jewett <-c()
samples <- seq(0,14.98, by=0.03211991)
for (i in 1:length(x)) {
if ((0.8< samples[x[i]]) & (samples[x[i]] <= 2) ) {
jewett[i]<- 1
}
else if (2< samples[x[i]] & (samples[x[i]] <= 3)) {
jewett[i]<- 2
}
else if (3 < samples[x[i]] & (samples[x[i]] <= 4)) {
jewett[i]<- 3
}
else if (4 < samples[x[i]] & (samples[x[i]] <= 4.9)) {
jewett[i]<- 4
}
else if (4.9< samples[x[i]] & (samples[x[i]] <= 6.8)) {
jewett[i]<- 5
}
else {
jewett[i]<- 0
}
}
return (jewett)
}
find_peaks <- function (x, m = 3, error = 5){
shape <- diff(sign(diff(x, na.pad = FALSE)))
p <- sapply(which(shape < 0), FUN = function(i){
z <- i - m + 1
z <- ifelse(z > 0, z, 1)
w <- i + m + 1
w <- ifelse(w < length(x), w, length(x))
if(all(x[c(z : i, (i + 2) : w)] <= x[i + 1])) return(i + 1) else return(numeric(0))
})
p <- unlist(p)
# Ypologismos twn peaks me sfalma error
pks<-c()
counter =1
for (j in 1:length(p)) {
if(j == length(p)){
pks[counter]<- p[j]
}
else if (abs(p[j+1] - p[j]) > error){
pks[counter] <- p[j]
counter = counter + 1
}
else {
pks[counter] <- p[j]
next
counter = counter + 1
}
}
pks
}
fix_jewett <- function(peaks1,jewett1,peaks2,jewett2) {
norm <- c(1.5,2.5,3.5,4.5,5.5)
samples <- seq(0,14.98, by=0.03211991)
################## & (length(is.na(jewett2[duplicated(jewett2)]))<1)
if (length(grep(TRUE,duplicated(jewett2,incomparables=NA)))>0 & length(grep(TRUE,duplicated(jewett1,incomparables=NA)))==0 ){
duplicate <- jewett2[duplicated(jewett2,incomparables=NA)]
dup_num <- as.numeric((jewett2[duplicated(jewett2,incomparables=NA)]))
dup_pos<-as.numeric(grep(TRUE,duplicated(jewett2,incomparables=NA)))
test_peak <- c(peaks2[dup_pos-1],peaks2[dup_pos])
if (grep(duplicate,jewett1)) {
corr_pos<-grep(duplicate,jewett1)
corr_peak <- peaks1[corr_pos]
if (abs(corr_peak-test_peak[1]) <= abs(corr_peak-test_peak[2])) {
peaks2 <- peaks2[!peaks2 %in% test_peak[2]]
jewett2 <- jewett2[-dup_pos]
}else{
peaks2 <- peaks2[!peaks2 %in% test_peak[1]]
jewett2 <- jewett2[-(dup_pos-1)]
}
}else {
if (abs(samples[test_peak[1]]-norm[dup_num]) < abs(samples[test_peak[2]]-norm[dup_num])){
peaks2 <- peaks2[!peaks2 %in% test_peak[2]]
jewett2 <- jewett2[-dup_pos]
}else{
peaks2 <- peaks2[!peaks2 %in% test_peak[1]]
jewett2 <- jewett2[-(dup_pos-1)]
}
}
} #1o else if
##################
else if (length(grep(TRUE,duplicated(jewett1,incomparables=NA)))>0 & length(grep(TRUE,duplicated(jewett2,incomparables=NA)))==0 ){
duplicate <- jewett1[duplicated(jewett1,incomparables=NA)]
dup_num <- as.numeric(jewett1[duplicated(jewett1,incomparables=NA)])
dup_pos<-as.numeric(grep(TRUE,duplicated(jewett1,incomparables=NA)))
test_peak <- c(peaks1[dup_pos-1],peaks1[dup_pos])
if (grep(duplicate,jewett2)) {
corr_pos<-grep(duplicate,jewett2)
corr_peak <- peaks2[corr_pos]
if(abs(corr_peak-test_peak[1]) <= abs(corr_peak-test_peak[2])) {
peaks1 <- peaks1[!peaks1 %in% test_peak[2]]
jewett1 <- jewett1[-dup_pos]
}else {
peaks1 <- peaks1[!peaks1 %in% test_peak[1]]
jewett1 <- jewett1[-(dup_pos-1)]
}
}else{
if (abs(samples[test_peak[1]]-norm[dup_num]) < abs(samples[test_peak[2]]-norm[dup_num])){
peaks1 <- peaks1[!peaks1 %in% test_peak[2]]
jewett1 <- jewett1[-dup_pos]
}else{
peaks1 <- peaks1[!peaks2 %in% test_peak[1]]
jewett1 <- jewett1[-(dup_pos-1)]
}
}
}# eksw else if
##################
else if (length(grep(TRUE,duplicated(jewett1,incomparables=NA)))>0 & length(grep(TRUE,duplicated(jewett2,incomparables=NA)))>0){
j1<-fix_jewett_single(peaks1,jewett1)
j2<-fix_jewett_single(peaks2,jewett2)
peaks1 <- j1[[1]]
peaks2 <- j2[[1]]
jewett1 <- j1[[2]]
jewett2 <- j2[[2]]
}
else{
peaks1<-peaks1
peaks2<-peaks2
jewett1<-jewett1
jewett2<-jewett2
}# teleutaio else
Jewetts <- list(peaks1,peaks2,jewett1,jewett2)
return(Jewetts)
}
fix_jewett_single <- function(peaks,jewett,avg) {
norm <- c(1.5,2.5,3.5,4.5,5.5)
samples <- seq(0,14.98, by=0.03211991)
if (length(grep(TRUE,duplicated(jewett,incomparables=NA)))>0) {
duplicate <- jewett[duplicated(jewett,incomparables=NA)]
dup_num <- as.numeric((jewett[duplicated(jewett,incomparables=NA)]))
dup_pos<-as.numeric(grep(TRUE,duplicated(jewett,incomparables=NA)))
test_peak <- c(peaks[dup_pos-1],peaks[dup_pos])
if ((abs(samples[test_peak[1]]-norm[dup_num]) < abs(samples[test_peak[2]]-norm[dup_num])) & (avg[test_peak[1]] > avg[test_peak[2]])) {
# afairw teleutaio
peaks <- peaks[!peaks %in% test_peak[2]]
jewett <- jewett[-dup_pos]
}
else if ((abs(samples[test_peak[1]]-norm[dup_num]) > abs(samples[test_peak[2]]-norm[dup_num])) & (avg[test_peak[1]] < avg[test_peak[2]])) {
peaks <- peaks[!peaks %in% test_peak[1]]
jewett <- jewett[-(dup_pos-1)]
}
else {
if (avg[test_peak[1]] > avg[test_peak[2]]) {
peaks <- peaks[!peaks %in% test_peak[2]]
jewett <- jewett[-dup_pos]
} else {
# afairw proteleutaio
peaks <- peaks[!peaks %in% test_peak[1]]
jewett <- jewett[-(dup_pos-1)]
}
}
}
else {
peaks<-peaks
jewett<-jewett
}
Jewetts <- list(peaks,jewett)
# Deletes NA and corresponding peaks from Jewett List
while(length(grep(TRUE,is.na(Jewetts[[2]])))>=1) {
na <- grep(TRUE,is.na(Jewetts[[2]]))
for (n in na){
Jewetts[[1]]<-Jewetts[[1]][-n]
Jewetts[[2]]<-Jewetts[[2]][-n]
}
}
# Deletes remaining peaks if there is a diff
if (length(Jewetts[[1]]) > length(Jewetts[[2]])) {
diff <- length(Jewetts[[1]]) - length(Jewetts[[2]])
for (d in diff) {
Jewetts[[1]] <- Jewetts[[1]][-(length(Jewetts[[1]]))]
}
}
# Fills 5 peaks with NA
for (i in 1:5) {
if (is.na(match(as.roman(i),Jewetts[[2]]))) {
x<-append(Jewetts[[1]],NA,after=(i-1))
y<-append(Jewetts[[2]],as.roman(i),after=(i-1))
}else {
x<- Jewetts[[1]]
y<- Jewetts[[2]]
}
Jewetts[[1]] <- x
Jewetts[[2]] <- as.roman(y)
}
return(Jewetts)
}
get_age <-function(arxeio2) {
file2 = xmlParse(arxeio2)
xml_data2 <- xmlToList(file2)
year<-strsplit(xml_data2[["Waveform"]][[".attrs"]][["DayOfBirth"]],"-")[[1]]
year<-as.numeric(year[1])
age<-2020-year
return(age)
}
get_intensity <- function(arxeio) {
file = xmlParse(arxeio)
xml_data <- xmlToList(file)
intensity <- xml_data[["Waveform"]][[".attrs"]][["Intensity"]]
return(intensity)
}
get_ms <- function(peak) {
samples <- seq(0,14.98, by=0.03211991)
ms <- samples[peak]
return(ms)
}
get_stim_side<- function(arxeio) {
file = xmlParse(arxeio)
xml_data <- xmlToList(file)
stimuli_side <- xml_data[["Waveform"]][[".attrs"]][["StimuliSide"]]
return(stimuli_side)
}
LIF_I <- function(intensity,latency_I) {
if (is.na(latency_I)) {
y1_min<-c(0.67,0.19,0.45,1.09)
y1_max<-c(2.05,1.96,1.927,2)
x<-c(70,80,90,100)
plot(1, type="n", xlim=c(60,110), ylim=c(0,10),
main = 'Wave I Latency / Intensity Function',
panel.first = grid(),
xlab="Intensity (dB nHL)", ylab="Time (ms)")
polygon(c(x,x[4:1]),c(y1_min,y1_max[4:1]),col='grey')
legend(x='topright',legend='Wave I peak is not present')
lat_I_plot<-recordPlot()
}else {
y1_min<-c(0.67,0.19,0.45,1.09)
y1_max<-c(2.05,1.96,1.927,2)
x<-c(70,80,90,100)
plot(1, type="n", xlim=c(60,110), ylim=c(0,10),
main = 'Wave I Latency / Intensity Function',
panel.first = grid(),
xlab="Intensity (dB nHL)", ylab="Time (ms)")
polygon(c(x,x[4:1]),c(y1_min,y1_max[4:1]),col='grey')
points(intensity,latency_I,pch=21,col='black',bg='red')
lat_I_plot<-recordPlot()
}
return(lat_I_plot)
}
LIF_III <- function(intensity,latency_III) {
if(is.na(latency_III)){
y3_min<-c(3.24,2.38,2.6,3.24)
y3_max<-c(3.8,3.95,3.9,3.8)
x<-c(70,80,90,100)
plot(1, type="n", xlim=c(60,110), ylim=c(0,10),
main = 'Wave III Latency / Intensity Function',
panel.first = grid(),
xlab="Intensity (dB nHL)", ylab="Time (ms)")
polygon(c(x,x[4:1]),c(y3_min,y3_max[4:1]),col='grey')
legend(x='topright',legend='Wave III peak is not present')
lat_III_plot<-recordPlot()
}else{
y3_min<-c(3.24,2.38,2.6,3.24)
y3_max<-c(3.8,3.95,3.9,3.8)
x<-c(70,80,90,100)
plot(1, type="n", xlim=c(60,110), ylim=c(0,10),
main = 'Wave III Latency / Intensity Function',
panel.first = grid(),
xlab="Intensity (dB nHL)", ylab="Time (ms)")
polygon(c(x,x[4:1]),c(y3_min,y3_max[4:1]),col='grey')
points(intensity,latency_III,pch=21,col='black',bg='red')
lat_III_plot<-recordPlot()
}
return(lat_III_plot)
}
LIF_V <- function(intensity,latency_V) {
if(is.na(latency_V)){
y5_min<-c(4.3,4.14,4.5,4.6)
y5_max<-c(6.2,6.8,6.97,7)
x<-c(70,80,90,100)
plot(1, type="n", xlim=c(60,110), ylim=c(0,10),
main = 'Wave V Latency / Intensity Function',
panel.first = grid(),
xlab="Intensity (dB nHL)", ylab="Time (ms)")
polygon(c(x,x[4:1]),c(y5_min,y5_max[4:1]),col='grey')
legend(x='topright',legend='Wave III peak is not present')
lat_V_plot<-recordPlot()
}else{
y5_min<-c(4.3,4.14,4.5,4.6)
y5_max<-c(6.2,6.8,6.97,7)
x<-c(70,80,90,100)
plot(1, type="n", xlim=c(60,110), ylim=c(0,10),
main = 'Wave V Latency / Intensity Function',
panel.first = grid(),
xlab="Intensity (dB nHL)", ylab="Time (ms)")
polygon(c(x,x[4:1]),c(y5_min,y5_max[4:1]),col='grey')
points(intensity,latency_V,pch=21,col='black',bg='red')
lat_V_plot<-recordPlot()
}
return(lat_V_plot)
}
aballas22/autoabr documentation built on June 15, 2020, 12:37 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions LIF_V LIF_III LIF_I get_stim_side get_ms get_intensity get_age fix_jewett_single fix_jewett find_amplitude dB_to_V create_area_analysis_df create_analysis_df change_jewetts analyze abr_test
abr_test <- function(patient){
#patient = list.files(pattern=paste("ID-",pat_id,sep=""))
patient <- sort(patient)
df <- data.frame("Files" = patient,stringsAsFactors = FALSE)
samples <- seq(0,14.98, by=0.03211991)
assign("samples",samples,envir = .GlobalEnv)
intensities <- c()
for(file in df$Files) {
intensity <- get_intensity(file)
intensities <- as.numeric(c(intensities, intensity))
}
df<-cbind(df,Intensity = intensities)
# Loop that find diff volumes
v = 1
volumes <- c()
for (i in 1:length(intensities)){
if (i==length(intensities)){
volumes[v] = intensities[i]
}
else {
if (intensities[i] == intensities[i+1]) {
v = v
}
else {
volumes[v] = intensities[i]
v = v + 1
}
}
}
valleys_list <- list()
peaks_list <- list()
avg_list <- list()
jewett_list <- list()
assign("volumes",volumes,envir = .GlobalEnv)
for(v in 1:length(volumes)) {
start <- head(grep(volumes[v],df$Intensity),n=1)
end <- tail(grep(volumes[v],df$Intensity),n=1)
amplitudes<-data.frame()
count = 1
sum = 0
for (i in start:end) {
arxeio<-df$Files[i]
amp<-find_amplitude(arxeio)
amplitudes<-rbind(amplitudes,amp)
sum = sum + amplitudes[count,]
count = count + 1
}
avg<-as.numeric(sum/(count-1))
peaks<-find_peaks(avg[1:220], m=4,error = 5)
valleys<-find_peaks(-avg[1:220], m=4)
#akrotata<-sort(c(peaks,valleys))
jewett<-as.roman(find_jewett(peaks))
#assign(paste("avg", v, sep = ""), avg, envir = .GlobalEnv)
#assign(paste("peaks", v, sep = ""), peaks, envir = .GlobalEnv)
#assign(paste("jewett", v, sep = ""), jewett, envir = .GlobalEnv)
#rm(amplitudes,sum,akrotata,amp,avg,count,peaks,valleys)
valleys_list[[v]]<-valleys
peaks_list[[v]] <- peaks
avg_list[[v]] <- avg
jewett_list[[v]] <- jewett
}
assign("valleys_list",valleys_list,envir = .GlobalEnv)
assign("peaks_list",peaks_list,envir = .GlobalEnv)
assign("avg_list",avg_list,envir = .GlobalEnv)
assign("jewett_list",jewett_list,envir = .GlobalEnv)
## Loop that fixes Jewett Labels
if (length(volumes) == 1) {
Jewett <- fix_jewett_single(peaks_list[[1]],jewett_list[[1]],avg_list[[1]])
}else if (length(volumes) == 2) {
Jewett <- fix_jewett(peaks_list[[1]],jewett_list[[1]],peaks_list[[2]],jewett_list[[2]])
}else if (length(volumes) == 3) {
Jewett <- list()
Jewett1 <- fix_jewett(peaks_list[[1]],jewett_list[[1]],peaks_list[[2]],jewett_list[[2]])
Jewett2 <- fix_jewett_single(peaks_list[[3]],jewett_list[[3]])
Jewett[[1]]<-Jewett1[[1]]
Jewett[[2]]<-Jewett1[[2]]
Jewett[[3]]<-Jewett2[[1]]
Jewett[[4]]<-Jewett1[[3]]
Jewett[[5]]<-Jewett1[[4]]
Jewett[[6]]<-Jewett2[[2]]
}else {
Jewett <- list()
peaks_temp <- list()
jewett_temp <-list()
for (v in 1:length(volumes)){
if (v == length(volumes)){
Jewett<-c(peaks_temp,jewett_temp)
}else {
j_temp <- fix_jewett(peaks_list[[v]],jewett_list[[v]],peaks_list[[v+1]],jewett_list[[v+1]])
peaks_temp[[v]] <- j_temp[[1]]
peaks_temp[[v+1]] <- j_temp[[2]]
jewett_temp[[v]] <- j_temp[[3]]
jewett_temp[[v+1]] <- j_temp[[4]]
}
}
}
if (length(volumes)<4) {
par(mfrow=c(length(volumes),1),xpd=TRUE)
}else{
par(mfrow=c(length(volumes)/2,2),xpd=TRUE)
}
assign("Jewett",Jewett,envir = .GlobalEnv)
samples <- seq(0,14.98, by=0.03211991)
## Creates ABR Plots
for (v in 1:length(volumes)){
l <- length(Jewett)
a <- avg_list[[v]]
p <- Jewett[[v]]
lab <- Jewett[[((l/2)+v)]]
plot(samples[1:225],a[1:225],type='l',yaxt='n',main = paste(volumes[v], 'nHL'),
col='red',panel.first = grid(),
xlab="Time (ms)", ylab="Amplitude (μV)")
points(samples[p],a[p])
plot_without = recordPlot()
assign("plot_without",plot_without,envir = .GlobalEnv)
text(samples[p],a[p],labels = lab, pos = 3)
}
plot = recordPlot()
assign("plot",plot,envir = .GlobalEnv)
return(plot)
}
analyze <- function(Jewett, avg_list, name) {
samples <- seq(0,14.98, by=0.03211991)
avg<-avg_list[[1]]
peaks <- Jewett[[1]]
name <- name
amplitudes <-c()
for(i in 1:length(peaks)) {
if (!is.na(peaks[i])) {
amp <- avg[peaks[i]]
}else {
amp <- NA
}
amplitudes <- c(amplitudes,amp)
}
latencies <- c()
for(i in 1:length(peaks)) {
if (!is.na(peaks[i])) {
lat <- samples[peaks[i]]
}else {
lat <- NA
}
latencies <- c(latencies,lat)
}
inter_latencies <-c()
for (i in 1:(length(latencies)-1)) {
for (j in (i+1):length(latencies)){
if (!is.na(latencies[i])) {
inter_lat <- latencies[j] - latencies[i]
}
else {
inter_lat <- NA
}
inter_latencies <- c(inter_latencies,inter_lat)
}
}
latency_analysis <- c(name,amplitudes,latencies,inter_latencies)
return(latency_analysis)
}
change_jewetts <- function(x) {
if (length(grep(x,c('Yes' , 'yes' , 'y' , 'Y' , 'YES')))!= 0) {
points <- dlg_input(message = "Enter num of labels you would like to mark (digit): ")
msgBox <- tkmessageBox(title = "Information",
message = paste("After you press OK, click max ", points$res," labels",sep=""), icon = "info", type = "ok")
replayPlot(plot_without)
new_peaks <- identify(samples[1:225],avg_list[[1]][1:225],n=as.numeric(points$res))
jewetts_new<-as.roman(find_jewett(new_peaks))
Jewett[[1]]<-new_peaks
Jewett[[2]]<-jewetts_new
# analyze latencies of changes labels
pat_id <- strsplit(patient[1]," ")[[1]][3]
name <- paste(pat_id,"_",volumes[1],"_dB",sep="")
latency_analysis <- analyze(Jewett,avg_list,name)
analysis_data[nrow(analysis_data) + 1,] <- latency_analysis
for (v in 1:length(volumes)){
l <- length(Jewett)
a <- avg_list[[v]]
p <- Jewett[[v]]
lab <- Jewett[[((l/2)+v)]]
plot(samples[1:225],a[1:225],type='l',yaxt='n',main = paste(volumes[v], 'nHL'),
col='red',panel.first = grid(),
xlab="Time (ms)", ylab="Amplitude (μV)")
points(samples[p],a[p])
text(samples[p],a[p],labels = lab, pos = 3)
}
final_plot = recordPlot()
}else {
# analyze latencies of unchanged labels
pat_id <- strsplit(patient[1]," ")[[1]][3]
name <- paste(pat_id,"_",volumes[1],"_dB",sep="")
latency_analysis <- analyze(Jewett,avg_list,name)
analysis_data[nrow(analysis_data) + 1,] <- latency_analysis
final_plot <- plot
}
}
create_analysis_df <- function(){
names <- c('Patient','I','II','III','IV','V','latency_I','latency_II','latency_III','latency_IV','latency_V','inter_I_II','inter_I_III','inter_I_IV','inter_I_V')
names <- c(names,'inter_II_III','inter_II_IV','inter_II_V','inter_III_IV','inter_III_V','inter_IV_V')
analysis_data<-data.frame(matrix(ncol=length(names)),nrow=0)
colnames(analysis_data)<-names
analysis_data <- analysis_data[-1,] # deletes 1st row
analysis_data <- analysis_data[-22] # deletes 22nd column
return(analysis_data)
}
create_area_analysis_df <- function() {
names<-c('Patient','Curve I','Curve II','Curve III','Curve IV','Curve V')
area_analysis <- data.frame(matrix(ncol=length(names)),nrow=0)
colnames(area_analysis)<-names
area_analysis <- area_analysis[-1,]
area_analysis <- area_analysis[-7]
return(area_analysis)
}
dB_to_V <- function(value) {
value <- 10^(value/20)
return(value)
}
fill_peaks <- function (Jewett) {
for (i in 1:5) {
if (is.na(match(as.roman(i),Jewett[[2]]))) {
x<-append(Jewett[[1]],NA,after=(i-1))
y<-append(Jewett[[2]],as.roman(i),after=(i-1))
}else {
x<- Jewett[[1]]
y<- Jewett[[2]]
}
Jewett[[1]] <- x
Jewett[[2]] <- as.roman(y)
}
return(Jewett)
}
find_amplitude <- function(arxeio){
file = xmlParse(arxeio)
xml_data <- xmlToList(file)
stimuli_side <- xml_data[["Waveform"]][[".attrs"]][["StimuliSide"]]
samples_num <- as.numeric(xml_data[["Waveform"]][["NumberOfSamples"]])
samples <- seq(0,14.98, by=0.03211991)
IPSI_A_Raw <- as.numeric(xml_data[["Waveform"]][["Response"]][["IPSI_A_Raw"]])
IPSI_B_Raw <- as.numeric(xml_data[["Waveform"]][["Response"]][["IPSI_B_Raw"]])
IPSI_Avg <- (IPSI_A_Raw + IPSI_B_Raw)/2
###### Metatropi Raw data se μV
gain <- strsplit(xml_data[["Waveform"]][["Gain"]]," ")[[1]]
gain <- as.numeric(gain[1])
gain <- dB_to_V(gain)
data_res <- (1.6/gain)/32768
amp <- (IPSI_Avg * data_res)*10^6
return(amp)
}
find_jewett <- function (x){
jewett <-c()
samples <- seq(0,14.98, by=0.03211991)
for (i in 1:length(x)) {
if ((0.8< samples[x[i]]) & (samples[x[i]] <= 2) ) {
jewett[i]<- 1
}
else if (2< samples[x[i]] & (samples[x[i]] <= 3)) {
jewett[i]<- 2
}
else if (3 < samples[x[i]] & (samples[x[i]] <= 4)) {
jewett[i]<- 3
}
else if (4 < samples[x[i]] & (samples[x[i]] <= 4.9)) {
jewett[i]<- 4
}
else if (4.9< samples[x[i]] & (samples[x[i]] <= 6.8)) {
jewett[i]<- 5
}
else {
jewett[i]<- 0
}
}
return (jewett)
}
find_peaks <- function (x, m = 3, error = 5){
shape <- diff(sign(diff(x, na.pad = FALSE)))
p <- sapply(which(shape < 0), FUN = function(i){
z <- i - m + 1
z <- ifelse(z > 0, z, 1)
w <- i + m + 1
w <- ifelse(w < length(x), w, length(x))
if(all(x[c(z : i, (i + 2) : w)] <= x[i + 1])) return(i + 1) else return(numeric(0))
})
p <- unlist(p)
# Ypologismos twn peaks me sfalma error
pks<-c()
counter =1
for (j in 1:length(p)) {
if(j == length(p)){
pks[counter]<- p[j]
}
else if (abs(p[j+1] - p[j]) > error){
pks[counter] <- p[j]
counter = counter + 1
}
else {
pks[counter] <- p[j]
next
counter = counter + 1
}
}
pks
}
fix_jewett <- function(peaks1,jewett1,peaks2,jewett2) {
norm <- c(1.5,2.5,3.5,4.5,5.5)
samples <- seq(0,14.98, by=0.03211991)
################## & (length(is.na(jewett2[duplicated(jewett2)]))<1)
if (length(grep(TRUE,duplicated(jewett2,incomparables=NA)))>0 & length(grep(TRUE,duplicated(jewett1,incomparables=NA)))==0 ){
duplicate <- jewett2[duplicated(jewett2,incomparables=NA)]
dup_num <- as.numeric((jewett2[duplicated(jewett2,incomparables=NA)]))
dup_pos<-as.numeric(grep(TRUE,duplicated(jewett2,incomparables=NA)))
test_peak <- c(peaks2[dup_pos-1],peaks2[dup_pos])
if (grep(duplicate,jewett1)) {
corr_pos<-grep(duplicate,jewett1)
corr_peak <- peaks1[corr_pos]
if (abs(corr_peak-test_peak[1]) <= abs(corr_peak-test_peak[2])) {
peaks2 <- peaks2[!peaks2 %in% test_peak[2]]
jewett2 <- jewett2[-dup_pos]
}else{
peaks2 <- peaks2[!peaks2 %in% test_peak[1]]
jewett2 <- jewett2[-(dup_pos-1)]
}
}else {
if (abs(samples[test_peak[1]]-norm[dup_num]) < abs(samples[test_peak[2]]-norm[dup_num])){
peaks2 <- peaks2[!peaks2 %in% test_peak[2]]
jewett2 <- jewett2[-dup_pos]
}else{
peaks2 <- peaks2[!peaks2 %in% test_peak[1]]
jewett2 <- jewett2[-(dup_pos-1)]
}
}
} #1o else if
##################
else if (length(grep(TRUE,duplicated(jewett1,incomparables=NA)))>0 & length(grep(TRUE,duplicated(jewett2,incomparables=NA)))==0 ){
duplicate <- jewett1[duplicated(jewett1,incomparables=NA)]
dup_num <- as.numeric(jewett1[duplicated(jewett1,incomparables=NA)])
dup_pos<-as.numeric(grep(TRUE,duplicated(jewett1,incomparables=NA)))
test_peak <- c(peaks1[dup_pos-1],peaks1[dup_pos])
if (grep(duplicate,jewett2)) {
corr_pos<-grep(duplicate,jewett2)
corr_peak <- peaks2[corr_pos]
if(abs(corr_peak-test_peak[1]) <= abs(corr_peak-test_peak[2])) {
peaks1 <- peaks1[!peaks1 %in% test_peak[2]]
jewett1 <- jewett1[-dup_pos]
}else {
peaks1 <- peaks1[!peaks1 %in% test_peak[1]]
jewett1 <- jewett1[-(dup_pos-1)]
}
}else{
if (abs(samples[test_peak[1]]-norm[dup_num]) < abs(samples[test_peak[2]]-norm[dup_num])){
peaks1 <- peaks1[!peaks1 %in% test_peak[2]]
jewett1 <- jewett1[-dup_pos]
}else{
peaks1 <- peaks1[!peaks2 %in% test_peak[1]]
jewett1 <- jewett1[-(dup_pos-1)]
}
}
}# eksw else if
##################
else if (length(grep(TRUE,duplicated(jewett1,incomparables=NA)))>0 & length(grep(TRUE,duplicated(jewett2,incomparables=NA)))>0){
j1<-fix_jewett_single(peaks1,jewett1)
j2<-fix_jewett_single(peaks2,jewett2)
peaks1 <- j1[[1]]
peaks2 <- j2[[1]]
jewett1 <- j1[[2]]
jewett2 <- j2[[2]]
}
else{
peaks1<-peaks1
peaks2<-peaks2
jewett1<-jewett1
jewett2<-jewett2
}# teleutaio else
Jewetts <- list(peaks1,peaks2,jewett1,jewett2)
return(Jewetts)
}
fix_jewett_single <- function(peaks,jewett,avg) {
norm <- c(1.5,2.5,3.5,4.5,5.5)
samples <- seq(0,14.98, by=0.03211991)
if (length(grep(TRUE,duplicated(jewett,incomparables=NA)))>0) {
duplicate <- jewett[duplicated(jewett,incomparables=NA)]
dup_num <- as.numeric((jewett[duplicated(jewett,incomparables=NA)]))
dup_pos<-as.numeric(grep(TRUE,duplicated(jewett,incomparables=NA)))
test_peak <- c(peaks[dup_pos-1],peaks[dup_pos])
if ((abs(samples[test_peak[1]]-norm[dup_num]) < abs(samples[test_peak[2]]-norm[dup_num])) & (avg[test_peak[1]] > avg[test_peak[2]])) {
# afairw teleutaio
peaks <- peaks[!peaks %in% test_peak[2]]
jewett <- jewett[-dup_pos]
}
else if ((abs(samples[test_peak[1]]-norm[dup_num]) > abs(samples[test_peak[2]]-norm[dup_num])) & (avg[test_peak[1]] < avg[test_peak[2]])) {
peaks <- peaks[!peaks %in% test_peak[1]]
jewett <- jewett[-(dup_pos-1)]
}
else {
if (avg[test_peak[1]] > avg[test_peak[2]]) {
peaks <- peaks[!peaks %in% test_peak[2]]
jewett <- jewett[-dup_pos]
} else {
# afairw proteleutaio
peaks <- peaks[!peaks %in% test_peak[1]]
jewett <- jewett[-(dup_pos-1)]
}
}
}
else {
peaks<-peaks
jewett<-jewett
}
Jewetts <- list(peaks,jewett)
# Deletes NA and corresponding peaks from Jewett List
while(length(grep(TRUE,is.na(Jewetts[[2]])))>=1) {
na <- grep(TRUE,is.na(Jewetts[[2]]))
for (n in na){
Jewetts[[1]]<-Jewetts[[1]][-n]
Jewetts[[2]]<-Jewetts[[2]][-n]
}
}
# Deletes remaining peaks if there is a diff
if (length(Jewetts[[1]]) > length(Jewetts[[2]])) {
diff <- length(Jewetts[[1]]) - length(Jewetts[[2]])
for (d in diff) {
Jewetts[[1]] <- Jewetts[[1]][-(length(Jewetts[[1]]))]
}
}
# Fills 5 peaks with NA
for (i in 1:5) {
if (is.na(match(as.roman(i),Jewetts[[2]]))) {
x<-append(Jewetts[[1]],NA,after=(i-1))
y<-append(Jewetts[[2]],as.roman(i),after=(i-1))
}else {
x<- Jewetts[[1]]
y<- Jewetts[[2]]
}
Jewetts[[1]] <- x
Jewetts[[2]] <- as.roman(y)
}
return(Jewetts)
}
get_age <-function(arxeio2) {
file2 = xmlParse(arxeio2)
xml_data2 <- xmlToList(file2)
year<-strsplit(xml_data2[["Waveform"]][[".attrs"]][["DayOfBirth"]],"-")[[1]]
year<-as.numeric(year[1])
age<-2020-year
return(age)
}
get_intensity <- function(arxeio) {
file = xmlParse(arxeio)
xml_data <- xmlToList(file)
intensity <- xml_data[["Waveform"]][[".attrs"]][["Intensity"]]
return(intensity)
}
get_ms <- function(peak) {
samples <- seq(0,14.98, by=0.03211991)
ms <- samples[peak]
return(ms)
}
get_stim_side<- function(arxeio) {
file = xmlParse(arxeio)
xml_data <- xmlToList(file)
stimuli_side <- xml_data[["Waveform"]][[".attrs"]][["StimuliSide"]]
return(stimuli_side)
}
LIF_I <- function(intensity,latency_I) {
if (is.na(latency_I)) {
y1_min<-c(0.67,0.19,0.45,1.09)
y1_max<-c(2.05,1.96,1.927,2)
x<-c(70,80,90,100)
plot(1, type="n", xlim=c(60,110), ylim=c(0,10),
main = 'Wave I Latency / Intensity Function',
panel.first = grid(),
xlab="Intensity (dB nHL)", ylab="Time (ms)")
polygon(c(x,x[4:1]),c(y1_min,y1_max[4:1]),col='grey')
legend(x='topright',legend='Wave I peak is not present')
lat_I_plot<-recordPlot()
}else {
y1_min<-c(0.67,0.19,0.45,1.09)
y1_max<-c(2.05,1.96,1.927,2)
x<-c(70,80,90,100)
plot(1, type="n", xlim=c(60,110), ylim=c(0,10),
main = 'Wave I Latency / Intensity Function',
panel.first = grid(),
xlab="Intensity (dB nHL)", ylab="Time (ms)")
polygon(c(x,x[4:1]),c(y1_min,y1_max[4:1]),col='grey')
points(intensity,latency_I,pch=21,col='black',bg='red')
lat_I_plot<-recordPlot()
}
return(lat_I_plot)
}
LIF_III <- function(intensity,latency_III) {
if(is.na(latency_III)){
y3_min<-c(3.24,2.38,2.6,3.24)
y3_max<-c(3.8,3.95,3.9,3.8)
x<-c(70,80,90,100)
plot(1, type="n", xlim=c(60,110), ylim=c(0,10),
main = 'Wave III Latency / Intensity Function',
panel.first = grid(),
xlab="Intensity (dB nHL)", ylab="Time (ms)")
polygon(c(x,x[4:1]),c(y3_min,y3_max[4:1]),col='grey')
legend(x='topright',legend='Wave III peak is not present')
lat_III_plot<-recordPlot()
}else{
y3_min<-c(3.24,2.38,2.6,3.24)
y3_max<-c(3.8,3.95,3.9,3.8)
x<-c(70,80,90,100)
plot(1, type="n", xlim=c(60,110), ylim=c(0,10),
main = 'Wave III Latency / Intensity Function',
panel.first = grid(),
xlab="Intensity (dB nHL)", ylab="Time (ms)")
polygon(c(x,x[4:1]),c(y3_min,y3_max[4:1]),col='grey')
points(intensity,latency_III,pch=21,col='black',bg='red')
lat_III_plot<-recordPlot()
}
return(lat_III_plot)
}
LIF_V <- function(intensity,latency_V) {
if(is.na(latency_V)){
y5_min<-c(4.3,4.14,4.5,4.6)
y5_max<-c(6.2,6.8,6.97,7)
x<-c(70,80,90,100)
plot(1, type="n", xlim=c(60,110), ylim=c(0,10),
main = 'Wave V Latency / Intensity Function',
panel.first = grid(),
xlab="Intensity (dB nHL)", ylab="Time (ms)")
polygon(c(x,x[4:1]),c(y5_min,y5_max[4:1]),col='grey')
legend(x='topright',legend='Wave III peak is not present')
lat_V_plot<-recordPlot()
}else{
y5_min<-c(4.3,4.14,4.5,4.6)
y5_max<-c(6.2,6.8,6.97,7)
x<-c(70,80,90,100)
plot(1, type="n", xlim=c(60,110), ylim=c(0,10),
main = 'Wave V Latency / Intensity Function',
panel.first = grid(),
xlab="Intensity (dB nHL)", ylab="Time (ms)")
polygon(c(x,x[4:1]),c(y5_min,y5_max[4:1]),col='grey')
points(intensity,latency_V,pch=21,col='black',bg='red')
lat_V_plot<-recordPlot()
}
return(lat_V_plot)
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.