In [6]:
# Load some needed packages
library(grid)
library(gridExtra)
library(ggplot2)
library(dplyr)
library(repr)
library(reshape2)
library(tidyr)
library(lme4)
library(utils)
library(multcomp)
library(xtable)
In [ ]:
timestamp_to_numeric <- function(x) {
times=unlist(strsplit(as.character(x), ":"))
out_time = NaN
if ( length(times) > 0 )
{
out_time <- as.numeric(times[1])*60.0
}
if ( length(times) > 1 )
{
out_time <- out_time + as.numeric(times[2])
}
if ( length(times) > 2 )
{
out_time <- out_time + as.numeric(times[3])/100.0
}
return(out_time)
}
In [189]:
abide_smri_scan_params<-read.csv("anat_abide_scan_params_bids.csv")
# reduce to just the Siemens systems
abide_smri_scan_params<-abide_smri_scan_params %>%
filter(Manufacturer == "Siemens") %>%
mutate(Scanner = paste(Manufacturer,ManufacturersModelName))
abide_smri_scan_params$Scanner<-factor(abide_smri_scan_params$Scanner)
abide_smri_scan_params<-abide_smri_scan_params %>%
rowwise() %>%
mutate(VoxelVolume = as.numeric(SliceThickness)*
prod(as.numeric(unlist(strsplit(as.character(PixelSpacing), "x")))))
abide_smri_scan_params<-abide_smri_scan_params %>%
rowwise() %>%
mutate(AcquisitionDuration = timestamp_to_numeric(AcquisitionDuration))
abide_smri_scan_params<-abide_smri_scan_params %>%
rowwise() %>%
mutate(AcquisitionDurationCalc = RepetitionTime * NumberOfSlices)
abide_smri_scan_params<-abide_smri_scan_params %>%
dplyr::select(Site,Scanner,EchoTime,RepetitionTime,InversionTime,FlipAngle,VoxelVolume) %>%
drop_na()
abide_smri_scan_params$Scanner<-droplevels(abide_smri_scan_params$Scanner)
In [190]:
abide_smri_scan_params %>% tbl_df()
In [191]:
abide_anat_spat_df<-read.csv("2016_05_ABIDE_qap_anatomical_spatial.csv")
id.vars=c('Participant','Site','Session','Series')
measure.vars=c('CNR','Cortical.Contrast','EFC','FBER','FWHM','Qi1','SNR')
abide_anat_spat_df<-abide_anat_spat_df[c(id.vars,measure.vars)]
# reduce to just session_1 and anat_1 and remove rows with missing values
abide_anat_spat_df <- abide_anat_spat_df %>%
filter(Session == "session_1" & Series == "anat_1") %>%
drop_na()
# make sure that participant is a factor
abide_anat_spat_df$Participant <- factor(abide_anat_spat_df$Participant)
# summary(abide_anat_spat_df)
# plots
qap_label_strings=c(CNR='CNR',
Cortical.Contrast='Cortical Contrast',
EFC='EFC',
FBER='FBER',
FWHM='Smoothness (FWHM)',
Qi1='Fraction of Artifact Voxels',
SNR='SNR')
# combine with parameters
abide_anat_df <- inner_join(abide_anat_spat_df, abide_smri_scan_params, by="Site")
In [192]:
abide_anat_df %>% tbl_df()
In [193]:
corr_smri_scan_params<-read.csv("anat_corr_scan_params.csv")
corr_smri_scan_params<-corr_smri_scan_params %>%
filter(Manufacturer == "Siemens") %>%
mutate(Scanner = paste(Manufacturer,ManufacturersModelName))
corr_smri_scan_params$Scanner<-factor(corr_smri_scan_params$Scanner)
corr_smri_scan_params<-corr_smri_scan_params %>%
rowwise() %>%
mutate(VoxelVolume = as.numeric(SliceThickness)*
prod(as.numeric(unlist(strsplit(as.character(PixelSpacing), "x")))))
corr_smri_scan_params<-corr_smri_scan_params %>%
filter(PulseSequenceType != "3D MP2RAGE") %>%
filter(PulseSequenceType != "3D MEMPR")
corr_smri_scan_params<-corr_smri_scan_params %>%
dplyr::select(Site,Scanner,EchoTime,RepetitionTime,InversionTime,FlipAngle,VoxelVolume) %>%
drop_na()
# tidy up columns, dropping unused levels and converting cols to numeric
for( i in colnames(corr_smri_scan_params)){
if (( i != "Site" ) & ( i != "Scanner" )){
if (! is.numeric(corr_smri_scan_params[[i]])){
corr_smri_scan_params[[i]]<-as.numeric(levels(corr_smri_scan_params[[i]]))[corr_smri_scan_params[[i]]]
}
}
}
# corr_smri_scan_params$Site<-droplevels(corr_smri_scan_params$Site)
# corr_smri_scan_params$Scanner<-droplevels(corr_smri_scan_params$Scanner)
corr_smri_scan_params %>% tbl_df()
In [15]:
corr_anat_spat_df<-read.csv("2016_05_CORR_qap_anatomical_spatial.csv")
id.vars=c('Participant','Site','Session','Series')
measure.vars=c('CNR','Cortical.Contrast','EFC','FBER','FWHM','Qi1','SNR')
corr_anat_spat_df<-corr_anat_spat_df[c(id.vars,measure.vars)]
# reduce to just session_1 and anat_1 and remove rows with missing values
corr_anat_spat_df <- corr_anat_spat_df %>%
filter(Session == "session_1" & Series == "anat_1") %>%
drop_na()
# make sure that participant is a factor
corr_anat_spat_df$Participant <- factor(corr_anat_spat_df$Participant)
# summary(abide_anat_spat_df)
# plots
qap_label_strings=c(CNR='CNR',
Cortical.Contrast='Cortical Contrast',
EFC='EFC',
FBER='FBER',
FWHM='Smoothness (FWHM)',
Qi1='Fraction of Artifact Voxels',
SNR='SNR')
corr_anat_spat_df <- inner_join(corr_anat_spat_df, corr_smri_scan_params, by="Site")
In [16]:
corr_anat_spat_df %>% tbl_df()
In [17]:
# Combine the twodatasets into one large dataframe
qap_param_df<-bind_rows(abide_anat_df,corr_anat_spat_df)
# mean center all of the covariates
qap_param_df<-qap_param_df %>% mutate_if(is.numeric, funs(scale(., scale=TRUE)))
In [194]:
meas_stats=data.frame()
for ( meas in measure.vars ){
print(meas)
full_formula=as.formula(
sprintf("%s ~ 0 + Scanner + EchoTime + InversionTime + RepetitionTime + FlipAngle + VoxelVolume",
meas))
reduced_formula=as.formula(
sprintf("%s ~ 0 + EchoTime + InversionTime + RepetitionTime + FlipAngle + VoxelVolume",
meas))
full_model<-lm(full_formula, qap_param_df)
reduced_model<-lm(reduced_formula, qap_param_df)
scanner_f=anova(full_model,reduced_model)
all_t=summary(full_model)
meas_stats <- rbind(meas_stats,
data.frame(Measure=meas,
covariate="Scanner",
stat=scanner_f$F[2],
val_p=scanner_f[2,"Pr(>F)"]))
for( covariate in c('EchoTime', 'InversionTime', 'RepetitionTime', 'FlipAngle', 'VoxelVolume')){
meas_stats <- rbind(meas_stats,
data.frame(Measure=meas,
covariate=covariate,
stat=all_t$coefficients[covariate,'t value'],
val_p=all_t$coefficients[covariate,'Pr(>|t|)']))
}
}
meas_stats$val_p <- p.adjust(meas_stats$val_p, "BH")
# smri_qap_param_model_red<-lm(SNR ~ 0 + EchoTime + InversionTime +
# RepetitionTime + FlipAngle + VoxelVolume,
# qap_param_df)
# # summary(smri_qap_param_model)
# # anova(smri_qap_param_model,smri_qap_param_model_red)
In [140]:
meas_stats = meas_stats %>%
gather(mes, value, -Measure, -covariate) %>%
mutate( cov_stat = paste(covariate, mes, sep="_")) %>%
dplyr::select(Measure, cov_stat, value) %>%
spread(cov_stat, value)
In [141]:
meas_stats %>%tbl_df()
In [142]:
print(xtable(meas_stats, digits=3, display=c('s','s','g','g','g','g','g','g','g','g','g','g','g','g')))
In [209]:
abide_fmri_scan_params<-read.csv("rest_abide_scan_params_bids_slicetimes.csv")
# reduce to just the Siemens systems
abide_fmri_scan_params<-abide_fmri_scan_params %>%
filter(Manufacturer == "Siemens") %>%
mutate(Scanner = paste(Manufacturer,ManufacturersModelName))
abide_fmri_scan_params$Scanner<-factor(abide_fmri_scan_params$Scanner)
abide_fmri_scan_params<-abide_fmri_scan_params %>%
rowwise() %>%
mutate(VoxelVolume = as.numeric(SliceThickness)*
prod(as.numeric(unlist(strsplit(as.character(PixelSpacing), "x")))))
abide_fmri_scan_params<-abide_fmri_scan_params %>%
rowwise() %>%
mutate(AcquisitionDuration = timestamp_to_numeric(AcquisitionDuration))
abide_fmri_scan_params<-abide_fmri_scan_params %>%
dplyr::select(Site, Scanner, AcquisitionDuration, EchoTime ,RepetitionTime, FlipAngle, VoxelVolume) %>%
drop_na()
abide_fmri_scan_params$Scanner<-droplevels(abide_fmri_scan_params$Scanner)
In [210]:
abide_fmri_scan_params %>% tbl_df()
In [211]:
abide_func_spat_df<-read.csv("2016_05_ABIDE_qap_functional_spatial.csv")
id.vars=c('Participant','Site','Session','Series')
measure.vars.spat=c('EFC','FBER','FWHM','Ghost_y','SNR')
# reduce to just the columns that we are interested in
abide_func_spat_df<-abide_func_spat_df[c(id.vars,measure.vars.spat)]
abide_func_temp_df<-read.csv("2016_05_ABIDE_qap_functional_temporal.csv")
id.vars=c('Participant','Site','Session','Series')
measure.vars.temp=c('Fraction.of.Outliers..Mean.','GCOR','Quality..Mean.',
'RMSD..Mean.', 'Std..DVARS..Mean.')
# reduce to just the columns that we are interested in
abide_func_temp_df<-abide_func_temp_df[c(id.vars, measure.vars.temp)]
# join the spatial and temporal measures
abide_func_df<-inner_join(abide_func_spat_df, abide_func_temp_df)
# reduce to just session_1 and anat_1 and remove rows with missing values
abide_func_df <- abide_func_df %>%
filter(Session == "session_1" & Series == "rest_1") %>%
droplevels() %>%
drop_na()
# make sure that participant is a factor
abide_func_df$Participant <- factor(abide_func_df$Participant)
# plots
qap_label_strings=c(Ghost_y='GSR',
EFC='EFC',
FBER='FBER',
FWHM='Smoothness (FWHM)',
SNR='SNR',
Fraction.of.Outliers..Mean.='Mean % Outliers',
GCOR='GCOR',
Quality..Mean.='Mean Quality',
RMSD..Mean.='Mean RMSD',
Std..DVARS..Mean.="Mean Std. Dvars.")
# combine with parameters
abide_func_df <- inner_join(abide_func_df, abide_fmri_scan_params, by="Site")
summary(abide_func_df)
In [214]:
corr_fmri_scan_params<-read.csv("func_corr_scan_params.csv")
# reduce to just the Siemens systems
corr_fmri_scan_params<-corr_fmri_scan_params %>%
filter(Manufacturer == "Siemens") %>%
mutate(Scanner = paste(Manufacturer,ManufacturersModelName))
corr_fmri_scan_params$Scanner<-factor(corr_fmri_scan_params$Scanner)
corr_fmri_scan_params<-corr_fmri_scan_params %>%
rowwise() %>%
mutate(VoxelVolume = as.numeric(SliceThickness)*
prod(as.numeric(unlist(strsplit(as.character(PixelSpacing), "x")))))
corr_fmri_scan_params<-corr_fmri_scan_params %>%
rowwise() %>%
mutate(AcquisitionDuration = timestamp_to_numeric(AcquisitionDuration))
corr_fmri_scan_params<-corr_fmri_scan_params %>%
dplyr::select(Site, Scanner, AcquisitionDuration, EchoTime ,RepetitionTime, FlipAngle, VoxelVolume) %>%
drop_na()
corr_fmri_scan_params$Scanner<-droplevels(corr_fmri_scan_params$Scanner)
In [215]:
summary(corr_fmri_scan_params)
In [216]:
corr_func_spat_df<-read.csv("2016_05_CoRR_qap_functional_spatial.csv")
id.vars=c('Participant','Site','Session','Series')
# reduce to just the columns that we are interested in
corr_func_spat_df<-corr_func_spat_df[c(id.vars,measure.vars.spat)]
corr_func_temp_df<-read.csv("2016_05_CORR_qap_functional_temporal.csv")
# reduce to just the columns that we are interested in
corr_func_temp_df<-corr_func_temp_df[c(id.vars, measure.vars.temp)]
# join the spatial and temporal measures
corr_func_df<-inner_join(corr_func_spat_df, corr_func_temp_df)
# reduce to just session_1 and anat_1 and remove rows with missing values
corr_func_df <- corr_func_df %>%
filter(Session == "session_1" & Series == "rest_1") %>%
droplevels() %>%
drop_na()
# make sure that participant is a factor
corr_func_df$Participant <- factor(corr_func_df$Participant)
# plots
qap_label_strings=c(Ghost_y='GSR',
EFC='EFC',
FBER='FBER',
FWHM='Smoothness (FWHM)',
SNR='SNR',
Fraction.of.Outliers..Mean.='Mean % Outliers',
GCOR='GCOR',
Quality..Mean.='Mean Quality',
RMSD..Mean.='Mean RMSD',
Std..DVARS..Mean.="Mean Std. Dvars.")
# combine with parameters
corr_func_df <- inner_join(corr_func_df, corr_fmri_scan_params, by="Site")
summary(corr_func_df)
In [217]:
# Combine the twodatasets into one large dataframe
qap_param_df<-bind_rows(abide_func_df,corr_func_spat_df)
# mean center all of the covariates
qap_param_df<-qap_param_df %>% mutate_if(is.numeric, funs(scale(., scale=TRUE)))
In [203]:
qap_param_df
In [218]:
meas_stats=data.frame()
measure.vars=c(measure.vars.spat, measure.vars.temp)
for ( meas in measure.vars ){
print(meas)
full_formula=as.formula(
sprintf("%s ~ 0 + Scanner + EchoTime + RepetitionTime + FlipAngle + VoxelVolume",
meas))
reduced_formula=as.formula(
sprintf("%s ~ 0 + EchoTime + RepetitionTime + FlipAngle + VoxelVolume",
meas))
full_model<-lm(full_formula, qap_param_df)
reduced_model<-lm(reduced_formula, qap_param_df)
scanner_f=anova(full_model,reduced_model)
all_t=summary(full_model)
meas_stats <- rbind(meas_stats,
data.frame(Measure=meas,
covariate="Scanner",
stat=scanner_f$F[2],
val_p=scanner_f[2,"Pr(>F)"]))
for( covariate in c('EchoTime', 'RepetitionTime', 'FlipAngle', 'VoxelVolume')){
meas_stats <- rbind(meas_stats,
data.frame(Measure=meas,
covariate=covariate,
stat=all_t$coefficients[covariate,'t value'],
val_p=all_t$coefficients[covariate,'Pr(>|t|)']))
}
}
meas_stats$val_p <- p.adjust(meas_stats$val_p, "BH")
# smri_qap_param_model_red<-lm(SNR ~ 0 + EchoTime + InversionTime +
# RepetitionTime + FlipAngle + VoxelVolume,
# qap_param_df)
# # summary(smri_qap_param_model)
# # anova(smri_qap_param_model,smri_qap_param_model_red)
In [219]:
meas_stats = meas_stats %>%
gather(mes, value, -Measure, -covariate) %>%
mutate( cov_stat = paste(covariate, mes, sep="_")) %>%
dplyr::select(Measure, cov_stat, value) %>%
spread(cov_stat, value)
In [220]:
meas_stats %>% tbl_df()
In [222]:
print(xtable(meas_stats, digits=3, display=c('s','s','g','g','g','g','g','g','g','g','g','g')))