Full Covariance Regression Applied to HCP Data

Goal

Apply full covariance regression model to preprocessed timeseries data from the Human Connectome Project.

library(CovRegFC)

## Loading required package: Rcpp

## Warning: replacing previous import 'cowplot::ggsave' by 'ggplot2::ggsave'
## when loading 'CovRegFC'

library(cowplot)

## Loading required package: ggplot2

## 
## Attaching package: 'cowplot'

## The following object is masked from 'package:ggplot2':
## 
##     ggsave

library(ggplot2)
library(reshape2)
library(magrittr)
library(stringr)
library(devtools)
library(rstan)

## Loading required package: StanHeaders

## rstan (Version 2.15.1, packaged: 2017-04-19 05:03:57 UTC, GitRev: 2e1f913d3ca3)

## For execution on a local, multicore CPU with excess RAM we recommend calling
## rstan_options(auto_write = TRUE)
## options(mc.cores = parallel::detectCores())

## 
## Attaching package: 'rstan'

## The following object is masked from 'package:magrittr':
## 
##     extract

library(abind)
library(BatchJobs)

## Loading required package: BBmisc

## Sourced 1 configuration files:

##   1: /Users/cseiler/Code/RLibraries/BatchJobs/etc/BatchJobs_global_config.R

## BatchJobs configuration:
##   cluster functions: Interactive
##   mail.from: 
##   mail.to: 
##   mail.start: none
##   mail.done: none
##   mail.error: none
##   default.resources: 
##   debug: FALSE
##   raise.warnings: FALSE
##   staged.queries: TRUE
##   max.concurrent.jobs: Inf
##   fs.timeout: NA
##   measure.mem: TRUE

library(reshape2)
library(rethinking)

## Loading required package: parallel

## rethinking (Version 1.59)

params

## $num_regions
## [1] "15"

Load HCP Data

Analysis based on data from the WU-Minn HCP 1200 Subjects Data Release with four complete rfMRI runs (with 100% of collected time points). Load subject information and behavioral measurements.

subject_info = read.csv("HCP_PTN820/sample_info.csv")
subject_info = subject_info[,c(
  "Subject",
  "Age",
  "Gender",
  "Acquisition",
  "PSQI_AmtSleep", # amount of sleep in hours
  "PSQI_Score" # Pittsburgh Sleep Quality Index (PSQI) Completed
  )] 

Classification of sleep duration from Hirshkowitz et al. (2015)¹:

• short sleepers: average equal or less than 6 hours each night
• conventional sleepers: average between 7 and 9 hours each night

sleep_duration = rep("undefined",nrow(subject_info))
sleep_duration[subject_info$PSQI_AmtSleep <= 6] = "short"
sleep_duration[(subject_info$PSQI_AmtSleep >= 7) & (subject_info$PSQI_AmtSleep <= 9)] = "conventional"
subject_info$sleep_duration = factor(sleep_duration)
table(subject_info$sleep_duration)

## 
## conventional        short    undefined 
##          710          366          130

str(subject_info)

## 'data.frame':    1206 obs. of  7 variables:
##  $ Subject       : int  100004 100206 100307 100408 100610 101006 101107 101208 101309 101410 ...
##  $ Age           : Factor w/ 4 levels "22-25","26-30",..: 1 2 2 3 2 3 1 3 2 2 ...
##  $ Gender        : Factor w/ 2 levels "F","M": 2 2 1 2 2 1 2 1 2 2 ...
##  $ Acquisition   : Factor w/ 13 levels "Q01","Q02","Q03",..: 6 11 1 3 8 6 6 12 6 6 ...
##  $ PSQI_AmtSleep : num  5 6 6.5 8 7 7 4 7 5.5 5.5 ...
##  $ PSQI_Score    : int  8 6 4 4 4 2 9 6 5 8 ...
##  $ sleep_duration: Factor w/ 3 levels "conventional",..: 2 2 3 1 1 1 2 1 2 2 ...

Load timeseries.

num_regions = as.integer(params$num_regions)
num_regions

## [1] 15

channel_names = paste0("R",1:num_regions)
path = paste0("HCP_PTN820/node_timeseries/3T_HCP820_MSMAll_d",num_regions,"_ts2")
file_names = list.files(path = path,pattern = ".txt")
file_subject_ids = strsplit(file_names,split = ".txt") %>% unlist
ts = lapply(file_subject_ids,function(Subject) {
  print(paste("reading subject:",Subject))
  full_path = paste0(path,"/",Subject,".txt")
  timeseries = read.csv(full_path,header = FALSE,sep = " ")
  timeseries$Subject = Subject
  timeseries
}) %>% do.call(rbind,.) %>% data.frame

## [1] "reading subject: 100206"
## [1] "reading subject: 100307"
## [1] "reading subject: 100408"
## [1] "reading subject: 100610"
## [1] "reading subject: 101006"
## [1] "reading subject: 101107"
## [1] "reading subject: 101309"
## [1] "reading subject: 101915"
## [1] "reading subject: 102008"
## [1] "reading subject: 102311"
## [1] "reading subject: 102513"
## [1] "reading subject: 102816"
## [1] "reading subject: 103111"
## [1] "reading subject: 103414"
## [1] "reading subject: 103515"
## [1] "reading subject: 103818"
## [1] "reading subject: 104012"
## [1] "reading subject: 104416"
## [1] "reading subject: 104820"
## [1] "reading subject: 105014"
## [1] "reading subject: 105115"
## [1] "reading subject: 105216"
## [1] "reading subject: 105620"
## [1] "reading subject: 105923"
## [1] "reading subject: 106016"
## [1] "reading subject: 106319"
## [1] "reading subject: 106521"
## [1] "reading subject: 107018"
## [1] "reading subject: 107321"
## [1] "reading subject: 107422"
## [1] "reading subject: 107725"
## [1] "reading subject: 108121"
## [1] "reading subject: 108222"
## [1] "reading subject: 108323"
## [1] "reading subject: 108525"
## [1] "reading subject: 108828"
## [1] "reading subject: 109123"
## [1] "reading subject: 109325"
## [1] "reading subject: 109830"
## [1] "reading subject: 110007"
## [1] "reading subject: 110411"
## [1] "reading subject: 110613"
## [1] "reading subject: 111009"
## [1] "reading subject: 111312"
## [1] "reading subject: 111413"
## [1] "reading subject: 111514"
## [1] "reading subject: 111716"
## [1] "reading subject: 112112"
## [1] "reading subject: 112314"
## [1] "reading subject: 112516"
## [1] "reading subject: 112920"
## [1] "reading subject: 113215"
## [1] "reading subject: 113619"
## [1] "reading subject: 113922"
## [1] "reading subject: 114217"
## [1] "reading subject: 114318"
## [1] "reading subject: 114419"
## [1] "reading subject: 114621"
## [1] "reading subject: 114823"
## [1] "reading subject: 114924"
## [1] "reading subject: 115017"
## [1] "reading subject: 115219"
## [1] "reading subject: 115320"
## [1] "reading subject: 115825"
## [1] "reading subject: 116524"
## [1] "reading subject: 116726"
## [1] "reading subject: 117122"
## [1] "reading subject: 117324"
## [1] "reading subject: 117930"
## [1] "reading subject: 118023"
## [1] "reading subject: 118124"
## [1] "reading subject: 118225"
## [1] "reading subject: 118528"
## [1] "reading subject: 118730"
## [1] "reading subject: 118932"
## [1] "reading subject: 119126"
## [1] "reading subject: 120111"
## [1] "reading subject: 120212"
## [1] "reading subject: 120515"
## [1] "reading subject: 120717"
## [1] "reading subject: 121416"
## [1] "reading subject: 121618"
## [1] "reading subject: 121921"
## [1] "reading subject: 122317"
## [1] "reading subject: 122620"
## [1] "reading subject: 122822"
## [1] "reading subject: 123117"
## [1] "reading subject: 123420"
## [1] "reading subject: 123521"
## [1] "reading subject: 123824"
## [1] "reading subject: 123925"
## [1] "reading subject: 124220"
## [1] "reading subject: 124422"
## [1] "reading subject: 124624"
## [1] "reading subject: 124826"
## [1] "reading subject: 125525"
## [1] "reading subject: 126325"
## [1] "reading subject: 126628"
## [1] "reading subject: 127327"
## [1] "reading subject: 127630"
## [1] "reading subject: 127933"
## [1] "reading subject: 128026"
## [1] "reading subject: 128127"
## [1] "reading subject: 128632"
## [1] "reading subject: 128935"
## [1] "reading subject: 129028"
## [1] "reading subject: 129129"
## [1] "reading subject: 129331"
## [1] "reading subject: 129634"
## [1] "reading subject: 130013"
## [1] "reading subject: 130316"
## [1] "reading subject: 130417"
## [1] "reading subject: 130619"
## [1] "reading subject: 130821"
## [1] "reading subject: 130922"
## [1] "reading subject: 131217"
## [1] "reading subject: 131419"
## [1] "reading subject: 131722"
## [1] "reading subject: 131823"
## [1] "reading subject: 131924"
## [1] "reading subject: 132017"
## [1] "reading subject: 132118"
## [1] "reading subject: 133019"
## [1] "reading subject: 133625"
## [1] "reading subject: 133827"
## [1] "reading subject: 133928"
## [1] "reading subject: 134021"
## [1] "reading subject: 134223"
## [1] "reading subject: 134324"
## [1] "reading subject: 134425"
## [1] "reading subject: 134728"
## [1] "reading subject: 134829"
## [1] "reading subject: 135225"
## [1] "reading subject: 135528"
## [1] "reading subject: 135730"
## [1] "reading subject: 135932"
## [1] "reading subject: 136227"
## [1] "reading subject: 136732"
## [1] "reading subject: 136833"
## [1] "reading subject: 137027"
## [1] "reading subject: 137128"
## [1] "reading subject: 137229"
## [1] "reading subject: 137633"
## [1] "reading subject: 137936"
## [1] "reading subject: 138231"
## [1] "reading subject: 138534"
## [1] "reading subject: 138837"
## [1] "reading subject: 139233"
## [1] "reading subject: 139637"
## [1] "reading subject: 139839"
## [1] "reading subject: 140117"
## [1] "reading subject: 140319"
## [1] "reading subject: 140824"
## [1] "reading subject: 140925"
## [1] "reading subject: 141119"
## [1] "reading subject: 141422"
## [1] "reading subject: 141826"
## [1] "reading subject: 142828"
## [1] "reading subject: 143325"
## [1] "reading subject: 143426"
## [1] "reading subject: 144125"
## [1] "reading subject: 144226"
## [1] "reading subject: 144428"
## [1] "reading subject: 144731"
## [1] "reading subject: 144832"
## [1] "reading subject: 145127"
## [1] "reading subject: 145834"
## [1] "reading subject: 146129"
## [1] "reading subject: 146331"
## [1] "reading subject: 146432"
## [1] "reading subject: 146533"
## [1] "reading subject: 146937"
## [1] "reading subject: 147030"
## [1] "reading subject: 147737"
## [1] "reading subject: 148032"
## [1] "reading subject: 148133"
## [1] "reading subject: 148335"
## [1] "reading subject: 148436"
## [1] "reading subject: 148840"
## [1] "reading subject: 148941"
## [1] "reading subject: 149236"
## [1] "reading subject: 149337"
## [1] "reading subject: 149539"
## [1] "reading subject: 149741"
## [1] "reading subject: 149842"
## [1] "reading subject: 150625"
## [1] "reading subject: 150726"
## [1] "reading subject: 150928"
## [1] "reading subject: 151223"
## [1] "reading subject: 151425"
## [1] "reading subject: 151526"
## [1] "reading subject: 151627"
## [1] "reading subject: 151728"
## [1] "reading subject: 151829"
## [1] "reading subject: 152831"
## [1] "reading subject: 153025"
## [1] "reading subject: 153227"
## [1] "reading subject: 153429"
## [1] "reading subject: 153631"
## [1] "reading subject: 153732"
## [1] "reading subject: 153833"
## [1] "reading subject: 154229"
## [1] "reading subject: 154431"
## [1] "reading subject: 154532"
## [1] "reading subject: 154734"
## [1] "reading subject: 154835"
## [1] "reading subject: 154936"
## [1] "reading subject: 155635"
## [1] "reading subject: 155938"
## [1] "reading subject: 156031"
## [1] "reading subject: 156233"
## [1] "reading subject: 156334"
## [1] "reading subject: 156435"
## [1] "reading subject: 156536"
## [1] "reading subject: 156637"
## [1] "reading subject: 157336"
## [1] "reading subject: 157437"
## [1] "reading subject: 157942"
## [1] "reading subject: 158035"
## [1] "reading subject: 158136"
## [1] "reading subject: 158338"
## [1] "reading subject: 158540"
## [1] "reading subject: 158843"
## [1] "reading subject: 159138"
## [1] "reading subject: 159239"
## [1] "reading subject: 159340"
## [1] "reading subject: 159441"
## [1] "reading subject: 159744"
## [1] "reading subject: 160123"
## [1] "reading subject: 160729"
## [1] "reading subject: 160830"
## [1] "reading subject: 161327"
## [1] "reading subject: 161630"
## [1] "reading subject: 161731"
## [1] "reading subject: 162026"
## [1] "reading subject: 162228"
## [1] "reading subject: 162329"
## [1] "reading subject: 162733"
## [1] "reading subject: 162935"
## [1] "reading subject: 163129"
## [1] "reading subject: 163331"
## [1] "reading subject: 163432"
## [1] "reading subject: 163836"
## [1] "reading subject: 164030"
## [1] "reading subject: 164131"
## [1] "reading subject: 164636"
## [1] "reading subject: 164939"
## [1] "reading subject: 165032"
## [1] "reading subject: 165638"
## [1] "reading subject: 165840"
## [1] "reading subject: 166438"
## [1] "reading subject: 166640"
## [1] "reading subject: 167036"
## [1] "reading subject: 167238"
## [1] "reading subject: 167743"
## [1] "reading subject: 168139"
## [1] "reading subject: 168240"
## [1] "reading subject: 168341"
## [1] "reading subject: 168745"
## [1] "reading subject: 169343"
## [1] "reading subject: 169444"
## [1] "reading subject: 169949"
## [1] "reading subject: 170631"
## [1] "reading subject: 171330"
## [1] "reading subject: 171532"
## [1] "reading subject: 171633"
## [1] "reading subject: 172029"
## [1] "reading subject: 172130"
## [1] "reading subject: 172332"
## [1] "reading subject: 172433"
## [1] "reading subject: 172534"
## [1] "reading subject: 172938"
## [1] "reading subject: 173334"
## [1] "reading subject: 173435"
## [1] "reading subject: 173536"
## [1] "reading subject: 173637"
## [1] "reading subject: 173738"
## [1] "reading subject: 173839"
## [1] "reading subject: 173940"
## [1] "reading subject: 174437"
## [1] "reading subject: 174841"
## [1] "reading subject: 175035"
## [1] "reading subject: 175237"
## [1] "reading subject: 175338"
## [1] "reading subject: 175439"
## [1] "reading subject: 175540"
## [1] "reading subject: 175742"
## [1] "reading subject: 176037"
## [1] "reading subject: 176239"
## [1] "reading subject: 176441"
## [1] "reading subject: 176542"
## [1] "reading subject: 176744"
## [1] "reading subject: 177241"
## [1] "reading subject: 177645"
## [1] "reading subject: 177746"
## [1] "reading subject: 178142"
## [1] "reading subject: 178243"
## [1] "reading subject: 178647"
## [1] "reading subject: 178748"
## [1] "reading subject: 178849"
## [1] "reading subject: 178950"
## [1] "reading subject: 179245"
## [1] "reading subject: 179346"
## [1] "reading subject: 180129"
## [1] "reading subject: 180432"
## [1] "reading subject: 180735"
## [1] "reading subject: 180836"
## [1] "reading subject: 180937"
## [1] "reading subject: 181131"
## [1] "reading subject: 181232"
## [1] "reading subject: 181636"
## [1] "reading subject: 182032"
## [1] "reading subject: 182436"
## [1] "reading subject: 182739"
## [1] "reading subject: 182840"
## [1] "reading subject: 183034"
## [1] "reading subject: 185139"
## [1] "reading subject: 185341"
## [1] "reading subject: 185442"
## [1] "reading subject: 185846"
## [1] "reading subject: 185947"
## [1] "reading subject: 186141"
## [1] "reading subject: 186444"
## [1] "reading subject: 187143"
## [1] "reading subject: 187345"
## [1] "reading subject: 187547"
## [1] "reading subject: 187850"
## [1] "reading subject: 188347"
## [1] "reading subject: 188448"
## [1] "reading subject: 188549"
## [1] "reading subject: 188751"
## [1] "reading subject: 189349"
## [1] "reading subject: 189450"
## [1] "reading subject: 190031"
## [1] "reading subject: 191033"
## [1] "reading subject: 191336"
## [1] "reading subject: 191437"
## [1] "reading subject: 191841"
## [1] "reading subject: 191942"
## [1] "reading subject: 192035"
## [1] "reading subject: 192136"
## [1] "reading subject: 192439"
## [1] "reading subject: 192540"
## [1] "reading subject: 192641"
## [1] "reading subject: 192843"
## [1] "reading subject: 193239"
## [1] "reading subject: 194140"
## [1] "reading subject: 194645"
## [1] "reading subject: 194746"
## [1] "reading subject: 194847"
## [1] "reading subject: 195041"
## [1] "reading subject: 195445"
## [1] "reading subject: 195647"
## [1] "reading subject: 195849"
## [1] "reading subject: 195950"
## [1] "reading subject: 196144"
## [1] "reading subject: 196346"
## [1] "reading subject: 196750"
## [1] "reading subject: 197348"
## [1] "reading subject: 197550"
## [1] "reading subject: 198249"
## [1] "reading subject: 198350"
## [1] "reading subject: 198451"
## [1] "reading subject: 198653"
## [1] "reading subject: 198855"
## [1] "reading subject: 199150"
## [1] "reading subject: 199251"
## [1] "reading subject: 199453"
## [1] "reading subject: 199655"
## [1] "reading subject: 199958"
## [1] "reading subject: 200008"
## [1] "reading subject: 200109"
## [1] "reading subject: 200311"
## [1] "reading subject: 200614"
## [1] "reading subject: 200917"
## [1] "reading subject: 201111"
## [1] "reading subject: 201414"
## [1] "reading subject: 201515"
## [1] "reading subject: 201818"
## [1] "reading subject: 202113"
## [1] "reading subject: 202719"
## [1] "reading subject: 203418"
## [1] "reading subject: 203923"
## [1] "reading subject: 204016"
## [1] "reading subject: 204319"
## [1] "reading subject: 204420"
## [1] "reading subject: 204521"
## [1] "reading subject: 204622"
## [1] "reading subject: 205119"
## [1] "reading subject: 205220"
## [1] "reading subject: 205725"
## [1] "reading subject: 205826"
## [1] "reading subject: 206222"
## [1] "reading subject: 207123"
## [1] "reading subject: 207426"
## [1] "reading subject: 208024"
## [1] "reading subject: 208125"
## [1] "reading subject: 208226"
## [1] "reading subject: 208327"
## [1] "reading subject: 209127"
## [1] "reading subject: 209228"
## [1] "reading subject: 209329"
## [1] "reading subject: 209834"
## [1] "reading subject: 209935"
## [1] "reading subject: 210011"
## [1] "reading subject: 210415"
## [1] "reading subject: 210617"
## [1] "reading subject: 211114"
## [1] "reading subject: 211215"
## [1] "reading subject: 211316"
## [1] "reading subject: 211417"
## [1] "reading subject: 211720"
## [1] "reading subject: 211922"
## [1] "reading subject: 212015"
## [1] "reading subject: 212116"
## [1] "reading subject: 212217"
## [1] "reading subject: 212318"
## [1] "reading subject: 212419"
## [1] "reading subject: 212823"
## [1] "reading subject: 213421"
## [1] "reading subject: 214019"
## [1] "reading subject: 214221"
## [1] "reading subject: 214423"
## [1] "reading subject: 214524"
## [1] "reading subject: 214625"
## [1] "reading subject: 214726"
## [1] "reading subject: 217126"
## [1] "reading subject: 217429"
## [1] "reading subject: 220721"
## [1] "reading subject: 221319"
## [1] "reading subject: 223929"
## [1] "reading subject: 224022"
## [1] "reading subject: 227432"
## [1] "reading subject: 228434"
## [1] "reading subject: 231928"
## [1] "reading subject: 233326"
## [1] "reading subject: 236130"
## [1] "reading subject: 237334"
## [1] "reading subject: 239944"
## [1] "reading subject: 245333"
## [1] "reading subject: 246133"
## [1] "reading subject: 248339"
## [1] "reading subject: 249947"
## [1] "reading subject: 250427"
## [1] "reading subject: 250932"
## [1] "reading subject: 251833"
## [1] "reading subject: 255639"
## [1] "reading subject: 256540"
## [1] "reading subject: 257542"
## [1] "reading subject: 257845"
## [1] "reading subject: 263436"
## [1] "reading subject: 268749"
## [1] "reading subject: 268850"
## [1] "reading subject: 270332"
## [1] "reading subject: 275645"
## [1] "reading subject: 280739"
## [1] "reading subject: 280941"
## [1] "reading subject: 283543"
## [1] "reading subject: 285345"
## [1] "reading subject: 285446"
## [1] "reading subject: 286650"
## [1] "reading subject: 287248"
## [1] "reading subject: 289555"
## [1] "reading subject: 290136"
## [1] "reading subject: 293748"
## [1] "reading subject: 295146"
## [1] "reading subject: 297655"
## [1] "reading subject: 298051"
## [1] "reading subject: 298455"
## [1] "reading subject: 299154"
## [1] "reading subject: 300618"
## [1] "reading subject: 303119"
## [1] "reading subject: 303624"
## [1] "reading subject: 304020"
## [1] "reading subject: 304727"
## [1] "reading subject: 305830"
## [1] "reading subject: 307127"
## [1] "reading subject: 308129"
## [1] "reading subject: 308331"
## [1] "reading subject: 309636"
## [1] "reading subject: 310621"
## [1] "reading subject: 311320"
## [1] "reading subject: 316633"
## [1] "reading subject: 316835"
## [1] "reading subject: 318637"
## [1] "reading subject: 320826"
## [1] "reading subject: 321323"
## [1] "reading subject: 322224"
## [1] "reading subject: 329440"
## [1] "reading subject: 330324"
## [1] "reading subject: 333330"
## [1] "reading subject: 334635"
## [1] "reading subject: 336841"
## [1] "reading subject: 339847"
## [1] "reading subject: 341834"
## [1] "reading subject: 346137"
## [1] "reading subject: 346945"
## [1] "reading subject: 348545"
## [1] "reading subject: 352132"
## [1] "reading subject: 352738"
## [1] "reading subject: 353740"
## [1] "reading subject: 356948"
## [1] "reading subject: 358144"
## [1] "reading subject: 361234"
## [1] "reading subject: 361941"
## [1] "reading subject: 365343"
## [1] "reading subject: 366042"
## [1] "reading subject: 366446"
## [1] "reading subject: 371843"
## [1] "reading subject: 377451"
## [1] "reading subject: 378857"
## [1] "reading subject: 379657"
## [1] "reading subject: 380036"
## [1] "reading subject: 381038"
## [1] "reading subject: 381543"
## [1] "reading subject: 382242"
## [1] "reading subject: 385450"
## [1] "reading subject: 386250"
## [1] "reading subject: 387959"
## [1] "reading subject: 389357"
## [1] "reading subject: 390645"
## [1] "reading subject: 391748"
## [1] "reading subject: 393247"
## [1] "reading subject: 393550"
## [1] "reading subject: 395251"
## [1] "reading subject: 395756"
## [1] "reading subject: 395958"
## [1] "reading subject: 397154"
## [1] "reading subject: 397760"
## [1] "reading subject: 397861"
## [1] "reading subject: 406432"
## [1] "reading subject: 406836"
## [1] "reading subject: 412528"
## [1] "reading subject: 414229"
## [1] "reading subject: 415837"
## [1] "reading subject: 422632"
## [1] "reading subject: 424939"
## [1] "reading subject: 429040"
## [1] "reading subject: 432332"
## [1] "reading subject: 433839"
## [1] "reading subject: 436239"
## [1] "reading subject: 436845"
## [1] "reading subject: 441939"
## [1] "reading subject: 445543"
## [1] "reading subject: 448347"
## [1] "reading subject: 449753"
## [1] "reading subject: 453441"
## [1] "reading subject: 456346"
## [1] "reading subject: 459453"
## [1] "reading subject: 465852"
## [1] "reading subject: 467351"
## [1] "reading subject: 475855"
## [1] "reading subject: 479762"
## [1] "reading subject: 480141"
## [1] "reading subject: 481951"
## [1] "reading subject: 485757"
## [1] "reading subject: 486759"
## [1] "reading subject: 495255"
## [1] "reading subject: 497865"
## [1] "reading subject: 499566"
## [1] "reading subject: 500222"
## [1] "reading subject: 506234"
## [1] "reading subject: 510326"
## [1] "reading subject: 512835"
## [1] "reading subject: 513736"
## [1] "reading subject: 517239"
## [1] "reading subject: 519950"
## [1] "reading subject: 520228"
## [1] "reading subject: 522434"
## [1] "reading subject: 523032"
## [1] "reading subject: 524135"
## [1] "reading subject: 525541"
## [1] "reading subject: 529549"
## [1] "reading subject: 529953"
## [1] "reading subject: 530635"
## [1] "reading subject: 531536"
## [1] "reading subject: 536647"
## [1] "reading subject: 540436"
## [1] "reading subject: 541943"
## [1] "reading subject: 545345"
## [1] "reading subject: 547046"
## [1] "reading subject: 548250"
## [1] "reading subject: 552544"
## [1] "reading subject: 553344"
## [1] "reading subject: 555348"
## [1] "reading subject: 555651"
## [1] "reading subject: 557857"
## [1] "reading subject: 559053"
## [1] "reading subject: 561242"
## [1] "reading subject: 561444"
## [1] "reading subject: 562345"
## [1] "reading subject: 562446"
## [1] "reading subject: 565452"
## [1] "reading subject: 566454"
## [1] "reading subject: 567052"
## [1] "reading subject: 567961"
## [1] "reading subject: 568963"
## [1] "reading subject: 570243"
## [1] "reading subject: 571144"
## [1] "reading subject: 572045"
## [1] "reading subject: 573249"
## [1] "reading subject: 573451"
## [1] "reading subject: 576255"
## [1] "reading subject: 579665"
## [1] "reading subject: 579867"
## [1] "reading subject: 580044"
## [1] "reading subject: 580347"
## [1] "reading subject: 580650"
## [1] "reading subject: 580751"
## [1] "reading subject: 581349"
## [1] "reading subject: 581450"
## [1] "reading subject: 583858"
## [1] "reading subject: 585256"
## [1] "reading subject: 585862"
## [1] "reading subject: 586460"
## [1] "reading subject: 587664"
## [1] "reading subject: 588565"
## [1] "reading subject: 592455"
## [1] "reading subject: 594156"
## [1] "reading subject: 597869"
## [1] "reading subject: 598568"
## [1] "reading subject: 599065"
## [1] "reading subject: 599469"
## [1] "reading subject: 599671"
## [1] "reading subject: 601127"
## [1] "reading subject: 604537"
## [1] "reading subject: 609143"
## [1] "reading subject: 611938"
## [1] "reading subject: 613538"
## [1] "reading subject: 614439"
## [1] "reading subject: 615744"
## [1] "reading subject: 616645"
## [1] "reading subject: 617748"
## [1] "reading subject: 618952"
## [1] "reading subject: 620434"
## [1] "reading subject: 622236"
## [1] "reading subject: 623844"
## [1] "reading subject: 626648"
## [1] "reading subject: 627549"
## [1] "reading subject: 627852"
## [1] "reading subject: 628248"
## [1] "reading subject: 633847"
## [1] "reading subject: 638049"
## [1] "reading subject: 645450"
## [1] "reading subject: 645551"
## [1] "reading subject: 647858"
## [1] "reading subject: 654350"
## [1] "reading subject: 654754"
## [1] "reading subject: 656253"
## [1] "reading subject: 656657"
## [1] "reading subject: 657659"
## [1] "reading subject: 660951"
## [1] "reading subject: 662551"
## [1] "reading subject: 663755"
## [1] "reading subject: 664757"
## [1] "reading subject: 665254"
## [1] "reading subject: 667056"
## [1] "reading subject: 668361"
## [1] "reading subject: 671855"
## [1] "reading subject: 672756"
## [1] "reading subject: 673455"
## [1] "reading subject: 677766"
## [1] "reading subject: 677968"
## [1] "reading subject: 679568"
## [1] "reading subject: 679770"
## [1] "reading subject: 680250"
## [1] "reading subject: 680957"
## [1] "reading subject: 683256"
## [1] "reading subject: 685058"
## [1] "reading subject: 686969"
## [1] "reading subject: 687163"
## [1] "reading subject: 690152"
## [1] "reading subject: 693764"
## [1] "reading subject: 695768"
## [1] "reading subject: 700634"
## [1] "reading subject: 702133"
## [1] "reading subject: 704238"
## [1] "reading subject: 705341"
## [1] "reading subject: 706040"
## [1] "reading subject: 707749"
## [1] "reading subject: 709551"
## [1] "reading subject: 713239"
## [1] "reading subject: 715041"
## [1] "reading subject: 715647"
## [1] "reading subject: 715950"
## [1] "reading subject: 720337"
## [1] "reading subject: 724446"
## [1] "reading subject: 725751"
## [1] "reading subject: 727553"
## [1] "reading subject: 727654"
## [1] "reading subject: 729254"
## [1] "reading subject: 729557"
## [1] "reading subject: 731140"
## [1] "reading subject: 732243"
## [1] "reading subject: 734045"
## [1] "reading subject: 735148"
## [1] "reading subject: 737960"
## [1] "reading subject: 742549"
## [1] "reading subject: 744553"
## [1] "reading subject: 748258"
## [1] "reading subject: 748662"
## [1] "reading subject: 749058"
## [1] "reading subject: 749361"
## [1] "reading subject: 751348"
## [1] "reading subject: 753150"
## [1] "reading subject: 753251"
## [1] "reading subject: 756055"
## [1] "reading subject: 759869"
## [1] "reading subject: 761957"
## [1] "reading subject: 765056"
## [1] "reading subject: 767464"
## [1] "reading subject: 769064"
## [1] "reading subject: 770352"
## [1] "reading subject: 771354"
## [1] "reading subject: 773257"
## [1] "reading subject: 779370"
## [1] "reading subject: 782561"
## [1] "reading subject: 783462"
## [1] "reading subject: 784565"
## [1] "reading subject: 788876"
## [1] "reading subject: 789373"
## [1] "reading subject: 792564"
## [1] "reading subject: 792766"
## [1] "reading subject: 792867"
## [1] "reading subject: 793465"
## [1] "reading subject: 800941"
## [1] "reading subject: 802844"
## [1] "reading subject: 803240"
## [1] "reading subject: 810843"
## [1] "reading subject: 812746"
## [1] "reading subject: 814649"
## [1] "reading subject: 816653"
## [1] "reading subject: 818859"
## [1] "reading subject: 820745"
## [1] "reading subject: 825048"
## [1] "reading subject: 826353"
## [1] "reading subject: 826454"
## [1] "reading subject: 833148"
## [1] "reading subject: 833249"
## [1] "reading subject: 835657"
## [1] "reading subject: 837560"
## [1] "reading subject: 837964"
## [1] "reading subject: 841349"
## [1] "reading subject: 843151"
## [1] "reading subject: 844961"
## [1] "reading subject: 845458"
## [1] "reading subject: 849264"
## [1] "reading subject: 849971"
## [1] "reading subject: 852455"
## [1] "reading subject: 856766"
## [1] "reading subject: 856968"
## [1] "reading subject: 857263"
## [1] "reading subject: 859671"
## [1] "reading subject: 861456"
## [1] "reading subject: 865363"
## [1] "reading subject: 867468"
## [1] "reading subject: 870861"
## [1] "reading subject: 871762"
## [1] "reading subject: 871964"
## [1] "reading subject: 872158"
## [1] "reading subject: 872562"
## [1] "reading subject: 872764"
## [1] "reading subject: 873968"
## [1] "reading subject: 877168"
## [1] "reading subject: 877269"
## [1] "reading subject: 878776"
## [1] "reading subject: 880157"
## [1] "reading subject: 882161"
## [1] "reading subject: 885975"
## [1] "reading subject: 887373"
## [1] "reading subject: 889579"
## [1] "reading subject: 891667"
## [1] "reading subject: 894067"
## [1] "reading subject: 894673"
## [1] "reading subject: 894774"
## [1] "reading subject: 896778"
## [1] "reading subject: 896879"
## [1] "reading subject: 898176"
## [1] "reading subject: 899885"
## [1] "reading subject: 901038"
## [1] "reading subject: 901139"
## [1] "reading subject: 901442"
## [1] "reading subject: 904044"
## [1] "reading subject: 907656"
## [1] "reading subject: 910241"
## [1] "reading subject: 910443"
## [1] "reading subject: 912447"
## [1] "reading subject: 917255"
## [1] "reading subject: 917558"
## [1] "reading subject: 919966"
## [1] "reading subject: 922854"
## [1] "reading subject: 923755"
## [1] "reading subject: 926862"
## [1] "reading subject: 927359"
## [1] "reading subject: 930449"
## [1] "reading subject: 932554"
## [1] "reading subject: 937160"
## [1] "reading subject: 942658"
## [1] "reading subject: 947668"
## [1] "reading subject: 951457"
## [1] "reading subject: 952863"
## [1] "reading subject: 955465"
## [1] "reading subject: 957974"
## [1] "reading subject: 958976"
## [1] "reading subject: 959574"
## [1] "reading subject: 965367"
## [1] "reading subject: 965771"
## [1] "reading subject: 966975"
## [1] "reading subject: 978578"
## [1] "reading subject: 979984"
## [1] "reading subject: 983773"
## [1] "reading subject: 984472"
## [1] "reading subject: 987983"
## [1] "reading subject: 990366"
## [1] "reading subject: 991267"
## [1] "reading subject: 992673"
## [1] "reading subject: 992774"
## [1] "reading subject: 993675"
## [1] "reading subject: 994273"
## [1] "reading subject: 996782"

names(ts)[1:num_regions] = channel_names

Merge timeseries and subject info into one data frame.

ts_subject_info = merge(ts,subject_info,by = "Subject")

Plot timeseries for subject 100206.

timeseries = subset(ts_subject_info,Subject == 100206)
timeseries$timepoint = rep(1:(nrow(timeseries)/4),4)
timeseries$run = lapply(paste0("run_",1:4),function(i) rep(i,nrow(timeseries)/4)) %>% unlist
timeseries_long = melt(timeseries,id.vars = c("timepoint","run",names(subject_info)))
ggplot(timeseries_long,aes(x = timepoint,y = value,color = variable)) + 
  geom_line() + facet_wrap(~run)

timeseries_long_subset = subset(timeseries_long,
                                timepoint < 51 & (variable == "R1" | variable == "R2" | variable == "R3"))
ggplot(timeseries_long_subset,aes(x = timepoint,y = value,color = variable)) + 
  geom_line() + facet_wrap(~run)

Take a sample of 20 subjects.

# set.seed(934321)
# ts_subject_info = subset(ts_subject_info,sleep_duration == "short" | sleep_duration == "conventional")
# subsample_size = 20
# sampled_subjects = sample(names(table(ts_subject_info$Subject)),size = subsample_size)
# ts_subject_info_subset = subset(ts_subject_info,Subject %in% sampled_subjects) %>% droplevels
# str(ts_subject_info_subset)
# table(ts_subject_info_subset$Subject)
# table(ts_subject_info_subset$sleep_duration)/table(ts_subject_info_subset$Subject)[1]

Compute average covariance matrices per task/condition.

subset_tasks = function(data,factor_name,level_name) {
  data_subgroup = data[which(data[factor_name] == level_name),] %>% droplevels
  region_ids = str_detect(names(ts_subject_info),"R") %>% which
  subject_ids = table(data_subgroup$Subject) %>% names
  scov = list()
  K = length(subject_ids)
  P = length(region_ids)
  scov$YY = array(dim = c(K,P,P))
  for(i in 1:K) {
    Y = subset(data_subgroup,Subject == subject_ids[i])[,region_ids]
    Y_c = scale(Y,center = TRUE,scale = FALSE)
    scov$YY[i,,] = t(Y_c) %*% Y_c
  }
  scov$N = table(data_subgroup$Subject)
  scov
}
#scov_short = subset_tasks(ts_subject_info_subset,"sleep_duration","short")
scov_short = subset_tasks(ts_subject_info,"sleep_duration","short")
#scov_conventional = subset_tasks(ts_subject_info_subset,"sleep_duration","conventional")
scov_conventional = subset_tasks(ts_subject_info,"sleep_duration","conventional")

Plot average covariance matrices per task/condition.

plot_cov = function(YY,limit_value) {
  rownames(YY) = colnames(YY) = channel_names
  get_upper_tri = function(mat) {
    mat[lower.tri(mat)] = NA
    mat
  }
  YY_long = melt(get_upper_tri(YY), na.rm = TRUE)
  ggplot(data = YY_long, aes(Var2, Var1, fill = value)) +
    geom_tile(color = "white") +
    scale_fill_gradient2(low = "blue", high = "red", mid = "white",
                         midpoint = 0, space = "Lab",
                         limit = c(-limit_value,limit_value),
                         name = "Covariance") +
    theme_minimal() +
    theme(axis.text.x = element_text(angle = 90, vjust = 1, hjust = 1),
          axis.title.x=element_blank(),
          axis.title.y=element_blank(),
          legend.position="right") +
    coord_fixed()
}
scov_short_mean = apply(scov_short$YY,MARGIN = c(2,3),FUN = mean)
scov_conventional_mean = apply(scov_conventional$YY,MARGIN = c(2,3),FUN = mean)
limit_value = max(abs(c(scov_short_mean[lower.tri(scov_short_mean)] %>% abs %>% max,
                        scov_conventional_mean[lower.tri(scov_conventional_mean)] %>% abs %>% max)))
plot_cov(scov_short_mean,limit_value)

plot_cov(scov_conventional_mean,limit_value)

Full Covariance Modeling

YY = abind(scov_conventional$YY,scov_short$YY,along = 1)
N = min(scov_conventional$N,scov_short$N)
short = c(rep(0,length(scov_conventional$N)),rep(1,length(scov_short$N)))
X = data.frame(intercept=rep(1,length(short)),short)
# run it in batch mode
log_name = paste0("HCP_full_",num_regions)
if(dir.exists(log_name)) unlink(log_name,recursive = TRUE) # cleanup
reg = makeRegistry(id = log_name, file.dir = log_name)

## Creating dir: HCP_full_15

## Warning in rsqlite_fetch(res@ptr, n = n): Don't need to call dbFetch() for
## statements, only for queries

## Warning in rsqlite_fetch(res@ptr, n = n): Don't need to call dbFetch() for
## statements, only for queries

## Warning in rsqlite_fetch(res@ptr, n = n): Don't need to call dbFetch() for
## statements, only for queries

## Warning in rsqlite_fetch(res@ptr, n = n): Don't need to call dbFetch() for
## statements, only for queries

## Saving registry: /Users/cseiler/CovRegFC_HCP/HCP_full_15/registry.RData

batchExport(reg, li = list(YY = YY,
                           X = X,
                           N = N,
                           scov_short_mean = scov_short_mean,
                           scov_conventional_mean = scov_conventional_mean),
            overwrite = TRUE)
fun = function(seed) {
  library(CovRegFC)
  CovRegFC::covreg_full(YY = YY,
                        X = X,
                        N = N,
                        scov_short_mean = scov_short_mean,
                        scov_conventional_mean = scov_conventional_mean,
                        seed = seed)
}
batchMap(reg, fun, 15324)

## Adding 1 jobs to DB.

## Warning in checkUserFunction(fun): The environment of provided function
## exceeds 10Mb.

## Warning in checkUserFunction(fun): Don't need to call dbFetch() for
## statements, only for queries

## Warning in checkUserFunction(fun): Don't need to call dbFetch() for
## statements, only for queries

submitJobs(reg)

## Saving conf: /Users/cseiler/CovRegFC_HCP/HCP_full_15/conf.RData

## Submitting 1 chunks / 1 jobs.

## Cluster functions: Interactive.

## Auto-mailer settings: start=none, done=none, error=none.

## Loading registry: /Users/cseiler/CovRegFC_HCP/HCP_full_15/registry.RData

## Loading RData files: N,scov_conventional_mean,scov_short_mean,X,YY

## Loading conf: /Users/cseiler/CovRegFC_HCP/HCP_full_15/conf.RData

## 2017-06-22 17:11:33: Starting job on node Christofs-MacBook-Pro-2.local.

## Auto-mailer settings: start=none, done=none, error=none.

## Setting work dir: /Users/cseiler/CovRegFC_HCP

## Loading jobs from file '/Users/cseiler/CovRegFC_HCP/HCP_full_15/jobs/01/1.rds'

## ########## Executing jid=1 ##########

## Timestamp: 2017-06-22 17:11:33

## Setting seed: 281998234

## Writing result file: /Users/cseiler/CovRegFC_HCP/HCP_full_15/jobs/01/1-result.RData

## 2017-06-22 18:19:13: All done.

## Setting work back to: /Users/cseiler/CovRegFC_HCP

## Memory usage according to gc:

## Sending 1 submit messages...
## Might take some time, do not interrupt this!

waitForJobs(reg,sleep = 360)

## Syncing registry ...

## Warning in rsqlite_fetch(res@ptr, n = n): Don't need to call dbFetch() for
## statements, only for queries

## Warning in rsqlite_fetch(res@ptr, n = n): Don't need to call dbFetch() for
## statements, only for queries

## Warning in rsqlite_fetch(res@ptr, n = n): Don't need to call dbFetch() for
## statements, only for queries

## Warning in rsqlite_fetch(res@ptr, n = n): Don't need to call dbFetch() for
## statements, only for queries

## [1] TRUE

res_list = reduceResultsList(reg)

## Reducing 1 results...

fit = res_list[[1]]
save(fit,file = paste0("fit_full_",num_regions,".Rdata"))

Convergence diagnostics.

CovRegFC::plot_diagnostics(fit,"sigma",num_par = 4)

CovRegFC::plot_diagnostics(fit,"Omega",num_par = 4)

CovRegFC::plot_diagnostics(fit,"nu",num_par = 4)

Plot differential correlation.

p_cov_da = CovRegFC::plot_post_mean(fit,channel_names)
p_cov_da

ggsave(p_cov_da,filename = paste0("cov_short_vs_conventional_",num_regions,".pdf"),height = 4,width = 6)

Differential analysis: shorter correlation stronger than conventional sleep duration correlations?

CovRegFC::plot_corr_da(fit,channel_names,"Differential correlations in short vs. conventional")

ggsave(filename = paste0("full_connectivity_",num_regions,".pdf"),height = 4,width = 5)

Simulations from Correlation Matrix Prior

Distribution of \(1000\) off-diagonal elements \(R_{12}\) extracted from correlation matrix drawn from LKJ prior.

Session Info

session_info()

## Session info -------------------------------------------------------------

##  setting  value                       
##  version  R version 3.3.2 (2016-10-31)
##  system   x86_64, darwin13.4.0        
##  ui       X11                         
##  language (EN)                        
##  collate  en_US.UTF-8                 
##  tz       America/Los_Angeles         
##  date     2017-06-22

## Packages -----------------------------------------------------------------

##  package      * version    date      
##  abind        * 1.4-5      2016-07-21
##  backports      1.1.0      2017-05-22
##  base         * 3.3.2      2016-10-31
##  base64enc      0.1-3      2015-07-28
##  BatchJobs    * 1.7        2016-10-01
##  BBmisc       * 1.11       2017-03-10
##  bit            1.1-12     2014-04-09
##  bit64          0.9-7      2017-05-08
##  blob           1.1.0      2017-06-17
##  brew           1.0-6      2011-04-13
##  checkmate      1.8.2      2016-11-02
##  coda           0.19-1     2016-12-08
##  codetools      0.2-15     2016-10-05
##  colorspace     1.3-2      2016-12-14
##  CovRegFC     * 0.1.0      2017-06-22
##  cowplot      * 0.7.0      2016-10-28
##  data.table     1.10.4     2017-02-01
##  datasets     * 3.3.2      2016-10-31
##  DBI            0.7        2017-06-18
##  devtools     * 1.13.2     2017-06-02
##  digest         0.6.12     2017-01-27
##  evaluate       0.10       2016-10-11
##  ggplot2      * 2.2.1      2016-12-30
##  graphics     * 3.3.2      2016-10-31
##  grDevices    * 3.3.2      2016-10-31
##  grid           3.3.2      2016-10-31
##  gridExtra      2.2.1      2016-02-29
##  gtable         0.2.0      2016-02-26
##  htmltools      0.3.6      2017-04-28
##  inline         0.3.14     2015-04-13
##  knitr          1.16       2017-05-18
##  labeling       0.3        2014-08-23
##  lattice        0.20-35    2017-03-25
##  lazyeval       0.2.0.9000 2017-05-09
##  loo            1.1.0      2017-03-27
##  magrittr     * 1.5        2014-11-22
##  MASS           7.3-47     2017-04-21
##  matrixStats    0.52.2     2017-04-14
##  memoise        1.1.0      2017-04-21
##  methods      * 3.3.2      2016-10-31
##  munsell        0.4.3      2016-02-13
##  mvtnorm        1.0-6      2017-03-02
##  parallel     * 3.3.2      2016-10-31
##  pkgconfig      2.0.1      2017-03-21
##  plyr           1.8.4      2016-06-08
##  RColorBrewer   1.1-2      2014-12-07
##  Rcpp         * 0.12.11    2017-05-22
##  reshape2     * 1.4.2      2016-10-22
##  rethinking   * 1.59       2017-02-01
##  rlang          0.1.1      2017-05-18
##  rmarkdown      1.6        2017-06-15
##  rprojroot      1.2        2017-01-16
##  RSQLite        2.0        2017-06-19
##  rstan        * 2.15.1     2017-04-19
##  scales         0.4.1      2016-11-09
##  sendmailR      1.2-1      2014-09-21
##  StanHeaders  * 2.15.0-1   2017-04-19
##  stats        * 3.3.2      2016-10-31
##  stats4         3.3.2      2016-10-31
##  stringi        1.1.5      2017-04-07
##  stringr      * 1.2.0      2017-02-18
##  tibble         1.3.3      2017-05-28
##  tools          3.3.2      2016-10-31
##  utils        * 3.3.2      2016-10-31
##  withr          1.0.2      2016-06-20
##  yaml           2.1.14     2016-11-12
##  source                                   
##  CRAN (R 3.2.5)                           
##  CRAN (R 3.3.2)                           
##  local                                    
##  CRAN (R 3.2.0)                           
##  Github (ChristofSeiler/BatchJobs@0e48846)
##  CRAN (R 3.3.2)                           
##  CRAN (R 3.2.0)                           
##  CRAN (R 3.3.2)                           
##  CRAN (R 3.3.2)                           
##  CRAN (R 3.2.0)                           
##  CRAN (R 3.3.0)                           
##  cran (@0.19-1)                           
##  CRAN (R 3.3.2)                           
##  CRAN (R 3.3.2)                           
##  local                                    
##  CRAN (R 3.3.0)                           
##  CRAN (R 3.3.2)                           
##  local                                    
##  CRAN (R 3.3.2)                           
##  CRAN (R 3.3.2)                           
##  cran (@0.6.12)                           
##  CRAN (R 3.2.3)                           
##  CRAN (R 3.3.2)                           
##  local                                    
##  local                                    
##  local                                    
##  CRAN (R 3.2.4)                           
##  CRAN (R 3.2.3)                           
##  cran (@0.3.6)                            
##  CRAN (R 3.2.0)                           
##  cran (@1.16)                             
##  CRAN (R 3.2.0)                           
##  CRAN (R 3.3.2)                           
##  Github (hadley/lazyeval@c155c3d)         
##  CRAN (R 3.3.2)                           
##  CRAN (R 3.2.0)                           
##  CRAN (R 3.3.2)                           
##  CRAN (R 3.3.2)                           
##  CRAN (R 3.3.2)                           
##  local                                    
##  CRAN (R 3.2.3)                           
##  CRAN (R 3.3.2)                           
##  local                                    
##  CRAN (R 3.3.2)                           
##  CRAN (R 3.2.5)                           
##  CRAN (R 3.2.0)                           
##  CRAN (R 3.3.2)                           
##  CRAN (R 3.3.0)                           
##  Github (rmcelreath/rethinking@a309712)   
##  CRAN (R 3.3.2)                           
##  CRAN (R 3.3.2)                           
##  CRAN (R 3.3.2)                           
##  CRAN (R 3.3.2)                           
##  CRAN (R 3.3.2)                           
##  CRAN (R 3.3.2)                           
##  CRAN (R 3.2.0)                           
##  CRAN (R 3.3.2)                           
##  local                                    
##  local                                    
##  cran (@1.1.5)                            
##  cran (@1.2.0)                            
##  CRAN (R 3.3.2)                           
##  local                                    
##  local                                    
##  CRAN (R 3.2.5)                           
##  cran (@2.1.14)

---

1. Hirshkowitz et al. (2015), National Sleep Foundation’s sleep recommendations: Methodology and results summary, Sleep Health.↩