#code by Alex O. Holcombe of University of Sydney, posted 5 October 2010 library (ggplot2) ##begin separate calculation for 2AFC #TwoAFCpCorrect <- function(dprime) { #proportion correct if d' = dprime # #probability that target response is greater than distracter strength # #assume target normally distributed with mean d' and variance 1 # #assume distracter normally distributed with mean 0 and variance 1 # #diff between two is distributed as normal variable with mean d' and variance 2, sd=sqrt(2) # #z-score point in distribution above which target greater than distracter is # # z = dprime/sqrt(2) want area under curve on thicker side of that line # p= pnorm(dprime/sqrt(2)) # p #} #ds=seq(.1,4,by=.1) #pCorrect= TwoAFCpCorrect(ds) #TwoAFC=data.frame(ds,pCorrect) #h<-ggplot(data=TwoAFC,aes(x=ds,y=pCorrect)) +geom_point() #xlabel('d-prime') #dprimeBetweenTargetAndDistracter = 2 ##simplification assuming target has no fluctuations #probDistracterExceedsTarget= 1- pnorm(dprimeBetweenTargetAndDistracter) #probNoDistracterExceeds = (1-probDistracterExceedsTarget)^m #proportion correct if d' = dprime #probability that target response is greater than distracter strength #assume target normally distributed with mean d' and variance 1 #assume distracter normally distributed with mean 0 and variance 1 #diff between two is distributed as normal variable with mean d' and variance 2, sd=sqrt(2) #z-score point in distribution above which target greater than distracter is # z = dprime/sqrt(2) want area under curve on thicker side of that line #p= pnorm(df\$dprime/sqrt(2)) #chance that will not exceed #Talking about not distribution of difference between one and other #chance that any one of them will exceed #Equal to 1- chance that none of them will exceed #here I'll assume independence, even though target distribution same for comparison with each distracter #p= p^(df\$m-1) #chance that in all target-distracter comparisons, target will come out on top #but this means your performance can go below chance. Because at least one of distracters #is likely to exceed target. But surely, individual distracter never has better chance of being #higher than target. What I calculated was the aggregate prob that #The chance that target is highest cannot go below 1/m #Let's say that chance individual distracter beats target is .5. Three, both beat=.25 #Chance that target beats all seems like .75*.75=.5625 . Only if target beats all do you #pick target? #Particularly obvious this is wrong if you consider the case wherein dprime=0 and the #prob that any distracter will exceed the target =.5. There, chance that target beats all #should be 1/m. But if I use .5*.5, chance that target beats distracter 1 is .5, and chance #that target beats distracter 2 is .5. Chance that beat both is .25, like getting 2 heads #With 2 coin flips, there are four possible outcomes. But with m alternatives, m possibilities #Comparing each distracter to a target. Distracter has 50% chance of beating target. #Target has 50% chance of beating distracter. .5*.5 =.25 chance of beating both distracters #Count the possibilities. Assign probability to each #Win Lose Lose, Lose Win Lose, Lose Lose Win #.5 .5 .5, #Target has 50% chance of beating distracter. .5*.5 =.25 chance of beating both distracters #should be .33333 # p(A-B) >0 and p(A-C)>0 makeIntegrandForDprime <- function(dprime,m) { integrand <- function(x) { dnorm(x-dprime) * pnorm(x)^(m-1) } return (integrand) } mAFCdprimeToPcorrect <- function(df) { #proportion correct if d' = dprime #probability that target response is greater than distracter strength #expects df\$dprime, df\$m integrand= makeIntegrandForDprime(df\$dprime,df\$m) pCorrect= integrate(integrand,lower= -Inf, upper= Inf) pCorrect = pCorrect\$value #because it's a structure with various values #print(pCorrect) df\$pCorrect = pCorrect return(df) #return(data.frame(pCorrect)) } mAFCdprimeToPcorrect <- function(df) { #proportion correct if d' = dprime #probability that target response is greater than distracter strength #expecting df\$dprime, df\$m, df\$lapseRate integrand= makeIntegrandForDprime(df\$dprime,df\$m) pCorrect= integrate(integrand,lower= -Inf, upper= Inf) pCorrect = pCorrect\$value #because it's a structure with various values #generate percent correct after incorporating lapse rate pCorrect = df\$lapseRate*(1/df\$m) + (1-df\$lapseRate)*pCorrect df\$pCorrect = pCorrect return(df) #return(data.frame(pCorrect)) } #make a plot showing curves that relate dprime to pCorrect as function of number of alternatives m=seq(2,10,by=1) lapseRate=.01 dfConds = expand.grid(m=m,dprime=seq(.5,4,by=.5),lapseRate=lapseRate) dfConds\$chance=1.0/dfConds\$m dfConds=ddply( dfConds, .(m,dprime), mAFCdprimeToPcorrect ) g<-ggplot( data=dfConds, aes(x=m, y= pCorrect,color=factor(dprime)) ) +ylim(0,1) +geom_point() +geom_line() g=g+ geom_line(aes(x=m,y=chance),color='black')+geom_point(aes(x=m,y=chance),color='black') #g=g+ geom_line(aes(x=m,y=(chance+ 1)/2),color='grey') #half-chance to ceiling was Eli's original suggestion g=g+ xlab('number of alternatives') +ylab('prob correct') +theme_bw() g=g+geom_text(aes(6,.2,label='chance'),color='black',size=3,angle=-20) g=g+ geom_hline(yintercept=.75,lty=2) #g=g+geom_text(aes(8,.6,label='halfChanceToCeiling'),color='grey',size=3,angle=-15) g ############################################################################################################ #to convert from percent correct to d-prime #need to do a search which minimizes discrepancy between predicted percent correct from mAFCdprimeToPcorrect and actual mAFCdprimeToPcorrectForOptim<- function(dprime,targetPctCorr,m,lapseRate) { #generate true percent correct for this dprime, and return discrepancy with target, squared temp=data.frame(m=m,dprime=dprime,lapseRate=lapseRate) pctCorr = mAFCdprimeToPcorrect(temp) # pctCorr= pctCorr\$pCorrect (targetPctCorr-pctCorr)^2 #squared error term } #test conversion pctCorr=.75; m=2; dprime=1 #fit = optim(dprimeGuess,mAFCpCorrectForOptim,m=m,pctCorr=pctCorr) #optim doesn't work for one-dimensional case fit = optimize(function(dprime) mAFCdprimeToPcorrectForOptim(dprime,pctCorr,m,lapseRate), interval=c(.0001,20) ) dprime = fit\$minimum #should be around .96 mAFCpCorrectToDprime<- function(pCorr,m,lapseRate) { #pCorrect = df\$lapseRate*(1/df\$m) + (1-df\$lapseRate)*pCorrect other direction #go from inclusion of lapses to exclusion pCorr = ( pCorr - lapseRate*(1.0/m) ) / (1-lapseRate) fit = optimize(function(dprime) mAFCdprimeToPcorrectForOptim(dprime,pCorr,m,lapseRate), interval=c(.0001,20) ) dprime= fit\$minimum dprime } mAFCpCorrectToDprime(pCorr=.95,m=2,lapseRate=.01) #test mAFCpCorrectToDprimeDf<- function(df) { dprime= mAFCpCorrectToDprime(df\$pCorr,df\$m,df\$lapseRate) } #for dataframes ############################################## #plot that relates a constant %corr (like Palmer used) to d-prime as a function of number of alternatives dfConds= expand.grid( m=seq(2,10,by=1), pCorr=seq(.5,.95,by=.05), lapseRate=lapseRate ) df=ddply( dfConds, .(m,pCorr,lapseRate), mAFCpCorrectToDprimeDf) names(df)[4]='dprime' quartz(); g<-ggplot( data=df, aes(x=m, y= dprime, color=factor(pCorr)) ) +geom_point() +geom_line() g=g+ xlab('number of alternatives') +theme_bw() g #now take Palmer's paper and try to determine dprime slope difference #he assumes the effect of having more items to process is to taking fewer samples of each. #so if you double set size from 1 to 2 items, then signal-to-noise ratio gets worse by factor of sqrt(2) #is 'signal-to-noise ratio' same as number of standard devations of separation? #Eli related proportion correct found at a particular m back to underlying d-prime. #So we saw how bad the effect of set size was on d-prime. #Now I'm trying to see how big an effect was found by Shaw and by Palmer #If Palmer reported percent correct, then I could calculate it directly. #However, Palmer only reports threshold line length, which then requires an assumption to relate back to dprime #He assumes it has the effect of shifting distribution to the right by some amount that's a constant w times the stim measure #e.g., line-length times the w. Since we don't know the constant, can't calculate absolute d-prime #instead, can only calculate change down to the multiplicative constant #perhaps generate plots for effect of set-size on some variable like TTC or line-length, assuming different constants #or generate plots relating line length threshold to d-prime, assuming different constants #so underlying dprime or d = t*w , where t is threshold and w is multiplicative factor #When w is large, then dprime corresponding to certain threshold is also large. Because doesn't take much stim diff #to improve difference between signal and distracters. #if I want to plot predicted t as function of m, rearrange t = d/w #parametrize for different dprimes, ws #plot that relates m to threshold stimulus difference, for different dprimes and ws threshPcorr =.75 dfConds= expand.grid( m=seq(2,10,by=1), dprime=seq(1,1,by=.5), w=seq(.5,2,by=.5), lapseRate=lapseRate ) mAFCpCorrectToThresh <- function(df) { #df must contain one m, dprime, and w #first calculate underlying dprime for this pct correct and m, then just multiply by constant w to get thresh dprime = mAFCpCorrectToDprime(threshPcorr,df\$m,df\$lapseRate) t = dprime / df\$w t } df=ddply( dfConds, .(m,dprime,w,lapseRate), mAFCpCorrectToThresh ) names(df)[5]='threshold' g<-ggplot( data=df, aes(x=m, y= threshold, color=factor(w)) ) +geom_point() +geom_line() g=g+ xlab('number of alternatives') +theme_bw() g h=g quartz() h h=h+scale_y_log10() #since curves are related by multiplicative factor, become parallel on semilog axis h h+scale_x_log10() #should just be able to understand how much d-prime may have elevated from the text of Palmer #I have to look at how much our dprime increased as function of set-size #Could simulate a threshold plot like Palmer from our data.. #Dprime should not increase more steeply than ... #Alternative approach is to plot prediction if can only process one item at a time, and compare to our data #Number of samples decreases by 1/m. Standard deviation is square root of number of samples, so dprime goes down #by sqrt(1/m), in other words decreases by half from 1 to 4 #So this applies to both distracters and targets. Their variance both inflates directly with m. sigma*m #Is that right? When processing m, each gets 1/m as many samples. Standard deviation = o/sqrt(1/m) #Dprime then equals original * 1/(sqrt(1/m)) #so, plot dprime elevation as function of m and compare to our results? #need some baseline number, and then can predict the elevation if m increases by given amount dprimeInflateCapacity<- function(baselineDprime, mMultiplier) { #to predict dprime for m, take ratio of increase of set size and old dprime dPrime = baselineDprime / sqrt(mMultipier) } #Eventually we shall specify processing capacity and show different curves for different capacity, for example #m on abscissa, predicted d-prime #Baseline data is proportion correct found for m=2 and corresponding calculated dprime #Then as function of m, predicted change in proportion correct according to infinite capacity, vs. according to #capacity=1. mAFCbaselinePcorrToHigherSetSize <- function(df) { baselinePcorr=df\$baselinePcorr; lapseRate=df\$lapseRate; baselineM=df\$baselineM; m=df\$m; capacity=df\$capacity; #first, convert baseline pcorr to dprime by way of lapse rate #adjust pCorr to underlying pCorr if hadn't had any lapses baselinePcorr = ( baselinePcorr - lapseRate*(1.0/baselineM) ) / (1-lapseRate) #calculate dprime baselineDprime= mAFCpCorrectToDprime(baselinePcorr,baselineM,lapseRate) newOldMratio = m/baselineM #calculate new d-prime that represents distance between target and a single distracter, given capacity #if capacity =1, samples previously devoted to reducing original sigma are now spread equally among all #so samples available to one are divided by ratio of new to old if (capacity == 1) newDprime = baselineDprime / sqrt(newOldMratio) else newDprime = baselineDprime #convert back to proportion correct dg=data.frame(dprime=newDprime,m=m,lapseRate=lapseRate) pCorr= mAFCdprimeToPcorrect(dg)\$pCorrect df\$pCorr = pCorr df } baselinePcorr=.90; lapseRate=.01 d=mAFCpCorrectToDprime(baselinePcorr,m=2,lapseRate=lapseRate) dfConds=expand.grid(m=seq(2,12,by=1),baselinePcorr=baselinePcorr, baselineM=2, lapseRate=lapseRate, capacity=c(1,100) ) df=ddply( dfConds, .(m,baselinePcorr,baselineM,lapseRate,capacity), mAFCbaselinePcorrToHigherSetSize ) quartz() g<-ggplot( data=df, aes(x=m, y= pCorr, color=factor(capacity)) ) +geom_point() +geom_line() +ylim(0,1) g=g+ xlab('number of alternatives') +theme_bw() g=g+ geom_line(aes(x=m,y=1/m),color='black')+geom_point(aes(x=m,y=1/m),color='black') g=g+geom_text(aes(6,.2,label='chance'),color='black',size=3,angle=-20) g #would be good to have capacity as a parameter. So, when processing 2 it gets as many samples as does 1. But when #processing 3, all three get fewer samples. #Other types of resources models would change the nature of the function relating #in theory could estimate w factor or function relating TTC threshold to sensitivity from single-shot experiment #Perhaps Palmer's method of plotting threshold can be thought of as equivalent to plotting dprime plus multiplicative factor? #In that case, I can recover dprime elevation?