## Reads in the fMRI data in brainm.dat (a 128 by 234 matrix) ## and arranges them into a list with 54 components ## data$L1S1, data$L1S2, ... data$L1S6, ... data$L9S1, data$L9S2, ... data$L9S6. ## L is for 'location', S is for 'stimulus' ## Each component is a matrix with one column (a time series of length 128) ## for each subject. Numbers of subjects varies: ## ## S1: "Awake-Brush" 5 subjects ## S2: "Awake-Heat" 4 subjects ## S3: "Awake-Shock" 5 subjects ## S4: "Low-Brush" 3 subjects ## S5: "Low-Heat" 5 subjects ## S6: "Low-Shock" 4 subjects ## Then Figure 7.1 in the text (mean responses per stimulus at location 1) is duplicated ## Similar figures are produced for the other 8 locations ## Read in the data - a 128 by 234 matrix z.all=ts(read.table("http://www.stat.ualberta.ca/~wiens/stat479/S&Sdatasets/brainm.dat")) # or: z.all = read.table("brainm.dat") ############## ## Arrange data in 9 groups of 26 columns; each group is a "location": L = NULL for (j in 1:9) { L = cbind(L, z.all[, j + 9*(0:25)]) } colnames(L) = NULL # L has all data at location 1 (first 26 columns), # then all data at location 2 (next 26 columns), ... # then all data at location 9 (final 26 columns) ############## ## Set the indicators of "location" where.loc = vector(length = 234) for (l in 1:9) where.loc[(1:26) + (l-1)*26] = l #where.loc # L[,where.loc == l] consists of the data from location "l" Location.names =c("Cortex 1", "Cortex 2", "Cortex 3", "Cortex 4", "Caudate", "Thalamus 1", "Thalamus 2", "Cerebellum 1", "Cerebellum 2") cat("Here are the labels of the locations:","\n") as.matrix(Location.names) ############## ## Set the indicators of "stimuli" stim = c(rep(1,5), rep(2,4), rep(3,5), rep(4,3), rep(5,5), rep(6,4)) where.stim = rep(stim,9) #where.stim # L[,where.stim == s] consists of the data from stimulus "s" Stimulus.names = c("Awake-Brush", "Awake-Heat", "Awake-Shock", "Low-Brush", "Low-Heat", "Low-Shock") cat("Here are the labels of the stimuli:","\n") as.matrix(Stimulus.names) ############## ## Create a list with components L1S1, L1S2, ... L1S6, L2S1, ..., L9S6 ## Each component is a matrix with one column (a series of length 128) per subject data = vector("list", 54) for (l in 1:9) {for(s in 1:6) { data[[(l-1)*6+s]] = L[, (where.loc == l & where.stim == s)] }} datanames = vector(length=54) for(l in 1:9) datanames[(l-1)*6+1:6]=paste("L",l,"S",1:6,sep="") names(data) = datanames #data # A check: #L5S3test = L[, (where.loc == 5 & where.stim == 3)] #L5S3test - data$L5S3 # Should be all zeros ############## ## Figure 7.1 in text; average reponses to stimuli at location 1: par(mfcol = c(3,2)) for(s in 1:6) plot(ts(rowMeans(data[[s]])), ylim = c(-1,1), ylab = Stimulus.names[s]) title(main = paste("\nMean responses at", Location.names[1], "to various stimuli"), outer = T) ############## ## Similar plots at other locations: for(l in 2:9) { win.graph() par(mfcol = c(3,2)) for(s in 1:6) plot(ts(rowMeans(data[[(l-1)*6+s]])), ylim = c(-1,1), ylab = Stimulus.names[s]) title(main = paste("\nMean responses at", Location.names[l], "to various stimuli"), outer = T) }