#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?