R/windowWDM.R

In akcochrane/ACmisc: Aaron Cochrane's miscellaneous functions

#' Find a vector of Wiener Diffusion Model parameters
#' 
#' Estimates a set of Wiener Diffusion Model parameters for a vector of formatted
#' data. Estimation uses \code{\link{RWiener::wdm}} and data must be formatted as a 
#' 1-dimensional vector according the to \code{RWiener} package (i.e., RT * 2 
#' * (boundary=='upper' - 0.5) ). 
#' Assumes that RTs are associated with equally-spaced and sorted timepoints (e.g., trials of a behavioral study).
#' First fits a \code{\link{RWiener::wdm}} for the entire RT vector, and applies some of these parameters
#' (NDT, bias, and possibly BS) to the entire RT vector. Next takes a set of windows, equally spaced and centered at
#' \code{windowSize/2}, and estimates the DR (and possibly BS) for each window. Averages overlapping windows. 
#' Returns a data frame of estimated parameters, including 95\% CI.
#'
#' @param dat Response time vector. Positive values indicates upper-boundary while negative values indicate lower-boundary.
#' @param windowSize Total width of each window for estimating parameters.
#' @param fit_BS Logical. If \code{FALSE} only drift rate is estimated, if \code{TRUE} both drift rate and boundary separation are estimated.
#'
#' @export
#'
#' @examples
#' \dontrun{
#' dat <- (rnorm(200,.2,.05) + rexp(200,1 / .2))*sample(c(-1,1),200,replace=T)
#'     m_drY_bsN <- windowWDM(dat,fit_BS = F)
#'     m_drY_bsY <- windowWDM(dat,fit_BS = T)
#'     ACmisc::pairsplot(cbind(trialNum = 1:200, RT = abs(dat),m_drY_bsY[,1:6]))
#'     }
windowWDM <- function(dat,windowSize=20,fit_BS=FALSE){
  
  fit_DR=T ## ideally this would be an argument, but only fitting BS always returns errors.
  
  if(F){  
    
    dat <- (rnorm(200,.2,.05) + rexp(200,1 / .2))*sample(c(-1,1),200,replace=T) 
  
    m_drY_bsN <- windowWDM(dat,fit_BS = F)
    m_drY_bsY <- windowWDM(dat,fit_BS = T)
    
    # m_drN_bsY <- windowWDM(dat,fit_DR = F, fit_BS = T)
    
    psych::pairs.panels(cbind(trialNum = 1:200, RT = abs(dat),m_drY_bsN[,1:3]))
    
    # psych::pairs.panels(cbind(1:200,m_drN_bsY[,4:6]))
    
    psych::pairs.panels(cbind(trialNum = 1:200, RT = abs(dat),m_drY_bsY[,1:6]))
  }
  
  require(RWiener)
  
  overallMod <- wdm(dat)
  overallPars <- overallMod$coefficients
  overallCI <- summary(overallMod)$cint
  
  if(fit_DR){drEst <- data.frame( Estimate = rep(NA,length(dat)), ci975 = NA , ci025 = NA)
  }else{drEst <- data.frame(Estimate = overallPars['delta'], ci025 = overallCI['delta',"2.5 %"], ci975 = overallCI['delta',"97.5 %"])}
  if(fit_BS){bsEst <- data.frame( Estimate = rep(NA,length(dat)), ci975 = NA , ci025 = NA)
  }else{bsEst <- data.frame(Estimate = overallPars['alpha'], ci025 = overallCI['alpha',"2.5 %"], ci975 = overallCI['alpha',"97.5 %"])}
  
  half_kernel <- round(windowSize/2)
  
  kernel_centers <- seq(half_kernel,length(dat) - half_kernel, half_kernel)
  
  for(curKernel in kernel_centers){
    
    curIndices <-  max(curKernel - half_kernel + 1, 1) : min(curKernel+half_kernel, length(dat))
    
    curDat <- dat[curIndices]
    
  try({ ## because optimization will inevitably fail for some windows, if the window size is too small
    
    switch(paste(fit_DR,fit_BS)
           
           ,'TRUE TRUE' = {
             curMod <- wdm(curDat,tau = overallPars['tau'],beta = overallPars['beta'])
             
             drEst[curIndices,'Estimate'] <- 
               rowMeans(cbind(drEst[curIndices,'Estimate'],curMod$coefficients['delta']),na.rm=T ) 
    drEst[curIndices,'ci025'] <- 
      rowMeans(cbind(drEst[curIndices,'ci025'],summary(curMod)$cint['delta',"2.5 %"]),na.rm=T )
    drEst[curIndices,'ci975'] <- 
      rowMeans(cbind(drEst[curIndices,'ci975'],summary(curMod)$cint['delta',"97.5 %"]),na.rm=T )
    
    bsEst[curIndices,'Estimate'] <- 
      rowMeans(cbind(bsEst[curIndices,'Estimate'],curMod$coefficients['alpha']),na.rm=T ) 
bsEst[curIndices,'ci025'] <- 
  rowMeans(cbind(bsEst[curIndices,'ci025'],summary(curMod)$cint['alpha',"2.5 %"]),na.rm=T )
bsEst[curIndices,'ci975'] <- 
  rowMeans(cbind(bsEst[curIndices,'ci975'],summary(curMod)$cint['alpha',"97.5 %"]),na.rm=T )
           }

,'TRUE FALSE' = {
  curMod <- wdm(curDat,alpha = overallPars['alpha'],tau = overallPars['tau'],beta = overallPars['beta'])
  
  drEst[curIndices,'Estimate'] <- 
    rowMeans(cbind(drEst[curIndices,'Estimate'],curMod$coefficients['delta']),na.rm=T ) 
drEst[curIndices,'ci025'] <- 
  rowMeans(cbind(drEst[curIndices,'ci025'],summary(curMod)$cint['delta',"2.5 %"]),na.rm=T )
drEst[curIndices,'ci975'] <- 
  rowMeans(cbind(drEst[curIndices,'ci975'],summary(curMod)$cint['delta',"97.5 %"]),na.rm=T )

}

,'FALSE TRUE' = {
  curMod <- wdm(curDat,delta = overallPars['delta'],tau = overallPars['tau'],beta = overallPars['beta'])
  
  bsEst[curIndices,'Estimate'] <- 
    rowMeans(cbind(bsEst[curIndices,'Estimate'],curMod$coefficients['alpha']),na.rm=T ) 
bsEst[curIndices,'ci025'] <- 
  rowMeans(cbind(bsEst[curIndices,'ci025'],summary(curMod)$cint['alpha',"2.5 %"]),na.rm=T )
bsEst[curIndices,'ci975'] <- 
  rowMeans(cbind(bsEst[curIndices,'ci975'],summary(curMod)$cint['alpha',"97.5 %"]),na.rm=T )

}
)
      
   },silent=T) # the kernel-fitting try()

} # close kernel loop
  
  colnames(drEst) <- paste0('dr_',colnames(drEst))
  colnames(bsEst) <- paste0('bs_',colnames(bsEst))
  
  suppressWarnings({
  outTable <- data.frame(
    drEst
    ,bsEst
    ,ndt_Estimate = overallPars['tau']
    ,ndt_c025 = overallCI['tau',"2.5 %"]
    ,ndt_c975 = overallCI['tau',"97.5 %"]
    ,bias_Estimate = overallPars['beta']
    ,bias_c025 = overallCI['beta',"2.5 %"]
    ,bias_c975 = overallCI['beta',"97.5 %"]
  )
  })
  
  if(any(apply(outTable,2,function(x){mean(is.na(x))>.25}))){
    warning('Estimation often failed. Using a larger window size is recommended.')
  }
  
  attr(outTable,'overall_model') <- overallMod
  attr(outTable,'overall_parameters') <- overallPars
  attr(outTable,'overall_CI') <- overallCI
  
  return(outTable)
}