data_dir <- file.path(getwd(), "data")Suboptimal sensory reweighting in car sick individuals: a postural study
Environment initialization
Specifying the paths for the data and scripts directories :
The data directory hosts the .brute files generated by our posturography platform.
.brute files store the posturography data fo each trial performed by each participants
The first part of the name of each files (i.e. s046) corresponds to the identification number of each participant.
The second part of the name corresponds to the number of the trial. Respectively cd1, cd2 and cd3 for the first, second and third trials.
The first column of these files (temps
) is the timestamp of each measure (every 10 ms during 180s).
The second column Consigne
correspond to the angle the platform must output to maintain sway-reference.
The third column Angle
corresponds of the current angle of the platform in space in the anteroposterior axis.
The four following columns (Capteur (Kg) (Nord Ouest)
, Capteur (Kg) (Nord Est)
, Capteur (Kg) (Sud Est)
, Capteur (Kg) (Sud Ouest)
), account for the force sensed on each of the four cardinally placed force sensors in Kg.
From these forces, the last two columns were computed to represent the Barycenter in mm (respectively on the x and y axis) of the participant on the plate.
scripts_dir <- file.path(getwd(), "scripts")Source the load-packages script and load/install (if necessary) needed packages :
# Specify the packages needed for the project
pkgs <- c("tidyverse",
"psych",
"signal",
"psd",
"tools",
"viridis",
"rmarkdown",
"writexl")
# Load the packages and install them if necessary
lpkg(pkgs)tidyverse psych signal psd tools viridis rmarkdown writexl
TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE Data setup
Loading the questionnaire data and compute the car-sickness score.
subj_quest <- read_csv(file.path(data_dir, "subj_quest.csv"))subj_quest <-
subj_quest %>%
mutate(CarSick = replace_na(CarMoun +
CarStraigh +
CarFront +
CarBack +
CarRead, 0)) Setting-up signal parameters :
# Sampling frequency (in Hz)
samp_freq <- 100
# Duration of each phase (static/sway-ref/static, in sec)
phase_dur <- 60
# Minimum and maximum frequencies of the spectral analysis
min_freq <- 0
max_freq <- 5
# FFT window duration (in sec)
window_dur <- 10
# Decimation ratio (not used since reasonable compute time)
dratio <- 1
# Rotation function for future plots
rotate <- function(x) t(apply(x, 2, rev))
# Normalization function for the spectra
sp_norm <- function(x) {
x / max(x)
}
# Number of samples in the analysis
nb_samp <- floor(3 * phase_dur * samp_freq)
# Width of the FFT analysis window
window_size <- round(window_dur * samp_freq)
# Window profile
window <- hamming(window_size)
# Initial and final sample indices
initial_idx <- max(c(
1, floor(window_size * min_freq / samp_freq)))
final_idx <- floor(window_size * max_freq / samp_freq)
# Number of steps for FFT analysis
nb_step <- nb_samp - window_sizeInitialization of the spectra storage :
# Temporary matrix for storing the spectra
spectra <- matrix(
NA, nrow = nb_step/dratio, ncol = final_idx - initial_idx + 1)
# List for the storing of results
rspec <- list()
rsptat <- list()
datastat <- tibble(
subject = character(),
condition = character(),
er1 = numeric(),
er2 = numeric(),
er3 = numeric(),
er4 = numeric(),
er5 = numeric(),
er2_er1 = numeric(),
er4_er3 = numeric(),
carsick = numeric(),
age = numeric(),
sexe = character(),
height = numeric(),
weight = numeric(),
validity = numeric())
cond <- c("cd1", "cd2", "cd3")Extracting participants weights from the COP files :
weightlist <- list()
# Looping through the subjects brute data and extracting weights :
for (s in seq(1, nrow(subj_quest))) {
subject <- subj_quest$Id[s]
filename <- paste0(
file_path_sans_ext(as.character(subject)),
"-",
"cd1",
".brute")
weightdf <- read.table(file.path(data_dir, filename),
sep = "\t",
skip = 1, header = TRUE)
weight <- weightdf[1,4]+
weightdf[1,5]+
weightdf[1,6]+
weightdf[1,7]
weightlist [subject] <- weight
}
#Creating a dataframe with weights and subj Id
weightdf <- as.data.frame(do.call(rbind,weightlist))
weightdf <- rownames_to_column(weightdf, var="Id")
colnames(weightdf)[1]="Id"
colnames(weightdf)[2]="weight"
subj_quest <- merge(subj_quest,weightdf,by="Id")
#Remove Invalid participans
subj_quest <- subj_quest[subj_quest$isvalid != 0,]Calculating PSDs sequences (spectra)
Since the spectra for each trial and each participant is relatively long to compute, they are saved in a rspec.dat in the data directory. For future runs this file will be used, if present, instead of computing the spectra again.
# Check if the rspec is already available in the directory
#and if not, compute it
if (file.exists(file.path(data_dir, "rspec.dat"))) {
load(file.path(data_dir, "rspec.dat"))
print("rspec.dat already exists, loading it from save...",
quote = F)
} else {
# Looping through the subjects
for (s in seq(1, nrow(subj_quest)))
{ # Loop through the subjects :
# Get the subject id
subject <- subj_quest$Id[s]
# Check the validity of each trial for this subject
valid1 <- subj_quest$valid1[subj_quest$Id == subject] == 1
valid2 <- subj_quest$valid2[subj_quest$Id == subject] == 1
valid3 <- subj_quest$valid3[subj_quest$Id == subject] == 1
# Creating vector mask for validity
vvec <- c(valid1, valid2, valid3)
# Applying the mask to the trials vector
trials <- cond[!!vvec]
for (j in seq(1, length(trials))) { # Loop through the trials :
# Get the conditions
trial <- trials[j]
perct <- round(((s/nrow(subj_quest))*100),2)
filename <- paste0(file_path_sans_ext(
as.character(subject))
, "-", trial,".brute")
cat(paste0(
"\r", " Computing the spectra for the file : ",
filename,
" ( ",
perct,
"% finished )"))
# Get the posturography data
data <- read.table(file.path(data_dir, filename),
sep = "\t",
skip = 1, header = TRUE)
data <- head(data, nb_samp)
y <- data$Barycentre.y..mm.
# Compute the spectra
for (i in seq(1, nb_step/dratio)) {
t <- i * dratio
idx <- seq(t, t + window_size - 1)
p <- psdcore(y [idx])
spec <- 10 * log10(p$spec)
spectra [i, ] <- spec [initial_idx : final_idx]
}
# Extract trial result in a list
rspec [[subject]] [[trial]] <- list(spectra = spectra)
}}
save (rspec, file = file.path (data_dir, "rspec.dat"))
}[1] rspec.dat already exists, loading it from save...cat("\r Done computing spectra ! \n Saved in data/rspec.dat")
Done computing spectra !
Saved in data/rspec.datComputing Energy Ratios
Setting the time and frequency bands of interest for the computing of energy ratios :
# Setting the two frequency bands
bf <- c(1:7) # get the columns corresponding
#to the 0.1 to 0.7 Hz range
hf <- c(7:13) # get the columns corresponding
#to the 0.7 to 1.3 Hz range
# Setting the five time windows
t1 <- c((5000/dratio):(6000/dratio)) # get the rows corresponding
#to the 50s to 60 sec range
t2 <- c((6000/dratio):(7000/dratio)) # get the rows corresponding
#to the 60s to 70 sec range
t3 <- c((11000/dratio):(12000/dratio)) # get the rows corresponding
#to the 110s to 120 sec range
t4 <- c((12000/dratio):(13000/dratio)) # get the rows corresponding
#to the 120s to 130 sec range
t5 <- c((16000/dratio):(17000/dratio)) # get the rows corresponding
#to the 160s to 170 sec rangeComputing Energy ratios for each trials of each participants and extracting them from a list to make a single data frame :
# Looping over each condition for all subjects
for (s in names(rspec)){
for ( c in names (rspec [[s]])) {
# Get low frequency band for each time window
#(Peterka : 0.1 - 0.7 Hz)
l1 <- rspec [[s]][[c]]$spectra[t1, bf]
l2 <- rspec [[s]][[c]]$spectra[t2, bf]
l3 <- rspec [[s]][[c]]$spectra[t3, bf]
l4 <- rspec [[s]][[c]]$spectra[t4, bf]
l5 <- rspec [[s]][[c]]$spectra[t5, bf]
# Get high frequency band for each time window
#(Peterka : 0.7 - 1.3 Hz)
h1 <- rspec [[s]][[c]]$spectra[t1, hf]
h2 <- rspec [[s]][[c]]$spectra[t2, hf]
h3 <- rspec [[s]][[c]]$spectra[t3, hf]
h4 <- rspec [[s]][[c]]$spectra[t4, hf]
h5 <- rspec [[s]][[c]]$spectra[t5, hf]
#Compute Energy Ratio
er1=(sum(h1)/sum(l1))
er2=(sum(h2)/sum(l2))
er3=(sum(h3)/sum(l3))
er4=(sum(h4)/sum(l4))
er5=(sum(h5)/sum(l5))
#Create a list with only these bands and energy ratios
rsptat [[s]] [[c]]<- list(
h1 = h1, # High frequencies over each time windows
h2 = h2,
h3 = h3,
h4 = h4,
h5 = h5,
l1 = l1, # Low frequencies over each time windows
l2 = l2,
l3 = l3,
l4 = l4,
l5 = l5,
er1 = er1, # Energy ratio over each time windows
er2 = er2,
er3 = er3,
er4 = er4,
er5 = er5,
er2_er1 = er2-er1,
er4_er3 = er4-er3,
carsick = subj_quest$CarSick[subj_quest$Id == s],
age = subj_quest$Age[subj_quest$Id == s],
sexe = subj_quest$Sexe[subj_quest$Id == s],
height = subj_quest$Height[subj_quest$Id == s],
weight = subj_quest$weight[subj_quest$Id == s],
validity = subj_quest$isvalid[subj_quest$Id == s])
}}for (s in names (rsptat))
# Loop over vision conditions
for (c in names (rsptat [[s]]))
# Loop over SR conditions
{
datatemp <- data.frame(
subject = s,
condition = c,
er1 = rsptat [[s]] [[c]]$er1,
er2 = rsptat [[s]] [[c]]$er2,
er3 = rsptat [[s]] [[c]]$er3,
er4 = rsptat [[s]] [[c]]$er4,
er5 = rsptat [[s]] [[c]]$er5,
er2_er1 = rsptat [[s]] [[c]]$er2_er1,
er4_er3 = rsptat [[s]] [[c]]$er4_er3,
carsick = rsptat [[s]] [[c]]$carsick,
age = rsptat [[s]] [[c]]$age,
sexe = rsptat [[s]] [[c]]$sexe,
height = rsptat [[s]] [[c]]$height,
weight = rsptat [[s]] [[c]]$weight,
validity = rsptat [[s]] [[c]]$validity)
datastat <- rbind(datastat,datatemp)
}Reformatting the dataframe by averaging ERs for every trials succeeded by each participants :
datamean <- aggregate(list(
er1 = datastat$er1,
er2 = datastat$er2,
er3 = datastat$er3,
er4 = datastat$er4,
er5 = datastat$er5,
er2_er1 = datastat$er2_er1,
er4_er3 = datastat$er4_er3,
carsick = datastat$carsick,
age= datastat$age,
weight = datastat$weight,
validity = datastat$validity
),
by = list(subject = datastat$subject),
FUN = mean)
names(subj_quest)[names(subj_quest) == 'Id'] <- "subject"
sexe <- subj_quest[,c("subject","Sexe")]
datamean <- merge(datamean,sexe,by="subject")Car sickness quartiles computation
Computing quartiles values for the car sickness score and adding quartile categorization in the dataframe :
# Extracting q1 and q3 (quartiles) values for car sickness
q1 <- quantile(datamean$carsick, 0.25)
q3 <- quantile(datamean$carsick, 0.75)datamean <- datamean %>%
mutate(carsick_quart =
ifelse(datamean$carsick < q1,-1,
ifelse(datamean$carsick > q3, 1, 0)),
.after = validity) Saving and displaying the final dataframe
write_xlsx(datamean,"datamean.xlsx")
paged_table(datamean)