In wadpac/mechanicalshakerexperiments: Mechanicalshakerexperiments
knitr::opts_chunk$set(echo = TRUE)
rm(list=ls())
graphics.off()
#==============================================================================
# Specify file paths and load data
shaker_experiments_folder = "/Users/annelindelettink/Documents/Work MacBook Pro Annelinde/Mechanical Shaker Machine"
datadir = paste0(shaker_experiments_folder, "/subsetted_raw_data/")
if (!dir.exists(datadir)) {
stop(paste0("Directory does not exist: ", datadir))
}
# Load mechanical shaker data for experiment 4
filename_experiment = paste0(datadir, "E4_hrmr.RData")
if (!file.exists(filename_experiment)) {
stop(paste0("File does not exist: ", filename_experiment))
}
load(filename_experiment)
# Load functions
#install.packages("mechanicalshakerexperiments")
#library(mechanicalshakerexperiments)
# source functions directly from file, to be replaced by package installation:
my_functions_folder = "/Users/annelindelettink/Documents/Work MacBook Pro Annelinde/Mechanical Shaker Machine/mechanicalshakerexperiments/R"
#my_functions_folder = "/home/vincent/projects/mechanicalshakerexperiments/R"
for (function_file in dir(my_functions_folder, full.names = T)) source(function_file) #load functions
# Select data with low sampling rate (100 Hz) only
index_high <- which(data$specifications$experiment == "ms_hrmr" & data$specifications$sampling_rate == "100")
data_high <- list(data = data$data[index_high],
specifications = data$specifications[index_high,])
rm(data)
outputdir <- "/Users/annelindelettink/Documents/Work MacBook Pro Annelinde/Mechanical Shaker Machine/analyses/E4_hrmr/"
Plot selected axis shaking direction
plot_list <- list()
for(signal in 1:length(data_high$data)) {
file_name <- paste0("acceleration data: ", names(data_high$data)[signal])
jpeg(paste0(paste0(outputdir, "plots"), paste0(file_name, ".jpeg")), width=600, height=500, res=120)
p <- ggplot2::ggplot() + ggplot2::geom_point(data = data_high$data[[signal]], ggplot2::aes(x = time, y = SD,) , alpha = 0.2, size = 0.25) + ggplot2::theme_bw() +
ggplot2::labs(title = file_name, y = expression(paste("acceleration (" ,italic("g"), ")"))) +
ggplot2::scale_color_manual(values = colors)
p
dev.off()
plot_list[[signal]] <- p
}
# Save and print the plot for each signal
for(plot in 1:length(plot_list)){
print(plot_list[[plot]])
}
rm(plot, plot_list, p)
Frequency domain analyses
# Plot parameters
plot = TRUE
xlabel = "Frequency (Hertz)"
ylabel = expression(paste("Spectral Density (" ,italic("g"), "2/Hertz)"))
XLIM = c(0, 5) # Limit frequency content to 5 Hz (as 250 rpm / 60 = 4.1667 Hz is the expected max)
# DERIVE RAW FREQUENCY SPECTRA (for visual inspection and save data for further analysis)
raw.spectra_high <- list()
for (file in 1:length(data_high$data)) { # for all accelerometer files
sampling_rate <- as.numeric(data_high$specifications$sampling_rate[file])
SD <- data_high$data[[file]]$SD #axis in shaking direction
rawSpectrum <- deriveSpectrum(SD, sampling_rate, file_name = names(data_high$data)[[file]], outputdir)
raw.spectra_high[[file]] <- rawSpectrum
}
rm(rawSpectrum, file)
names(raw.spectra_high) <- data_high$specifications$label
mean.freqspec_high <- deriveMeanSpectrum(raw.spectra_high)
raw_spectra <- list(spectra = raw.spectra_high, mean_spectrum = mean.freqspec_high, specifications = data_high$specifications) # Save derived raw spectra
save(raw_spectra, file = paste0(outputdir, "E4_frequency_spectra.RData"))
rm(raw.spectra_high, SD, sampling_rate)
#load(paste0(outputdir, "E4_frequency_spectra.RData"))
Mean frequency spectra
plot(raw_spectra$mean_spectrum$x, raw_spectra$mean_spectrum$y, xlab = xlabel, ylab = ylabel, type = "l", xlim = XLIM, bty= "l", main = "Mean frequency spectrum high sampling rate and mixed dynamic range")
rm(mean.freqspec_high)
Derive frequency bins
peaks_high <- derivePeaks(raw_spectra$mean_spectrum, 100, XLIM, plot = TRUE)
bins_high <- deriveBins(peaks_high, raw_spectra$mean_spectrum, XLIM = c(0, 5), plot = TRUE)
Frequency spectra plots
Overlay all frequency spectra
plot <- plot(raw_spectra$spectra[[1]]$specx, raw_spectra$spectra[[1]]$specy, main = "Overlay all frequency spectra",
xlab = xlabel, ylab = ylabel, type = "l", xlim = XLIM, bty= "l")
for(i in 1:length(bins_high)){
plot <- plot + abline(v=bins_high[i], col = "red")
}
for (spec in 2:length(raw_spectra$spectra)) {
plot + lines(raw_spectra$spectra[[spec]]$specx, raw_spectra$spectra[[spec]]$specy, type = "l")
}
rm(spec, i)
Seperate plots for frequency spectra
for (spec in 1:length(raw_spectra$spectra)) {
spectrum <- raw_spectra$spectra[[spec]]
title <- paste("Frequency spectrum plot: ", raw_spectra$specifications$label[[spec]])
plot(spectrum$specx, spectrum$specy, main = title,
xlab = xlabel, ylab = ylabel, type = "l", xlim = XLIM, bty= "l")
for(i in 1:length(bins_high)){
abline(v=bins_high[i], col = "red")
}
}
rm(spec, title, spectrum, i)
comparison_values <- list()
for(spectrum in 1:length(raw_spectra$spectra)){
comparison_values[[spectrum]] <- deriveComparisonValues(raw_spectra$spectra[[spectrum]], bins_high)
}
names(comparison_values) <- names(raw_spectra$spectra)
#cat("\nSaving data...")
save(comparison_values, file = paste0(outputdir, "E4_comparison_values.RData"))
rm(spectrum)
# Wide format: for both outcomes domFreq and meanPSD
df_wide <- data.frame()
id <- names(comparison_values)
df_wide <- cbind(id, data_high$specifications$brand, data_high$specifications$sampling_rate,
data_high$specifications$dynamic_range)
domFreq_wide <- data.frame()
meanPSD_wide <- data.frame()
for (accelerometer in 1:length(comparison_values)) {
domFreq_wide <- rbind(domFreq_wide, comparison_values[[accelerometer]]$domFreq)
meanPSD_wide <- rbind(meanPSD_wide, comparison_values[[accelerometer]]$meanPSD)
}
domFreq_ms_wide <- cbind(df_wide, domFreq_wide)
colnames(domFreq_ms_wide) <- c("serialnumber", "brand", "sampling_rate", "dynamic_range",
"1", "2", "3", "4", "5", "6", "7", "8", "9", "10", "11", "12", "13", "14")
meanPSD_ms_wide <- cbind(df_wide, meanPSD_wide)
colnames(meanPSD_ms_wide) <- c("serialnumber", "brand", "sampling_rate", "dynamic_range",
"1", "2", "3", "4", "5", "6", "7", "8", "9", "10", "11", "12", "13", "14")
rm(df_wide, id, domFreq_wide, meanPSD_wide, accelerometer)
# Combining domFreq and meanPSD (gather columns for outcome variables)
library(tidyverse)
domFreq_long <- domFreq_ms_wide %>%
gather(key = "bin", value = "domFreq", "1", "2", "3", "4", "5", "6", "7", "8", "9", "10", "11", "12", "13", "14") %>%
rstatix::convert_as_factor(serialnumber, bin, brand, sampling_rate, dynamic_range)
meanPSD_long <- meanPSD_ms_wide %>%
gather(key = "bin", value = "meanPSD", "1", "2", "3", "4", "5", "6", "7", "8", "9", "10", "11", "12", "13", "14") %>%
rstatix::convert_as_factor(serialnumber, bin, brand, sampling_rate, dynamic_range)
df_long <- merge(domFreq_long, meanPSD_long) # Join both data.frame objects to derive a single dataset
rm(domFreq_long, domFreq_ms_wide, meanPSD_long, meanPSD_ms_wide)
df_long <- df_long[order(df_long$serialnumber),]
#cat("\nSaving data...")
save(df_long, file = paste0(outputdir, "E4_comparison_values_long.RData"))
rm(peaks_high, index_high)
#load(outputdir, "E4_comparison_values_long.RData")
Statistical comparison: Mixed model analyses
Mean power spectral density: meanPSD ~ dynamic_range + bin + dynamic_range * bin
meanPSD_interaction_id <- lme4::lmer(meanPSD ~ bin*dynamic_range + (1 | serialnumber), REML = TRUE, data = df_long)
# Evaluate if adding a random intercept for frequency bin is necessary
meanPSD_interaction_id_bin <- lme4::lmer(meanPSD ~ bin*dynamic_range + (1 | serialnumber) + (1 | bin), REML = TRUE, data = df_long) # No convergence, so use the simpler model
anova(meanPSD_interaction_id, meanPSD_interaction_id_bin, refit = FALSE) # This is not necessary, so use simplest model
rm(meanPSD_interaction_id_bin)
Model assumption: normal distribution of residuals
ggpubr::ggqqplot(df_long, "meanPSD", facet.by = "bin", title = "QQ-plot mean power spectral density")
residuals_meanPSD <- resid(meanPSD_interaction_id)
#summary(residuals_meanPSD)
hist(residuals_meanPSD)
rm(residuals_meanPSD)
Model summary
car::Anova(meanPSD_interaction_id, type="II", test.statistic="F")
require(lmerTest)
test.meanPSD <- as(meanPSD_interaction_id, "merModLmerTest")
print(summary(test.meanPSD, ddf="Kenward-Roger"), correlation = FALSE)
#summary(meanPSD_interaction_id)
vcov(meanPSD_interaction_id, type = "random")
rm(test.meanPSD)
#report::report(meanPSD_interaction_id)
Pairwise comparisons
Bin
maineffect_bin <- emmeans::emmeans(meanPSD_interaction_id, comparison="pairwise",
specs= ~ bin, adjust="sidak", data = df_long)
maineffect_bin
emmeans::contrast(maineffect_bin, "pairwise")
rm(maineffect_bin)
Dominant frequency: domFreq ~ dynamic_range + bin + dynamic_range*bin
domFreq_interaction_id <- lme4::lmer(domFreq ~ bin*dynamic_range + (1 | serialnumber), REML = TRUE, data = df_long)
# Evaluate if adding a random intercept for frequency bin is necessary
domFreq_interaction_id__bin <- lme4::lmer(domFreq ~ bin*dynamic_range + (1 | serialnumber) + (1 | bin), REML = TRUE, data = df_long)
anova(domFreq_interaction_id, domFreq_interaction_id__bin, refit = FALSE) # This is not necessary, so use simplest model
rm(domFreq_interaction_id__bin)
Model assumption: normal distribution of residuals
ggpubr::ggqqplot(df_long, "domFreq", facet.by = "bin", title = "QQ-plot dominant frequency")
residuals_domFreq <- resid(domFreq_interaction_id)
#summary(residuals_domFreq)
hist(residuals_domFreq) # approximately normally distributed
rm(residuals_domFreq)
Model summary
car::Anova(domFreq_interaction_id, type="II", test.statistic="F")
require(lmerTest)
test.domFreq <- as(domFreq_interaction_id, "merModLmerTest")
print(summary(test.domFreq, ddf="Kenward-Roger"), correlation = FALSE)
#summary(domFreq_interaction_id)
vcov(domFreq_interaction_id, type = "random")
#report::report(domFreq_interaction_id)
rm(test.domFreq)
Pairwise comparisons
Bin
maineffect_bin <- emmeans::emmeans(domFreq_interaction_id, comparison="pairwise",
specs= ~ bin, adjust="sidak", data = df_long)
maineffect_bin
emmeans::contrast(maineffect_bin, "pairwise")
rm(maineffect_bin)
Boxplot
bxp_psd <- createBoxplot(df_long, bins_high, outcome = "meanPSD", group = "dynamic_range")
bxp_domFreq <- createBoxplot(df_long, bins_high, outcome = "domFreq", group = "dynamic_range")
ggplot2::ggsave(file=paste0(outputdir, "plots/boxplot_E4_meanPSD.png"), bxp_psd, width = 10, height = 8, dpi = 900) #saves g
ggplot2::ggsave(file=paste0(outputdir, "plots/boxplot_E4_domFreq.png"), bxp_domFreq, width = 10, height = 8, dpi = 900) #saves g
plot(bxp_psd)
plot(bxp_domFreq)
Time domain analyses
Cross-correlation coefficients
correlation_matrices <- computeCrossCorrelations(data_high$data, plot = FALSE)
save(correlation_matrices, file = paste0(outputdir, "E4_cross_correlations_overall.RData"))
correlations_perShakerCondition <- computeCrossCorrelationsPerShakerCondition(data_high$data)
save(correlations_perShakerCondition, file = paste0(outputdir, "E4_cross_correlations_perShakerCondition.RData"))
#load(paste0(outputdir, "E4_cross_correlations_perShakerCondition.RData"))
Overall heatmap
fullmatrix <- correlation_matrices$correlationMatrix
# Add the dynamic range to the labels
labels <- rownames(fullmatrix)
dr_labels <- c()
for(lab in 1:length(labels)){
dr_labels <- c(dr_labels, paste(data_high$specifications[data_high$specifications$label == labels[lab], 6], labels[lab], sep = "_"))
}
colnames(fullmatrix) <- dr_labels
rownames(fullmatrix) <- dr_labels
heatmap <- corrplot::corrplot(pmax(fullmatrix, t(fullmatrix), na.rm = TRUE),
type = 'lower', method = 'color', order = "alphabet",
number.cex = 1.25, cl.cex = 2, tl.cex = 2,
col.lim = c(0,1), addCoef.col = 'black', tl.srt = 45,
diag = TRUE, tl.col = "black")
rm(fullmatrix, heatmap)
Heatmap per shaker condition
fullmatrices <- correlations_perShakerCondition
# Add the dynamic range to the labels
for(matrix in 1:length(fullmatrices)){
colnames(fullmatrices[[matrix]]$correlationMatrix) <- dr_labels
rownames(fullmatrices[[matrix]]$correlationMatrix) <- dr_labels
}
heatmaps <- list()
for (condition in 1:length(fullmatrices)){
png(height=1800, width=1800, file=paste0(outputdir, paste0("heatmaps/heatmap_E4_", condition, ".png")))
heatmaps[[condition]] <- corrplot::corrplot(pmax(fullmatrices[[condition]]$correlationMatrix,
t(fullmatrices[[condition]]$correlationMatrix), na.rm = TRUE),
title = paste0("shaker condition: ", names(fullmatrices)[condition]),
type = 'lower', method = 'color', order = "alphabet",
number.cex = 1.25, cl.cex = 2, tl.cex = 2,
col.lim = c(-1,1), addCoef.col = 'black', tl.srt = 45, mar=c(0,0,2,0),
diag = TRUE, tl.col = "black")
dev.off()
}
rm(fullmatrices, heatmaps)
Comparison of cross-correlation coefficients between dynamic ranges
avg_crosscor_between <- list()
avg_crosscor_within <- list()
for (condition in 1:length(correlations_perShakerCondition)){
avg_crosscor <- dataframeShapeComparisonBeween(correlations_perShakerCondition[[condition]], data_high$specifications, between = "dynamic_range")
avg_crosscor_between[[condition]] <- avg_crosscor[avg_crosscor$same_group == "0",]
avg_crosscor_within[[condition]] <- avg_crosscor[avg_crosscor$same_group == "1",]
}
rm(avg_crosscor, condition)
names(avg_crosscor_between) <- names(correlations_perShakerCondition)
names(avg_crosscor_within) <- names(correlations_perShakerCondition)
table_between <- list()
for(condition in 1:length(avg_crosscor_between)){
df_correlations <- avg_crosscor_between[[condition]]
df_correlations$correlation <- as.numeric(df_correlations$correlation)
cor24 <- rbind(df_correlations[which(df_correlations$between1 == 2 & df_correlations$between2 == 4) ,],
df_correlations[which(df_correlations$between1 == 4 & df_correlations$between2 == 2) ,])
cor28 <- rbind(df_correlations[which(df_correlations$between1 == 2 & df_correlations$between2 == 8) ,],
df_correlations[which(df_correlations$between1 == 8 & df_correlations$between2 == 2) ,])
cor216 <- rbind(df_correlations[which(df_correlations$between1 == 2 & df_correlations$between2 == 16) ,],
df_correlations[which(df_correlations$between1 == 16 & df_correlations$between2 == 2) ,])
cor48 <- rbind(df_correlations[which(df_correlations$between1 == 4 & df_correlations$between2 == 8) ,],
df_correlations[which(df_correlations$between1 == 8 & df_correlations$between2 == 4) ,])
cor416 <- rbind(df_correlations[which(df_correlations$between1 == 4 & df_correlations$between2 == 16) ,],
df_correlations[which(df_correlations$between1 == 16 & df_correlations$between2 == 4) ,])
cor816 <- rbind(df_correlations[which(df_correlations$between1 == 8 & df_correlations$between2 == 16) ,],
df_correlations[which(df_correlations$between1 == 16 & df_correlations$between2 == 8) ,])
tab <- rbind(summary(cor24$correlation), summary(cor28$correlation), summary(cor216$correlation),
summary(cor48$correlation), summary(cor416$correlation), summary(cor816$correlation))
sd <- c(sd(cor24$correlation), sd(cor28$correlation), sd(cor216$correlation),
sd(cor48$correlation), sd(cor416$correlation), sd(cor816$correlation))
tab <- cbind(tab, sd)
groups <- c("2_vs_4", "2_vs_8", "2_vs_18", "4_vs_8", "4_vs_16", "8_vs_16")
tab <- tab[,-1]
tab <- tab[,-5]
tab <- cbind(groups, tab)
table_between[[condition]] <- tab
}
rm(tab, cor24, cor28, cor216, cor48, cor416, cor816, sd, groups, condition)
names(table_between) <- names(correlations_perShakerCondition)
save(table_between, file = paste0(outputdir, "E4_line_graph_data.RData"))
for(condition in 1:length(table_between)){
print(paste0("Shaker condition: ", names(table_between)[condition]))
print(table_between[[condition]])
}
shaker_frequencies <- c(30, 40, 50, 62, 75, 87, 100, 112, 125, 137, 150, 160, 175, 187, 200, 212, 225, 237, 250)
groups <- c("2_vs_4", "2_vs_8", "2_vs_18", "4_vs_8", "4_vs_16", "8_vs_16")
comparison <- rep(groups, length(shaker_frequencies))
correlations <- c()
se <- c()
shaker_frequency <- c()
for(condition in 1:length(table_between)){
correlations <- c(correlations, table_between[[condition]][,4])
se <- c(se, as.numeric(table_between[[condition]][,6])*1.96)
shaker_frequency <- c(shaker_frequency, rep(shaker_frequencies[condition], length(table_between[[condition]][,4])))
}
overview_tab_between <- cbind(comparison, correlations, shaker_frequency, se)
overview_tab_between <- as.data.frame(overview_tab_between)
line_plot <- ggplot2::ggplot(data = overview_tab_between, ggplot2::aes(x = as.factor(shaker_frequency), y = as.numeric(correlations), group = as.factor(comparison))) +
ggplot2::geom_line(ggplot2::aes(color = as.factor(comparison), linetype = as.factor(comparison)))+
ggplot2::geom_point(ggplot2::aes(color = as.factor(comparison), shape = as.factor(comparison))) +
ggplot2::labs(x="Shaker frequency (rpm)", y = "Cross-correlation coefficient", color = "Dynamic range comparison") +
ggplot2::scale_x_discrete(limits = c("30", "40", "50", "62", "75", "87", "100", "112", "125", "137", "150", "160", "175", "187", "200", "212", "225", "237", "250")) +
ggplot2::theme_classic() +
ggplot2::scale_color_manual(values = viridis::viridis(6)) +
ggplot2::theme(legend.position="bottom") +
ggplot2::geom_hline(yintercept=0.25, linetype="dashed") +
ggplot2::geom_hline(yintercept=0.5, linetype="dashed") +
ggplot2::geom_hline(yintercept=0.75, linetype="dashed") +
ggplot2::geom_errorbar(ggplot2::aes(ymin = as.numeric(correlations) - as.numeric(se),
ymax = as.numeric(correlations) + as.numeric(se),
color = as.factor(comparison)),
width=0.1) +
ggplot2::labs(color = "Dynamic range comparison", linetype = "Dynamic range comparison", shape = "Dynamic range comparison")
ggplot2::ggsave(file=paste0(outputdir, "plots/linegraph_E4_crosscorrelations.png"), line_plot, width = 10, height = 8, dpi = 900) #saves g
plot(line_plot)
wadpac/mechanicalshakerexperiments documentation built on July 2, 2024, 11:49 p.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.
knitr::opts_chunk$set(echo = TRUE) rm(list=ls()) graphics.off() #============================================================================== # Specify file paths and load data shaker_experiments_folder = "/Users/annelindelettink/Documents/Work MacBook Pro Annelinde/Mechanical Shaker Machine" datadir = paste0(shaker_experiments_folder, "/subsetted_raw_data/") if (!dir.exists(datadir)) { stop(paste0("Directory does not exist: ", datadir)) } # Load mechanical shaker data for experiment 4 filename_experiment = paste0(datadir, "E4_hrmr.RData") if (!file.exists(filename_experiment)) { stop(paste0("File does not exist: ", filename_experiment)) } load(filename_experiment) # Load functions #install.packages("mechanicalshakerexperiments") #library(mechanicalshakerexperiments) # source functions directly from file, to be replaced by package installation: my_functions_folder = "/Users/annelindelettink/Documents/Work MacBook Pro Annelinde/Mechanical Shaker Machine/mechanicalshakerexperiments/R" #my_functions_folder = "/home/vincent/projects/mechanicalshakerexperiments/R" for (function_file in dir(my_functions_folder, full.names = T)) source(function_file) #load functions # Select data with low sampling rate (100 Hz) only index_high <- which(data$specifications$experiment == "ms_hrmr" & data$specifications$sampling_rate == "100") data_high <- list(data = data$data[index_high], specifications = data$specifications[index_high,]) rm(data) outputdir <- "/Users/annelindelettink/Documents/Work MacBook Pro Annelinde/Mechanical Shaker Machine/analyses/E4_hrmr/"
Plot selected axis shaking direction
plot_list <- list() for(signal in 1:length(data_high$data)) { file_name <- paste0("acceleration data: ", names(data_high$data)[signal]) jpeg(paste0(paste0(outputdir, "plots"), paste0(file_name, ".jpeg")), width=600, height=500, res=120) p <- ggplot2::ggplot() + ggplot2::geom_point(data = data_high$data[[signal]], ggplot2::aes(x = time, y = SD,) , alpha = 0.2, size = 0.25) + ggplot2::theme_bw() + ggplot2::labs(title = file_name, y = expression(paste("acceleration (" ,italic("g"), ")"))) + ggplot2::scale_color_manual(values = colors) p dev.off() plot_list[[signal]] <- p }
# Save and print the plot for each signal for(plot in 1:length(plot_list)){ print(plot_list[[plot]]) } rm(plot, plot_list, p)
Frequency domain analyses
# Plot parameters plot = TRUE xlabel = "Frequency (Hertz)" ylabel = expression(paste("Spectral Density (" ,italic("g"), "2/Hertz)")) XLIM = c(0, 5) # Limit frequency content to 5 Hz (as 250 rpm / 60 = 4.1667 Hz is the expected max)
# DERIVE RAW FREQUENCY SPECTRA (for visual inspection and save data for further analysis) raw.spectra_high <- list() for (file in 1:length(data_high$data)) { # for all accelerometer files sampling_rate <- as.numeric(data_high$specifications$sampling_rate[file]) SD <- data_high$data[[file]]$SD #axis in shaking direction rawSpectrum <- deriveSpectrum(SD, sampling_rate, file_name = names(data_high$data)[[file]], outputdir) raw.spectra_high[[file]] <- rawSpectrum } rm(rawSpectrum, file) names(raw.spectra_high) <- data_high$specifications$label mean.freqspec_high <- deriveMeanSpectrum(raw.spectra_high)
raw_spectra <- list(spectra = raw.spectra_high, mean_spectrum = mean.freqspec_high, specifications = data_high$specifications) # Save derived raw spectra save(raw_spectra, file = paste0(outputdir, "E4_frequency_spectra.RData")) rm(raw.spectra_high, SD, sampling_rate) #load(paste0(outputdir, "E4_frequency_spectra.RData"))
Mean frequency spectra
plot(raw_spectra$mean_spectrum$x, raw_spectra$mean_spectrum$y, xlab = xlabel, ylab = ylabel, type = "l", xlim = XLIM, bty= "l", main = "Mean frequency spectrum high sampling rate and mixed dynamic range") rm(mean.freqspec_high)
Derive frequency bins
peaks_high <- derivePeaks(raw_spectra$mean_spectrum, 100, XLIM, plot = TRUE) bins_high <- deriveBins(peaks_high, raw_spectra$mean_spectrum, XLIM = c(0, 5), plot = TRUE)
Frequency spectra plots
Overlay all frequency spectra
plot <- plot(raw_spectra$spectra[[1]]$specx, raw_spectra$spectra[[1]]$specy, main = "Overlay all frequency spectra", xlab = xlabel, ylab = ylabel, type = "l", xlim = XLIM, bty= "l") for(i in 1:length(bins_high)){ plot <- plot + abline(v=bins_high[i], col = "red") } for (spec in 2:length(raw_spectra$spectra)) { plot + lines(raw_spectra$spectra[[spec]]$specx, raw_spectra$spectra[[spec]]$specy, type = "l") } rm(spec, i)
Seperate plots for frequency spectra
for (spec in 1:length(raw_spectra$spectra)) { spectrum <- raw_spectra$spectra[[spec]] title <- paste("Frequency spectrum plot: ", raw_spectra$specifications$label[[spec]]) plot(spectrum$specx, spectrum$specy, main = title, xlab = xlabel, ylab = ylabel, type = "l", xlim = XLIM, bty= "l") for(i in 1:length(bins_high)){ abline(v=bins_high[i], col = "red") } } rm(spec, title, spectrum, i)
comparison_values <- list() for(spectrum in 1:length(raw_spectra$spectra)){ comparison_values[[spectrum]] <- deriveComparisonValues(raw_spectra$spectra[[spectrum]], bins_high) } names(comparison_values) <- names(raw_spectra$spectra) #cat("\nSaving data...") save(comparison_values, file = paste0(outputdir, "E4_comparison_values.RData")) rm(spectrum)
# Wide format: for both outcomes domFreq and meanPSD df_wide <- data.frame() id <- names(comparison_values) df_wide <- cbind(id, data_high$specifications$brand, data_high$specifications$sampling_rate, data_high$specifications$dynamic_range) domFreq_wide <- data.frame() meanPSD_wide <- data.frame() for (accelerometer in 1:length(comparison_values)) { domFreq_wide <- rbind(domFreq_wide, comparison_values[[accelerometer]]$domFreq) meanPSD_wide <- rbind(meanPSD_wide, comparison_values[[accelerometer]]$meanPSD) } domFreq_ms_wide <- cbind(df_wide, domFreq_wide) colnames(domFreq_ms_wide) <- c("serialnumber", "brand", "sampling_rate", "dynamic_range", "1", "2", "3", "4", "5", "6", "7", "8", "9", "10", "11", "12", "13", "14") meanPSD_ms_wide <- cbind(df_wide, meanPSD_wide) colnames(meanPSD_ms_wide) <- c("serialnumber", "brand", "sampling_rate", "dynamic_range", "1", "2", "3", "4", "5", "6", "7", "8", "9", "10", "11", "12", "13", "14") rm(df_wide, id, domFreq_wide, meanPSD_wide, accelerometer)
# Combining domFreq and meanPSD (gather columns for outcome variables) library(tidyverse) domFreq_long <- domFreq_ms_wide %>% gather(key = "bin", value = "domFreq", "1", "2", "3", "4", "5", "6", "7", "8", "9", "10", "11", "12", "13", "14") %>% rstatix::convert_as_factor(serialnumber, bin, brand, sampling_rate, dynamic_range) meanPSD_long <- meanPSD_ms_wide %>% gather(key = "bin", value = "meanPSD", "1", "2", "3", "4", "5", "6", "7", "8", "9", "10", "11", "12", "13", "14") %>% rstatix::convert_as_factor(serialnumber, bin, brand, sampling_rate, dynamic_range) df_long <- merge(domFreq_long, meanPSD_long) # Join both data.frame objects to derive a single dataset rm(domFreq_long, domFreq_ms_wide, meanPSD_long, meanPSD_ms_wide) df_long <- df_long[order(df_long$serialnumber),] #cat("\nSaving data...") save(df_long, file = paste0(outputdir, "E4_comparison_values_long.RData")) rm(peaks_high, index_high) #load(outputdir, "E4_comparison_values_long.RData")
Statistical comparison: Mixed model analyses
Mean power spectral density: meanPSD ~ dynamic_range + bin + dynamic_range * bin
meanPSD_interaction_id <- lme4::lmer(meanPSD ~ bin*dynamic_range + (1 | serialnumber), REML = TRUE, data = df_long) # Evaluate if adding a random intercept for frequency bin is necessary meanPSD_interaction_id_bin <- lme4::lmer(meanPSD ~ bin*dynamic_range + (1 | serialnumber) + (1 | bin), REML = TRUE, data = df_long) # No convergence, so use the simpler model anova(meanPSD_interaction_id, meanPSD_interaction_id_bin, refit = FALSE) # This is not necessary, so use simplest model rm(meanPSD_interaction_id_bin)
Model assumption: normal distribution of residuals
ggpubr::ggqqplot(df_long, "meanPSD", facet.by = "bin", title = "QQ-plot mean power spectral density") residuals_meanPSD <- resid(meanPSD_interaction_id) #summary(residuals_meanPSD) hist(residuals_meanPSD) rm(residuals_meanPSD)
Model summary
car::Anova(meanPSD_interaction_id, type="II", test.statistic="F") require(lmerTest) test.meanPSD <- as(meanPSD_interaction_id, "merModLmerTest") print(summary(test.meanPSD, ddf="Kenward-Roger"), correlation = FALSE) #summary(meanPSD_interaction_id) vcov(meanPSD_interaction_id, type = "random") rm(test.meanPSD) #report::report(meanPSD_interaction_id)
Pairwise comparisons
Bin
maineffect_bin <- emmeans::emmeans(meanPSD_interaction_id, comparison="pairwise", specs= ~ bin, adjust="sidak", data = df_long) maineffect_bin emmeans::contrast(maineffect_bin, "pairwise") rm(maineffect_bin)
Dominant frequency: domFreq ~ dynamic_range + bin + dynamic_range*bin
domFreq_interaction_id <- lme4::lmer(domFreq ~ bin*dynamic_range + (1 | serialnumber), REML = TRUE, data = df_long) # Evaluate if adding a random intercept for frequency bin is necessary domFreq_interaction_id__bin <- lme4::lmer(domFreq ~ bin*dynamic_range + (1 | serialnumber) + (1 | bin), REML = TRUE, data = df_long) anova(domFreq_interaction_id, domFreq_interaction_id__bin, refit = FALSE) # This is not necessary, so use simplest model rm(domFreq_interaction_id__bin)
Model assumption: normal distribution of residuals
ggpubr::ggqqplot(df_long, "domFreq", facet.by = "bin", title = "QQ-plot dominant frequency") residuals_domFreq <- resid(domFreq_interaction_id) #summary(residuals_domFreq) hist(residuals_domFreq) # approximately normally distributed rm(residuals_domFreq)
Model summary
car::Anova(domFreq_interaction_id, type="II", test.statistic="F") require(lmerTest) test.domFreq <- as(domFreq_interaction_id, "merModLmerTest") print(summary(test.domFreq, ddf="Kenward-Roger"), correlation = FALSE) #summary(domFreq_interaction_id) vcov(domFreq_interaction_id, type = "random") #report::report(domFreq_interaction_id) rm(test.domFreq)
Pairwise comparisons
Bin
maineffect_bin <- emmeans::emmeans(domFreq_interaction_id, comparison="pairwise", specs= ~ bin, adjust="sidak", data = df_long) maineffect_bin emmeans::contrast(maineffect_bin, "pairwise") rm(maineffect_bin)
Boxplot
bxp_psd <- createBoxplot(df_long, bins_high, outcome = "meanPSD", group = "dynamic_range") bxp_domFreq <- createBoxplot(df_long, bins_high, outcome = "domFreq", group = "dynamic_range") ggplot2::ggsave(file=paste0(outputdir, "plots/boxplot_E4_meanPSD.png"), bxp_psd, width = 10, height = 8, dpi = 900) #saves g ggplot2::ggsave(file=paste0(outputdir, "plots/boxplot_E4_domFreq.png"), bxp_domFreq, width = 10, height = 8, dpi = 900) #saves g
plot(bxp_psd)
plot(bxp_domFreq)
Time domain analyses
Cross-correlation coefficients
correlation_matrices <- computeCrossCorrelations(data_high$data, plot = FALSE) save(correlation_matrices, file = paste0(outputdir, "E4_cross_correlations_overall.RData")) correlations_perShakerCondition <- computeCrossCorrelationsPerShakerCondition(data_high$data) save(correlations_perShakerCondition, file = paste0(outputdir, "E4_cross_correlations_perShakerCondition.RData")) #load(paste0(outputdir, "E4_cross_correlations_perShakerCondition.RData"))
Overall heatmap
fullmatrix <- correlation_matrices$correlationMatrix # Add the dynamic range to the labels labels <- rownames(fullmatrix) dr_labels <- c() for(lab in 1:length(labels)){ dr_labels <- c(dr_labels, paste(data_high$specifications[data_high$specifications$label == labels[lab], 6], labels[lab], sep = "_")) } colnames(fullmatrix) <- dr_labels rownames(fullmatrix) <- dr_labels heatmap <- corrplot::corrplot(pmax(fullmatrix, t(fullmatrix), na.rm = TRUE), type = 'lower', method = 'color', order = "alphabet", number.cex = 1.25, cl.cex = 2, tl.cex = 2, col.lim = c(0,1), addCoef.col = 'black', tl.srt = 45, diag = TRUE, tl.col = "black") rm(fullmatrix, heatmap)
Heatmap per shaker condition
fullmatrices <- correlations_perShakerCondition # Add the dynamic range to the labels for(matrix in 1:length(fullmatrices)){ colnames(fullmatrices[[matrix]]$correlationMatrix) <- dr_labels rownames(fullmatrices[[matrix]]$correlationMatrix) <- dr_labels } heatmaps <- list() for (condition in 1:length(fullmatrices)){ png(height=1800, width=1800, file=paste0(outputdir, paste0("heatmaps/heatmap_E4_", condition, ".png"))) heatmaps[[condition]] <- corrplot::corrplot(pmax(fullmatrices[[condition]]$correlationMatrix, t(fullmatrices[[condition]]$correlationMatrix), na.rm = TRUE), title = paste0("shaker condition: ", names(fullmatrices)[condition]), type = 'lower', method = 'color', order = "alphabet", number.cex = 1.25, cl.cex = 2, tl.cex = 2, col.lim = c(-1,1), addCoef.col = 'black', tl.srt = 45, mar=c(0,0,2,0), diag = TRUE, tl.col = "black") dev.off() } rm(fullmatrices, heatmaps)
Comparison of cross-correlation coefficients between dynamic ranges
avg_crosscor_between <- list() avg_crosscor_within <- list() for (condition in 1:length(correlations_perShakerCondition)){ avg_crosscor <- dataframeShapeComparisonBeween(correlations_perShakerCondition[[condition]], data_high$specifications, between = "dynamic_range") avg_crosscor_between[[condition]] <- avg_crosscor[avg_crosscor$same_group == "0",] avg_crosscor_within[[condition]] <- avg_crosscor[avg_crosscor$same_group == "1",] } rm(avg_crosscor, condition) names(avg_crosscor_between) <- names(correlations_perShakerCondition) names(avg_crosscor_within) <- names(correlations_perShakerCondition)
table_between <- list() for(condition in 1:length(avg_crosscor_between)){ df_correlations <- avg_crosscor_between[[condition]] df_correlations$correlation <- as.numeric(df_correlations$correlation) cor24 <- rbind(df_correlations[which(df_correlations$between1 == 2 & df_correlations$between2 == 4) ,], df_correlations[which(df_correlations$between1 == 4 & df_correlations$between2 == 2) ,]) cor28 <- rbind(df_correlations[which(df_correlations$between1 == 2 & df_correlations$between2 == 8) ,], df_correlations[which(df_correlations$between1 == 8 & df_correlations$between2 == 2) ,]) cor216 <- rbind(df_correlations[which(df_correlations$between1 == 2 & df_correlations$between2 == 16) ,], df_correlations[which(df_correlations$between1 == 16 & df_correlations$between2 == 2) ,]) cor48 <- rbind(df_correlations[which(df_correlations$between1 == 4 & df_correlations$between2 == 8) ,], df_correlations[which(df_correlations$between1 == 8 & df_correlations$between2 == 4) ,]) cor416 <- rbind(df_correlations[which(df_correlations$between1 == 4 & df_correlations$between2 == 16) ,], df_correlations[which(df_correlations$between1 == 16 & df_correlations$between2 == 4) ,]) cor816 <- rbind(df_correlations[which(df_correlations$between1 == 8 & df_correlations$between2 == 16) ,], df_correlations[which(df_correlations$between1 == 16 & df_correlations$between2 == 8) ,]) tab <- rbind(summary(cor24$correlation), summary(cor28$correlation), summary(cor216$correlation), summary(cor48$correlation), summary(cor416$correlation), summary(cor816$correlation)) sd <- c(sd(cor24$correlation), sd(cor28$correlation), sd(cor216$correlation), sd(cor48$correlation), sd(cor416$correlation), sd(cor816$correlation)) tab <- cbind(tab, sd) groups <- c("2_vs_4", "2_vs_8", "2_vs_18", "4_vs_8", "4_vs_16", "8_vs_16") tab <- tab[,-1] tab <- tab[,-5] tab <- cbind(groups, tab) table_between[[condition]] <- tab } rm(tab, cor24, cor28, cor216, cor48, cor416, cor816, sd, groups, condition) names(table_between) <- names(correlations_perShakerCondition) save(table_between, file = paste0(outputdir, "E4_line_graph_data.RData"))
for(condition in 1:length(table_between)){ print(paste0("Shaker condition: ", names(table_between)[condition])) print(table_between[[condition]]) }
shaker_frequencies <- c(30, 40, 50, 62, 75, 87, 100, 112, 125, 137, 150, 160, 175, 187, 200, 212, 225, 237, 250) groups <- c("2_vs_4", "2_vs_8", "2_vs_18", "4_vs_8", "4_vs_16", "8_vs_16") comparison <- rep(groups, length(shaker_frequencies)) correlations <- c() se <- c() shaker_frequency <- c() for(condition in 1:length(table_between)){ correlations <- c(correlations, table_between[[condition]][,4]) se <- c(se, as.numeric(table_between[[condition]][,6])*1.96) shaker_frequency <- c(shaker_frequency, rep(shaker_frequencies[condition], length(table_between[[condition]][,4]))) } overview_tab_between <- cbind(comparison, correlations, shaker_frequency, se) overview_tab_between <- as.data.frame(overview_tab_between) line_plot <- ggplot2::ggplot(data = overview_tab_between, ggplot2::aes(x = as.factor(shaker_frequency), y = as.numeric(correlations), group = as.factor(comparison))) + ggplot2::geom_line(ggplot2::aes(color = as.factor(comparison), linetype = as.factor(comparison)))+ ggplot2::geom_point(ggplot2::aes(color = as.factor(comparison), shape = as.factor(comparison))) + ggplot2::labs(x="Shaker frequency (rpm)", y = "Cross-correlation coefficient", color = "Dynamic range comparison") + ggplot2::scale_x_discrete(limits = c("30", "40", "50", "62", "75", "87", "100", "112", "125", "137", "150", "160", "175", "187", "200", "212", "225", "237", "250")) + ggplot2::theme_classic() + ggplot2::scale_color_manual(values = viridis::viridis(6)) + ggplot2::theme(legend.position="bottom") + ggplot2::geom_hline(yintercept=0.25, linetype="dashed") + ggplot2::geom_hline(yintercept=0.5, linetype="dashed") + ggplot2::geom_hline(yintercept=0.75, linetype="dashed") + ggplot2::geom_errorbar(ggplot2::aes(ymin = as.numeric(correlations) - as.numeric(se), ymax = as.numeric(correlations) + as.numeric(se), color = as.factor(comparison)), width=0.1) + ggplot2::labs(color = "Dynamic range comparison", linetype = "Dynamic range comparison", shape = "Dynamic range comparison") ggplot2::ggsave(file=paste0(outputdir, "plots/linegraph_E4_crosscorrelations.png"), line_plot, width = 10, height = 8, dpi = 900) #saves g
plot(line_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.