Compare the Double-wavelet transform approach and Mean voxel approach for Tasked-Induces fMRI data

This is the R shinyapp I made for my dissertation work “Using double-wavelet transform to take spatial-temporal correlation in tasked-induced fMRI data analysis”. The hightlights of this work as follows:

Feel free to play it and let me know your thought! Have fun! Link: http://statcomp2.vanderbilt.edu:37212/block_simu/

shinyUI

# ui.R
source("block_functions.R")

shinyUI(navbarPage("Simulation for comparing Double Wavelet Transform and Mean Voxel Approach",
                   tabPanel("Block Design",
                            fluidPage(
                              fluidRow(
                                column(3,
                                       h4("Simulation Setting") ,
                                       numericInput("N.sim", 
                                                    label = h5("Number of Simulation"),
                                                    value = 100),
                                       numericInput("N.subj",
                                                    label = h5("Number of subjects per Simulation"),
                                                    value = 5),
                                       hr(),
                                       numericInput("N.dim1", 
                                                    label = h5("Roi Length"), 
                                                    value = 10),
                                       numericInput("N.dim2", 
                                                    label = h5("Roi Width"), 
                                                    value = 10),
                                       h6("Note: Number of voxels = Length * Width"),
                                       hr(),
                                       numericInput("effect_size", 
                                                    label = h5("Effect size for Activate ROI"), 
                                                    value = 0.2),
                                       selectInput("N.time",h5("Length of Time series"),
                                                   choices = c(128), selected = 128)
                                       ),
                                column(4,
                                       actionButton("go", "Start Simulation"), 
                                       hr(),
                                       tableOutput("table"),
                                       verbatimTextOutput("value")
                                       
                                ),
      
                                # correlation size 
                                column(2,
                                h4("Double Wavelet"),
                                # double wavelet option
                                selectInput("waveP",h5("Spatial wavelet"),
                                            choices = allwave, selected = "db3"),
                                selectInput("waveT",h5("Temporal wavelet"),
                                            choices = allwave, selected = "sym8")
                                
                                ),
      
                              column(3,
                                     h4("Temporal correlation"),
                                     
                                     sliderInput("phi", label = h5("Autoregressive One parameter"), 
                                                 min = 0, max = 1, value = 0.6),
                                     
                                     selectInput("spat_cor",h4("Spatial Correlation"),
                                                 choices = c( "Independent" ,"Exponential", "Identical"),
                                                 selected = "Exponential"),
                        
                                     numericInput("randomsigma", 
                                                  label = h5("Random Variance"), 
                                                  value = 3), 
                                     
                                     
                              numericInput("spat_phi", 
                                           label = h5("Exponential Decaying parameter"), 
                                           value = 2),
                              
                              numericInput("phi_sigma", 
                                           label = h5("Exponential Decaying Sigma"), 
                                           value = 1.5),
                              
                              numericInput("GauSigma", 
                                           label = h5("Gaussian Smoothing Sigma"), 
                                           value = 1.5)
                              
                              )
                              
      ) 
    )
  )
))

ShinyServer

# server.R

source("block_functions.R")

shinyServer(function(input, output) {
  
  analysis <- function(input=input,X=X,con=con){
    shiny::isolate({
      source("block_functions.R")
      
    if (input$spat_cor ==  "Independent" ){
      genfunc <- gen2roi_error_ind
    } else if (input$spat_cor ==  "Exponential" ){
      genfunc <- gen2roi_error_exp
    } else if (input$spat_cor ==  "Identical" ){
      genfunc <- gen2roi_error_same
    }
      
      data <- genfunc(input,X)
      data <- normalseries(data)
      data_smooth <- spatial_smoothing(data=data, input)
      result_mean <- block_mean(data_smooth, X=X, con=con)
      result_dw <- block_dw(data,input, X=X, con=con, level=1)
      result <- rbind( c(result_mean, "Mean") , c(result_dw, "DW") )
    }
    )
      return(result)
    
  }
  
  output$value <- renderPrint({ print(paste0("Simulation ", input$go))  })
  
  resulttable <- eventReactive(input$go, {

    # Parallel computing
    cl <- makeCluster(8, 'PSOCK')
    clusterExport(
    cl, varlist=c("input","analysis"),
    envir=environment())
    
    mat <- parSapply(cl, 1:(input$N.sim*input$N.subj), function(i) {
      source("block_functions.R")
      analysis(input, X, con)
    },simplify = FALSE)
    stopCluster(cl)
    

    data <- as.data.frame(do.call("rbind", mat))
    colnames(data) <- c("roi1", "roi2", "method")
    data$roi1 <- as.numeric(as.character(data$roi1)) 
    data$roi2 <- as.numeric(as.character(data$roi2)) 
    
    meandata <- data[ data$method == "Mean" ,]
    dwdata <- data[ data$method == "DW" ,]
    
    errorrate <- function(testdata, Nsimu = input$N.sim, Nsub = input$N.subj){
      resultP <- rep(NA,Nsimu)
      
      for (i in 1: Nsimu){
        test1 <- t.test( testdata[ (1:Nsub) + (i-1)*Nsub  ] )
        resultP[i] <- test1$p.value
      }
      return(resultP)
    }
    
    meantypeI <- mean(errorrate(meandata$roi2) < 0.05)
    meantypeII <- mean(errorrate(meandata$roi1) > 0.05)
    
    dwtypeI <- mean(errorrate(dwdata$roi2) < 0.05)
    dwtypeII <- mean(errorrate(dwdata$roi1) > 0.05)
    
    result <- rbind( c(meantypeI , meantypeII ), c(dwtypeI, dwtypeII) )
    colnames(result) <- c("Type I Error", "Type II Error")
    rownames(result) <- c("Mean-Voxel", "Double-Wavelet")
    
    result
      })
  
  output$table <- renderTable({
    resulttable()
    }, rownames = TRUE)

})