Permalink
Switch branches/tags
Nothing to show
Find file Copy path
Fetching contributors…
Cannot retrieve contributors at this time
873 lines (593 sloc) 28 KB
title author date output
Timeseries analyses for MyConnectome Project
Russ Poldrack
November 1, 2014
html_document

Timeseries analyses for MyConnectome dataset

Code available at: https://github.com/poldrack/myconnectome/blob/master/myconnectome/timeseries/timeseries_analyses.Rmd

This notebook runs all of the timeseries analyses presented in the MyConnectome paper. The tables list all tests that survived at FDR-corrected p < 0.1, correcting within each set of tests.

The variables listed in each table are:

  • xvar:name of X variable
  • yvar: name of Y variable
  • cor.val: Pearson correlation between X and Y variables
  • t.arima: t value from ARIMA model for X regressor
  • arima.p: uncorrected p value from ARIMA model
  • t.drift: t value from ARIMA model for drift parameter
  • drift.p: p value from ARIMA model for drift parameter
  • AR: autoregressive model order selected by automatic ARIMA model selection
  • MA: moving average model order selected by automatic ARIMA model selection
  • nobs: number of observations contributing to test (differs across tests due to missing data)
  • pval_bh: p-value after Benajmini-Hochberg false discovery rate correction
library(forecast)
library(knitr)
library(zoo)


basedir=Sys.getenv('MYCONNECTOME_DIR')

thresh=0.1  # FDR threshold
save_latex=TRUE
OUTPUT_DIR=sprintf('%s/timeseries',basedir)

behav=load_behav_data()
wincorr=load_fmri_data('wincorr')
bwcorr=load_fmri_data('bwcorr')
netdat=load_network_data()
fd=load_fd_data()
food=load_food_data()
rnaseq_wgcna=load_rnaseq_data(limit_ME_to_enriched=FALSE)
metab=load_metab_data()
pindex=load_participation_index()

data_names=c()
for (i in 1:(dim(metab)[2]-1)) {
    data_names=rbind(data_names,sprintf('C%d:%s',i,names(metab)[i]))
  	}
data_names=rbind(data_names,'date')
names(metab)=data_names


fullmetab=load_metab_data(use_clustered_data=FALSE)
data_names=c()
for (i in 1:(dim(fullmetab)[2]-1)) {
    n=gsub('_NIST','',names(fullmetab)[i])
    n=gsub('_',' ',n)

    data_names=rbind(data_names,n)
    }
data_names=rbind(data_names,'date')
names(fullmetab)=data_names


immport=load_ImmPort_data()

# Create the subset of data that we are interested in for the behavioral analysis - generate a version without the tuesday/thursday variable for use with the blood data (which are only on tuesdays)

xvars=c('panas.positive','panas.negative','panas.fatigue','afterscan.Anxietyduringscan','afterscan.diastolic','afterscan.pulse','afterscan.systolic','morning.Sleepquality','morning.Soreness','prevevening.Alcohol','prevevening.Guthealth','prevevening.Psoriasisseverity','prevevening.Stress', 'prevevening.Timespentoutdoors','TuesThurs', 'temp.mean',"email.LIWCcdi","email.LIWCnegemo","email.LIWCposemo",'zeo.zq')

behav_keep=subset(behav,select=c(xvars,'date'))
behav_keep_no_tth=subset(behav_keep,select=-TuesThurs)

Behavioral variables vs. each other:

vnames=c('behav','behav')
if (file.exists(paste(OUTPUT_DIR,sprintf('out.dat.%s_%s.txt',vnames[1],vnames[2]),sep='/'))) {
  #print("loading existing file...")
  d=read.table(paste(OUTPUT_DIR,sprintf('out.dat.%s_%s.txt',vnames[1],vnames[2]),sep='/'))
  } else {
d=est_bivariate_arima_model(behav_keep, behav_keep,spacing='1 day',skip_ident=TRUE)
}
kable_wrap(d)
output_results(d,vnames,OUTPUT_DIR=OUTPUT_DIR)

Within-module connectivity vs. behavior

vnames=c('wincorr','behav')

if (file.exists(paste(OUTPUT_DIR,sprintf('out.dat.%s_%s.txt',vnames[1],vnames[2]),sep='/'))) {
  #print("loading existing file...")
  d=read.table(paste(OUTPUT_DIR,sprintf('out.dat.%s_%s.txt',vnames[1],vnames[2]),sep='/'))
  } else {
    d=est_bivariate_arima_model(wincorr,behav_keep)
    }

kable_wrap(d)
output_results(d,vnames,OUTPUT_DIR=OUTPUT_DIR)

Between-module connectivity vs. behavior

vnames=c('bwcorr','behav')

if (file.exists(paste(OUTPUT_DIR,sprintf('out.dat.%s_%s.txt',vnames[1],vnames[2]),sep='/'))) {
  #print("loading existing file...")
  d=read.table(paste(OUTPUT_DIR,sprintf('out.dat.%s_%s.txt',vnames[1],vnames[2]),sep='/'))
  } else {d=est_bivariate_arima_model(bwcorr,behav_keep)}

kable_wrap(d)
output_results(d,vnames,OUTPUT_DIR=OUTPUT_DIR)

Network measures vs. behavior

vnames=c('netdat','behav')

if (file.exists(paste(OUTPUT_DIR,sprintf('out.dat.%s_%s.txt',vnames[1],vnames[2]),sep='/'))) {
  #print("loading existing file...")
  d=read.table(paste(OUTPUT_DIR,sprintf('out.dat.%s_%s.txt',vnames[1],vnames[2]),sep='/'))
  } else {d=est_bivariate_arima_model(netdat,behav_keep)}

kable_wrap(d)
output_results(d,vnames,OUTPUT_DIR=OUTPUT_DIR)

Participation index vs. behavior

vnames=c('pindex','behav')

if (file.exists(paste(OUTPUT_DIR,sprintf('out.dat.%s_%s.txt',vnames[1],vnames[2]),sep='/'))) {
  #print("loading existing file...")
  d=read.table(paste(OUTPUT_DIR,sprintf('out.dat.%s_%s.txt',vnames[1],vnames[2]),sep='/'))
  } else {d=est_bivariate_arima_model(pindex,behav_keep)}
print('skipping printout')

output_results(d,vnames,OUTPUT_DIR=OUTPUT_DIR)

Behavior vs. RNA-seq:

vnames=c('wgcna','behav')

if (file.exists(paste(OUTPUT_DIR,sprintf('out.dat.%s_%s.txt',vnames[1],vnames[2]),sep='/'))) {
  #print("loading existing file...")
  d=read.table(paste(OUTPUT_DIR,sprintf('out.dat.%s_%s.txt',vnames[1],vnames[2]),sep='/'))
  } else {d=est_bivariate_arima_model(rnaseq_wgcna,behav_keep_no_tth,spacing='1 week')}

kable_wrap(d)
output_results(d,vnames,OUTPUT_DIR=OUTPUT_DIR)


Within-module connectivity vs. RNA-seq

vnames=c('wgcna','wincorr')

if (file.exists(paste(OUTPUT_DIR,sprintf('out.dat.%s_%s.txt',vnames[1],vnames[2]),sep='/'))) {
  #print("loading existing file...")
  d=read.table(paste(OUTPUT_DIR,sprintf('out.dat.%s_%s.txt',vnames[1],vnames[2]),sep='/'))
  } else {d=est_bivariate_arima_model(rnaseq_wgcna,wincorr,spacing='1 week')}

kable_wrap(d)
output_results(d,vnames,OUTPUT_DIR=OUTPUT_DIR)

Between-module connectivity vs. RNA-seq

vnames=c('wgcna','bwcorr')

if (file.exists(paste(OUTPUT_DIR,sprintf('out.dat.%s_%s.txt',vnames[1],vnames[2]),sep='/'))) {
  #print("loading existing file...")
  d=read.table(paste(OUTPUT_DIR,sprintf('out.dat.%s_%s.txt',vnames[1],vnames[2]),sep='/'))
  } else {d=est_bivariate_arima_model(rnaseq_wgcna,bwcorr,spacing='1 week')}

kable_wrap(d)
output_results(d,vnames,OUTPUT_DIR=OUTPUT_DIR)

Network measures vs. RNA-seq

vnames=c('wgcna','netdat')

if (file.exists(paste(OUTPUT_DIR,sprintf('out.dat.%s_%s.txt',vnames[1],vnames[2]),sep='/'))) {
  #print("loading existing file...")
  d=read.table(paste(OUTPUT_DIR,sprintf('out.dat.%s_%s.txt',vnames[1],vnames[2]),sep='/'))
  } else {d=est_bivariate_arima_model(rnaseq_wgcna,netdat,spacing='1 week')}

kable_wrap(d)
output_results(d,vnames,OUTPUT_DIR=OUTPUT_DIR)

ImmPort immune genes vs behavior

vnames=c('immport','behav')

if (file.exists(paste(OUTPUT_DIR,sprintf('out.dat.%s_%s.txt',vnames[1],vnames[2]),sep='/'))) {
  #print("loading existing file...")
  d=read.table(paste(OUTPUT_DIR,sprintf('out.dat.%s_%s.txt',vnames[1],vnames[2]),sep='/'))
  } else {d=est_bivariate_arima_model(immport,behav_keep_no_tth,spacing='1 week')}

kable_wrap(d)
output_results(d,vnames,OUTPUT_DIR=OUTPUT_DIR)

Test specific hypothesis regarding psoriasis

vnames=c('psoriasis','immport')
psoriasis=subset(behav,select=c(date,prevevening.Psoriasisseverity))
if (file.exists(paste(OUTPUT_DIR,sprintf('out.dat.%s_%s.txt',vnames[1],vnames[2]),sep='/'))) {
  #print("loading existing file...")
  d=read.table(paste(OUTPUT_DIR,sprintf('out.dat.%s_%s.txt',vnames[1],vnames[2]),sep='/'))
  } else {d=est_bivariate_arima_model(immport,psoriasis,spacing='1 week')}

kable_wrap(d)
output_results(d,vnames,OUTPUT_DIR=OUTPUT_DIR)

Test specific hypothesis regarding mood

vnames=c('mood','immport')
mood=subset(behav,select=c(date,panas.positive,panas.negative))
if (file.exists(paste(OUTPUT_DIR,sprintf('out.dat.%s_%s.txt',vnames[1],vnames[2]),sep='/'))) {
  #print("loading existing file...")
  d=read.table(paste(OUTPUT_DIR,sprintf('out.dat.%s_%s.txt',vnames[1],vnames[2]),sep='/'))
  } else {d=est_bivariate_arima_model(immport,mood,spacing='1 week')}

kable_wrap(d)
output_results(d,vnames,OUTPUT_DIR=OUTPUT_DIR)

ImmPort immune genes vs within-module connectivity

vnames=c('immport','wincorr')

if (file.exists(paste(OUTPUT_DIR,sprintf('out.dat.%s_%s.txt',vnames[1],vnames[2]),sep='/'))) {
  #print("loading existing file...")
  d=read.table(paste(OUTPUT_DIR,sprintf('out.dat.%s_%s.txt',vnames[1],vnames[2]),sep='/'))
  } else {d=est_bivariate_arima_model(immport,wincorr,spacing='1 week')}

kable_wrap(d)
output_results(d,vnames,OUTPUT_DIR=OUTPUT_DIR)

ImmPort immune genes vs between-module connectivity

vnames=c('immport','bwcorr')

if (file.exists(paste(OUTPUT_DIR,sprintf('out.dat.%s_%s.txt',vnames[1],vnames[2]),sep='/'))) {
  #print("loading existing file...")
  d=read.table(paste(OUTPUT_DIR,sprintf('out.dat.%s_%s.txt',vnames[1],vnames[2]),sep='/'))
  } else {d=est_bivariate_arima_model(immport,bwcorr,spacing='1 week')}

kable_wrap(d)
output_results(d,vnames,OUTPUT_DIR=OUTPUT_DIR)

ImmPort immune genes vs network measures

vnames=c('immport','netdat')

if (file.exists(paste(OUTPUT_DIR,sprintf('out.dat.%s_%s.txt',vnames[1],vnames[2]),sep='/'))) {
  #print("loading existing file...")
  d=read.table(paste(OUTPUT_DIR,sprintf('out.dat.%s_%s.txt',vnames[1],vnames[2]),sep='/'))
  } else {d=est_bivariate_arima_model(immport,netdat,spacing='1 week')}

kable_wrap(d)
output_results(d,vnames,OUTPUT_DIR=OUTPUT_DIR)

ImmPort immune genes vs food

vnames=c('food','immport')

if (file.exists(paste(OUTPUT_DIR,sprintf('out.dat.%s_%s.txt',vnames[1],vnames[2]),sep='/'))) {
  #print("loading existing file...")
  d=read.table(paste(OUTPUT_DIR,sprintf('out.dat.%s_%s.txt',vnames[1],vnames[2]),sep='/'))
  } else {d=est_bivariate_arima_model(food,immport,spacing='1 week')}

kable_wrap(d)
output_results(d,vnames,OUTPUT_DIR=OUTPUT_DIR)

ImmPort immune genes vs metabolism

vnames=c('immport','metab')

if (file.exists(paste(OUTPUT_DIR,sprintf('out.dat.%s_%s.txt',vnames[1],vnames[2]),sep='/'))) {
  #print("loading existing file...")
  d=read.table(paste(OUTPUT_DIR,sprintf('out.dat.%s_%s.txt',vnames[1],vnames[2]),sep='/'))
  } else {d=est_bivariate_arima_model(immport,metab,spacing='1 week')}

kable_wrap(d)
output_results(d,vnames,OUTPUT_DIR=OUTPUT_DIR)

ImmPort immune genes vs metabolism (full set)

vnames=c('immport','fullmetab')

if (file.exists(paste(OUTPUT_DIR,sprintf('out.dat.%s_%s.txt',vnames[1],vnames[2]),sep='/'))) {
  #print("loading existing file...")
  d=read.table(paste(OUTPUT_DIR,sprintf('out.dat.%s_%s.txt',vnames[1],vnames[2]),sep='/'))
  } else {d=est_bivariate_arima_model(immport,fullmetab,spacing='1 week')}

kable_wrap(d)
output_results(d,vnames,OUTPUT_DIR=OUTPUT_DIR)

food vs. metabolism

note: food needs to be X variable because auto.arima doesn't handle binary Y variable

vnames=c('food','metab')

if (file.exists(paste(OUTPUT_DIR,sprintf('out.dat.%s_%s.txt',vnames[1],vnames[2]),sep='/'))) {
  #print("loading existing file...")
  d=read.table(paste(OUTPUT_DIR,sprintf('out.dat.%s_%s.txt',vnames[1],vnames[2]),sep='/'))
  } else {d=est_bivariate_arima_model(food,metab,spacing='1 week')}

kable_wrap(d)
output_results(d,vnames,OUTPUT_DIR=OUTPUT_DIR)

food vs. metabolism (individual metabolites)

note: food needs to be X variable because auto.arima doesn't handle binary Y variable


vnames=c('food','fullmetab')

if (file.exists(paste(OUTPUT_DIR,sprintf('out.dat.%s_%s.txt',vnames[1],vnames[2]),sep='/'))) {
  #print("loading existing file...")
  d=read.table(paste(OUTPUT_DIR,sprintf('out.dat.%s_%s.txt',vnames[1],vnames[2]),sep='/'))
  } else {d=est_bivariate_arima_model(food,fullmetab,spacing='1 week')}

kable_wrap(d)
output_results(d,vnames,OUTPUT_DIR=OUTPUT_DIR)

behav vs. metabolsm


vnames=c('behav','metab')

if (file.exists(paste(OUTPUT_DIR,sprintf('out.dat.%s_%s.txt',vnames[1],vnames[2]),sep='/'))) {
  #print("loading existing file...")
  d=read.table(paste(OUTPUT_DIR,sprintf('out.dat.%s_%s.txt',vnames[1],vnames[2]),sep='/'))
  } else {d=est_bivariate_arima_model(behav_keep_no_tth,metab,spacing='1 week')}

kable_wrap(d)
output_results(d,vnames,OUTPUT_DIR=OUTPUT_DIR)

behav vs. metabolsm (full set)

vnames=c('behav','fullmetab')

if (file.exists(paste(OUTPUT_DIR,sprintf('out.dat.%s_%s.txt',vnames[1],vnames[2]),sep='/'))) {
  #print("loading existing file...")
  d=read.table(paste(OUTPUT_DIR,sprintf('out.dat.%s_%s.txt',vnames[1],vnames[2]),sep='/'))
  } else {d=est_bivariate_arima_model(behav_keep_no_tth,fullmetab,spacing='1 week')}

kable_wrap(d)
output_results(d,vnames,OUTPUT_DIR=OUTPUT_DIR)

metabolism vs. within-network connectivity


vnames=c('metab','wincorr')

if (file.exists(paste(OUTPUT_DIR,sprintf('out.dat.%s_%s.txt',vnames[1],vnames[2]),sep='/'))) {
  #print("loading existing file...")
  d=read.table(paste(OUTPUT_DIR,sprintf('out.dat.%s_%s.txt',vnames[1],vnames[2]),sep='/'))
  } else {d=est_bivariate_arima_model(metab,wincorr,spacing='1 week')}

kable_wrap(d)
output_results(d,vnames,OUTPUT_DIR=OUTPUT_DIR)

metabolism vs. within-network connectivity (full set of metabolites)


vnames=c('fullmetab','wincorr')

if (file.exists(paste(OUTPUT_DIR,sprintf('out.dat.%s_%s.txt',vnames[1],vnames[2]),sep='/'))) {
  #print("loading existing file...")
  d=read.table(paste(OUTPUT_DIR,sprintf('out.dat.%s_%s.txt',vnames[1],vnames[2]),sep='/'))
  } else {d=est_bivariate_arima_model(fullmetab,wincorr,spacing='1 week')}

kable_wrap(d)
output_results(d,vnames,OUTPUT_DIR=OUTPUT_DIR)

metabolism vs. between-network connectivity


vnames=c('metab','bwcorr')

if (file.exists(paste(OUTPUT_DIR,sprintf('out.dat.%s_%s.txt',vnames[1],vnames[2]),sep='/'))) {
  #print("loading existing file...")
  d=read.table(paste(OUTPUT_DIR,sprintf('out.dat.%s_%s.txt',vnames[1],vnames[2]),sep='/'))
  } else {d=est_bivariate_arima_model(metab,bwcorr,spacing='1 week')}

kable_wrap(d)
output_results(d,vnames,OUTPUT_DIR=OUTPUT_DIR)

metabolism vs. within-network connectivity (full set of metabolites)


vnames=c('fullmetab','bwcorr')

if (file.exists(paste(OUTPUT_DIR,sprintf('out.dat.%s_%s.txt',vnames[1],vnames[2]),sep='/'))) {
  #print("loading existing file...")
  d=read.table(paste(OUTPUT_DIR,sprintf('out.dat.%s_%s.txt',vnames[1],vnames[2]),sep='/'))
  } else {d=est_bivariate_arima_model(fullmetab,bwcorr,spacing='1 week')}

kable_wrap(d)
output_results(d,vnames,OUTPUT_DIR=OUTPUT_DIR)

metabolism vs. metabolism


vnames=c('metab','metab')

if (file.exists(paste(OUTPUT_DIR,sprintf('out.dat.%s_%s.txt',vnames[1],vnames[2]),sep='/'))) {
  #print("loading existing file...")
  d=read.table(paste(OUTPUT_DIR,sprintf('out.dat.%s_%s.txt',vnames[1],vnames[2]),sep='/'))
  } else {d=est_bivariate_arima_model(metab,metab,skip_ident=TRUE,spacing='1 week')}

kable_wrap(d)
output_results(d,vnames,OUTPUT_DIR=OUTPUT_DIR)

within-network connectivity

vnames=c('wincorr','wincorr')

if (file.exists(paste(OUTPUT_DIR,sprintf('out.dat.%s_%s.txt',vnames[1],vnames[2]),sep='/'))) {
  #print("loading existing file...")
  d=read.table(paste(OUTPUT_DIR,sprintf('out.dat.%s_%s.txt',vnames[1],vnames[2]),sep='/'))
  } else {d=est_bivariate_arima_model(wincorr,wincorr,skip_ident=TRUE)}

kable_wrap(d)
output_results(d,vnames,OUTPUT_DIR=OUTPUT_DIR)

between-network connectivity

vnames=c('bwcorr','bwcorr')

if (file.exists(paste(OUTPUT_DIR,sprintf('out.dat.%s_%s.txt',vnames[1],vnames[2]),sep='/'))) {
  #print("loading existing file...")
  d=read.table(paste(OUTPUT_DIR,sprintf('out.dat.%s_%s.txt',vnames[1],vnames[2]),sep='/'))
  } else {d=est_bivariate_arima_model(bwcorr,bwcorr,skip_ident=TRUE)}
print('skipping printout')
output_results(d,vnames,OUTPUT_DIR=OUTPUT_DIR)

within-network connectivity vs. network measures

vnames=c('wincorr','netdat')

if (file.exists(paste(OUTPUT_DIR,sprintf('out.dat.%s_%s.txt',vnames[1],vnames[2]),sep='/'))) {
  #print("loading existing file...")
  d=read.table(paste(OUTPUT_DIR,sprintf('out.dat.%s_%s.txt',vnames[1],vnames[2]),sep='/'))
  } else {d=est_bivariate_arima_model(wincorr,netdat)}

kable_wrap(d)
output_results(d,vnames,OUTPUT_DIR=OUTPUT_DIR)

between-network connectivity vs. network measures

vnames=c('bwcorr','netdat')

if (file.exists(paste(OUTPUT_DIR,sprintf('out.dat.%s_%s.txt',vnames[1],vnames[2]),sep='/'))) {
  #print("loading existing file...")
  d=read.table(paste(OUTPUT_DIR,sprintf('out.dat.%s_%s.txt',vnames[1],vnames[2]),sep='/'))
  } else {d=est_bivariate_arima_model(bwcorr,netdat)}

kable_wrap(d)
output_results(d,vnames,OUTPUT_DIR=OUTPUT_DIR)

network measures

vnames=c('netdat','netdat')

if (file.exists(paste(OUTPUT_DIR,sprintf('out.dat.%s_%s.txt',vnames[1],vnames[2]),sep='/'))) {
  #print("loading existing file...")
  d=read.table(paste(OUTPUT_DIR,sprintf('out.dat.%s_%s.txt',vnames[1],vnames[2]),sep='/'))
  } else {d=est_bivariate_arima_model(netdat,netdat,skip_ident=TRUE)}

kable_wrap(d)
output_results(d,vnames,OUTPUT_DIR=OUTPUT_DIR)

food vs. within-network connectivity

vnames=c('food','wincorr')

if (file.exists(paste(OUTPUT_DIR,sprintf('out.dat.%s_%s.txt',vnames[1],vnames[2]),sep='/'))) {
  #print("loading existing file...")
  d=read.table(paste(OUTPUT_DIR,sprintf('out.dat.%s_%s.txt',vnames[1],vnames[2]),sep='/'))
  } else {d=est_bivariate_arima_model(food,wincorr)}

kable_wrap(d)
output_results(d,vnames,OUTPUT_DIR=OUTPUT_DIR)

food vs. behavior

Remove food variables that have too many zeros or ones, otherwise will fail

food$apples=NULL  
food$blueberries=NULL
behav_keep_no_tth$zeo.zq=NULL

vnames=c('food','behav')

if (file.exists(paste(OUTPUT_DIR,sprintf('out.dat.%s_%s.txt',vnames[1],vnames[2]),sep='/'))) {
  #print("loading existing file...")
  d=read.table(paste(OUTPUT_DIR,sprintf('out.dat.%s_%s.txt',vnames[1],vnames[2]),sep='/'))
  } else {d=est_bivariate_arima_model(food,behav_keep_no_tth)}

kable_wrap(d)
output_results(d,vnames,OUTPUT_DIR=OUTPUT_DIR)

compare head motion to other measures

FD vs. behavior


vnames=c('fd','behav')

if (file.exists(paste(OUTPUT_DIR,sprintf('out.dat.%s_%s.txt',vnames[1],vnames[2]),sep='/'))) {
  #print("loading existing file...")
  d=read.table(paste(OUTPUT_DIR,sprintf('out.dat.%s_%s.txt',vnames[1],vnames[2]),sep='/'))
  } else {d=est_bivariate_arima_model(fd,behav_keep)}

kable_wrap(d)
output_results(d,vnames,OUTPUT_DIR=OUTPUT_DIR)

FD vs. RNA-seq

vnames=c('fd','wgcna')

if (file.exists(paste(OUTPUT_DIR,sprintf('out.dat.%s_%s.txt',vnames[1],vnames[2]),sep='/'))) {
  #print("loading existing file...")
  d=read.table(paste(OUTPUT_DIR,sprintf('out.dat.%s_%s.txt',vnames[1],vnames[2]),sep='/'))
  } else {d=est_bivariate_arima_model(fd,rnaseq_wgcna,spacing='1 week')}

kable_wrap(d)
output_results(d,vnames,OUTPUT_DIR=OUTPUT_DIR)

FD vs. metabolites


vnames=c('fd','metab')

if (file.exists(paste(OUTPUT_DIR,sprintf('out.dat.%s_%s.txt',vnames[1],vnames[2]),sep='/'))) {
  #print("loading existing file...")
  d=read.table(paste(OUTPUT_DIR,sprintf('out.dat.%s_%s.txt',vnames[1],vnames[2]),sep='/'))
  } else {d=est_bivariate_arima_model(fd,metab,spacing='1 week')}

kable_wrap(d)
output_results(d,vnames,OUTPUT_DIR=OUTPUT_DIR)

FD vs. wincorr


vnames=c('fd','wincorr')

if (file.exists(paste(OUTPUT_DIR,sprintf('out.dat.%s_%s.txt',vnames[1],vnames[2]),sep='/'))) {
  #print("loading existing file...")
  d=read.table(paste(OUTPUT_DIR,sprintf('out.dat.%s_%s.txt',vnames[1],vnames[2]),sep='/'))
  } else {d=est_bivariate_arima_model(fd,wincorr)}

kable_wrap(d)
output_results(d,vnames,OUTPUT_DIR=OUTPUT_DIR)

FD vs. netdat


vnames=c('fd','netdat')

if (file.exists(paste(OUTPUT_DIR,sprintf('out.dat.%s_%s.txt',vnames[1],vnames[2]),sep='/'))) {
  #print("loading existing file...")
  d=read.table(paste(OUTPUT_DIR,sprintf('out.dat.%s_%s.txt',vnames[1],vnames[2]),sep='/'))
  } else {d=est_bivariate_arima_model(fd,netdat)}

kable_wrap(d)
output_results(d,vnames,OUTPUT_DIR=OUTPUT_DIR)

RNA-seq


vnames=c('wgcna','wgcna')

if (file.exists(paste(OUTPUT_DIR,sprintf('out.dat.%s_%s.txt',vnames[1],vnames[2]),sep='/'))) {
  #print("loading existing file...")
  d=read.table(paste(OUTPUT_DIR,sprintf('out.dat.%s_%s.txt',vnames[1],vnames[2]),sep='/'))
  } else {d=est_bivariate_arima_model(rnaseq_wgcna, rnaseq_wgcna,spacing='1 week',skip_ident=TRUE)}
# turning off output because there are too many significant results here - uncomment to turn it back on
#kable_wrap(d)
output_results(d,vnames,OUTPUT_DIR=OUTPUT_DIR)

RNA-seq vs. metabolism

vnames=c('wgcna','metab')

if (file.exists(paste(OUTPUT_DIR,sprintf('out.dat.%s_%s.txt',vnames[1],vnames[2]),sep='/'))) {
  #print("loading existing file...")
  d=read.table(paste(OUTPUT_DIR,sprintf('out.dat.%s_%s.txt',vnames[1],vnames[2]),sep='/'))
  } else {d=est_bivariate_arima_model(rnaseq_wgcna,metab,spacing='1 week')}

kable_wrap(d)
output_results(d,vnames,OUTPUT_DIR=OUTPUT_DIR)

RNA-seq vs. metabolism (full set)


vnames=c('wgcna','fullmetab')

if (file.exists(paste(OUTPUT_DIR,sprintf('out.dat.%s_%s.txt',vnames[1],vnames[2]),sep='/'))) {
  #print("loading existing file...")
  d=read.table(paste(OUTPUT_DIR,sprintf('out.dat.%s_%s.txt',vnames[1],vnames[2]),sep='/'))
  } else {d=est_bivariate_arima_model(rnaseq_wgcna,fullmetab,spacing='1 week')}

kable_wrap(d)
output_results(d,vnames,OUTPUT_DIR=OUTPUT_DIR)

RNA-seq vs. food


vnames=c('food','wgcna')

if (file.exists(paste(OUTPUT_DIR,sprintf('out.dat.%s_%s.txt',vnames[1],vnames[2]),sep='/'))) {
  #print("loading existing file...")
  d=read.table(paste(OUTPUT_DIR,sprintf('out.dat.%s_%s.txt',vnames[1],vnames[2]),sep='/'))
  } else {d=est_bivariate_arima_model(food,rnaseq_wgcna,spacing='1 week')}

kable_wrap(d)
output_results(d,vnames,OUTPUT_DIR=OUTPUT_DIR)