Simulating fMRI Data

The code here attempts to first simulate the ‘neuronal’ response with a particular graphical structure and then the response at each node is convolved with a hemodynamic resopnse. Finally, I test the results using simple correlations and the Patel’s tau (a causal modeling approach based on ascendency). # Generate Data

Simulate Neuronal Response

We’ll first simulate the ‘neuronal’ response as random data S and then simulate the connectivity by multiplying matrix W by S.

I end up not using this but thought I’d leave it here.

library(MASS)
# Specify the covariance between the nodes
W <- matrix(c(
  8,2,4,
  1,4,7,
  2,4,6
), 3, 3, byrow=T)
S <- matrix(rnorm(1000*3), 1000, 3)
A <- t(W %*% t(S))

Since I don’t think what I did above is right, I’m trying a function from the the pcalg (Methods for Graphical Models and Causal Inference) package. This will simulate the data based on a directed acyclic graph (DAG).

library(pcalg)
library(igraph)
## generate random DAG
# p = number of nodes
# prob = probability of connecting a node to another node
# lB and uB are lower/upper limits of edge weights
p <- 20
rDAG <- randomDAG(p, prob = 0.2, lB=0.1, uB=1)
## generate 1000 samples of DAG using standard normal error distribution
n <- 2000
d.normMat <- rmvDAG(n, rDAG, errDist="normal")

We can plot the DAG followed by a sample of the simulated time-series. I only plot the first 200 time-points for the first 4 nodes.

## plot (not sure how to plot this to show that it is top-down)
g    <- graph.adjacency(t(wgtMatrix(rDAG)), mode="directed", weighted=T, diag=F)
plot.igraph(g)

## alternative plot
# library(Rgraphviz)
# plot(rDAG, main = "randomDAG(20, prob = 0.2, ..)")
## plot some of the data
plot.ts(d.normMat[1:200,1:4], main="Sample Time-Series from 4 Nodes")

Convolve with Hemodynamic Response

We use the canonical hemodynamic respnse (mainly because it’s faster), which is a double-gamma function. Note that commented out is the ballon function and that I’m scaling the peak response to 1.

library(neuRosim)

## Loading required package: deSolve
## This is neuRosim 0.2-12
## neuRosim is BETA software! Please report any bugs.

accuracy <- 0.1 # convolution to be done at this time resolution
totaltime <- n*accuracy
#s.conv <- balloon(d.normMat[,1], totaltime, accuracy)
#s.conv <- s.conv/max(s.conv)
## get the hrf
tpts   <- seq(accuracy, totaltime, accuracy)
chrf   <- canonicalHRF(tpts, verbose = FALSE)
conv.normMat <- sapply(1:ncol(d.normMat), function(i) {
  #s.conv <- balloon(d.normMat[,i], totaltime, accuracy)
  s.conv <- convolve(chrf, rev(d.normMat[,i]))
  s.conv <- s.conv/max(s.conv)
  s.conv
})

We should also add some noise to each node’s time-series using functions from the neuRosim package. I have chosen to add white noise temporally correlated, and low noise with an SNR of 10. See simTSrestingstate for more info since that’s the function I based the below code on.

SNR   <- 10
nscan <- nrow(conv.normMat)
simMat <- matrix(0, nrow(conv.normMat), ncol(conv.normMat))
for (i in 1:ncol(conv.normMat)) {
  sigma <- abs(mean(conv.normMat[,i]))/SNR
  # white noise
  n.white <- c(systemnoise(dim = c(1), sigma = sigma, nscan = nscan, 
                           type = "rician", verbose = F, vee = 1))
  # auto-regressive noise
  n.temp <- c(temporalnoise(dim = c(1), sigma = sigma, nscan = nscan, 
                            rho = 0.2, verbose = F))
  # low-frequency noise
  n.low <- c(lowfreqdrift(dim = c(1), freq = 128, nscan = nscan, 
                          TR = accuracy, verbose = F))
  # physiological noise (heart rate and respiratory rate)
  n.phys <- c(physnoise(dim = c(1), sigma = sigma, 
                        nscan = nscan, TR = accuracy, freq.heart = 1.17, 
                        freq.resp = 0.2, verbose = F))
  # combine all the noise together with some arbitrary amount of weights
  w <- c(0.2, 0.2, 0.1, 0.3)
  n.noise <- (w[1]*n.white + w[2]*n.temp + w[3]*n.low + w[4]*n.phys)/sqrt(sum(w^2))
  # add in the noise
  simMat[,i] <- conv.normMat[,i] + (n.noise - mean(n.noise))
}

We can plot the HRF used in the convolution with the simulated neuronal response followed by plotting the actual convolved neural response for the first 4 nodes (as before).

plot.ts(conv.normMat[1:200,1:4], main="Convolved Time-Series from 4 Nodes")

plot.ts(simMat[1:200,1:4], main="Convolved Time-Series with Noise from 4 Nodes")

We will now downsample the data to get at what we might really collect with fMRI. We take every 10th time-point to get 1s TR.

simMat2 <- simMat[seq(1,length(tpts),by=10),]

Testing the Network

We can try to recover the causal network using a variety of measures. The measures that are good would be those that can succesfully recover the network.

For now, I will try the Patel’s tau (cite), which we had previously discussed. From my understanding, this measure assumes A causes B (A -> B) if A is ascendent to B. For instance, in the case of cloud -> rain, we know that when there are clouds there is sometimes rain and sometimes not but when there is rain there are always clouds.

I test the patel’s tau with simulated neuronal data and then with the convoled data. It doesn’t appear to do very well with either data (near chance at recovering the directionality).

patel.tau <- function(x, S=NULL, to.scale=F, verbose=F) {
  if (to.scale) x <- scale(x)
  if (is.null(S)) S <- matrix(1, ncol(x), ncol(x))
  
  # 1. Map the data into the interval [0,1]:
  if (verbose) cat("map data into the interval [0,1]\n")
  X10 <- apply(x, 2, quantile, 0.1) # cutoff for the 10th percentile
  X90 <- apply(x, 2, quantile, 0.9) # cutoff for the 90th percentile
  a   <- sweep(sweep(x, 2, X10), 2, X90-X10, FUN="/") # (X-X10)/(X90-X10)
  x2  <- apply(a, c(1,2), function(v) max(min(v,1),0)) # keep in [0,1] interval
  
  # 2. theta1 = X2*Y2 or x2*t(x2)
  if (verbose) cat("theta1\n")
  theta1 <- crossprod(x2)
  
  # 3. theta2 = X2 * (1-Y2)
  if (verbose) cat("theta2\n")
  theta2 <- crossprod(x2, 1-x2)
  
  # 4. theta3 = (1-X2)*Y2
  if (verbose) cat("theta3\n")
  theta3 <- crossprod(1-x2,x2)
  
  # 5. Set tau
  if (verbose) cat("tau\n")
  tau  <- matrix(0, ncol(x), ncol(x))
  inds <- theta2 > theta3
  tau[inds] <- 1 - (theta1[inds] + theta3[inds])/(theta1[inds] + theta2[inds])
  tau[!inds] <- (theta1[!inds] + theta2[!inds])/(theta1[!inds] + theta3[!inds]) - 1
  
  # 6. Directed unweighted graph
  if (verbose) cat("dag\n")
  dag <- S * (tau<0)
  
  list(dag=dag, tau=tau)
}

nconns <- sum(wgtMatrix(rDAG)!=0)
symMat <- (wgtMatrix(rDAG)!=0) + (wgtMatrix(rDAG,F)!=0)

cat("Testing patel's tau with non-noise neuronal data\n")
res1 <- patel.tau(d.normMat, symMat)
table(patel=res1$dag, orig=(wgtMatrix(rDAG)!=0)*1)
#res1 <- patel.tau(d.normMat[seq(1,length(tpts),by=10),], symMat)
#table(patel=res1$dag, orig=(wgtMatrix(rDAG)!=0)*1)

cat("Convolved data\n")
res1 <- patel.tau(simMat2, symMat)
table(patel=res1$dag, orig=(wgtMatrix(rDAG)!=0)*1)

## Testing patel's tau with non-noise neuronal data

##      orig
## patel   0   1
##     0 338  22
##     1  22  18

## Convolved data

##      orig
## patel   0   1
##     0 344  16
##     1  16  24

Scratch

The below is for future reference when using the neuRosim package and I want to generate some simulated task data

totaltime <- 200
onsets <- seq(1, 200, 40)
dur <- 20
s <- stimfunction(totaltime = totaltime, onsets = onsets,
                  durations = dur, accuracy = 0.1)
canonical <- specifydesign(totaltime = 200, onsets = list(onsets),
                           durations = list(dur), effectsize = 1, TR = 2,
                           conv = "double-gamma")
balloon <- specifydesign(totaltime = 200, onsets = list(onsets),
                         durations = list(dur), effectsize = 1, TR = 2, 
                         conv = "Balloon")

## Warning: Default parameter values are used. See ?balloon.fnc for more
## details.

plot.ts(s, main="BoxCar Simulated Duration of Task Events")

plot.ts(canonical, main="Convolved with Double-Gamma")

plot.ts(balloon, main="Convolved with Balloon")