Function to fit unequal variance SDT model to confidence ratings data

ratingsROC_MLE.R

ratingsROC <- function(Fcount =c(2,8,37,23,30), Hcount = c(61,15,15,5,4)) { 
  #R script to fit ROC data from confidence rating task
  #CMR + WJH Mar 26 2015

  #data from Macmillan & Creelman 2005 Example 3b
  #Low Frequency words, confidence rating recognition task
  # Counts are ordered from 'sure old' to 'sure new' response for each stim class
  
  fa_probs <- cumsum(Fcount)/sum(Fcount)
  h_probs <- cumsum(Hcount)/sum(Hcount)
  
  gSDT<-function(p,obs,useMLE =TRUE) {  
    #this function fits data from a 5-level ratings condition with the SDT model
    preds<-matrix(0,ncol=2,nrow=nrow(obs),dimnames = list(NULL,names(obs))) #this matrix will hold predicted counts
    n <-sum(obs$FA) #total number of lure trials = total number of target trials in this example
  
    # Make sure c1 > c2 > c3 > c4 or return a redonkulus likelihood
    crit <- p[grep('c',names(p),ignore.case=TRUE)]
    if ( !all(diff(crit)<0) ) {
      return(1000000)
    }
    d <- p['d']
    s <- p['s']
  
    preds[1:nrow(preds)-1,'FA'] <- pnorm(crit,lower.tail=F)
    preds[2:(nrow(preds)-1),'FA'] <-diff(preds[1:(nrow(preds)-1),'FA'])
    preds[,'FA'] <- n*preds[,'FA']
    preds[nrow(preds),'FA'] <- n-sum(preds[,'FA'])
  
    preds[1:nrow(preds)-1,'H'] <- pnorm(crit, mean = d, sd = 1/s,lower.tail=F)
    preds[2:(nrow(preds)-1),'H'] <-diff(preds[1:(nrow(preds)-1),'H'])
    preds[,'H'] <- n*preds[,'H']
    preds[nrow(preds),'H'] <- n-sum(preds[,'H'])
  
    if (useMLE) {
      ll <- as.matrix(obs)*log(preds/n)
      ll[is.infinite(ll) | is.nan(ll)] <- 0
      return(-sum(ll))
    } else {
      #compute G^2
      g2 <- 2*as.matrix(obs)*(log(as.matrix(obs))-log(preds))
      g2[is.infinite(g2) | is.nan(g2)] <- 0  
      return(sum(g2))
    }

  }
  nIntervals <- length(Hcount)-1
  p_names <- c("s", "d", paste('c',1:nIntervals,sep=''))
  start <- rep(0,length(p_names))
  names(start) <- p_names
  # max difference in s2/s1 ratio is 2, max dprime is 3
  # most conservative criteron maxes out a 3, other more liberal ones max out at 2
  max_vals <- c(2,4,seq(3,3-.01*(nIntervals-1), by=-.01) )
  # smallest difference in s2/s1 ratio is .1, min dprime is dprime is 0
  # most conservative criteron maxes out at -3, other more liberal ones max out at .01
  min_vals <- c(.1,0,seq(-3,-3 -.01*(nIntervals-1), by=-.01))
  err<-100000.0 #start with large value, which we then seek to minimize  
  
  for (m in 1:100) {  #do it 100 time with different starting parameters) 
    start[] <- c(runif(2),sort(runif(4),dec=TRUE))
    # fit rating data, starting with 6 initial parameter values sampled from a uniform distribution
    fit<-nlminb(start,gSDT,control=list(eval.max=1500,iter.max=1500),
               lower=min_vals,upper=max_vals, obs = data.frame(H=Hcount,FA = Fcount), useMLE=FALSE)
    if(fit$objective<err)  {
      err<-fit$objective
      bestfit <- fit
    }
  }
  
  # Known best fit
  # derp <- structure(c(0.560440499277371, 2.53406429320426, 2.03904738869317, 
  #             1.27366387425955, 0.0876087175909641, -0.53067977664763), 
  #             .Names = c("s", "d", "c1", "c2", "c3", "c4"))
  #   knownfit <-nlminb(derp,gSDT,control=list(eval.max=1500,iter.max=1500),
  #              lower=min_vals,upper=max_vals, obs = data.frame(H=Hcount,FA = Fcount))
  
  slope<-bestfit$par['s'] #parameters are listed in order:  slope, d'1, c1,c2, c3, c4
  d2<-slope*bestfit$par['d']  #equation 3.1
  da<-sqrt(2/(1+slope^2))*d2 #equation 3.4
    
  Fcurve<-seq(.001,.999,.01) #range of false alarms
  Hcurve<-pnorm(da/sqrt(2/(1+slope^2)) + slope*qnorm(Fcurve))  #what is the hit rate given parameters and F?
   
  plot(fa_probs, h_probs,type="p",pch=2,cex=1.5,xlab="False Alarm Rate",ylab="Hit Rate",xlim=c(0,1),ylim=c(0,1),main="M&C Exp 3b LF words",las=1,cex.lab=1.2)
  text(.8,.2,paste("da = ",round(da,3),sep=""),cex=1.2)
  text(.8,.1,paste("slope = ",round(slope,3),sep=""),cex=1.2)
  points(Fcurve,Hcurve,type="l",lwd=2)
  return(bestfit)
}