Below, some examples of how to calculate and plot psychometric functions. I use ggplot2 to plot the graphs. The data (datLopezMoliner.txt) come from the paper López-Moliner & Linares 2006.

Plotting logistic functions using glm (function that implements generalized linear models in R).

library("ggplot2")

dat<-read.table("datLopezMoliner.txt",header=TRUE)


fit<-function(df) {

model<-glm(response~phase,data=df,family=binomial(link=logit))

grid <- data.frame(phase = seq(min(dat$phase),max(dat$phase), length = 200))

grid$response <- predict(model, grid, type = "response")

grid        

}

fits<-ddply(dat,.(subject,cond,speed),fit)


p<-ggplot(dat, aes(x=phase, y=response, colour=cond)) + facet_grid(speed~subject)

p<-p+stat_summary(fun.data = "mean_cl_boot",position=position_jitter(width=0.04))

p<-p+geom_line(data=fits)

p

Plotting cumulative gaussians and calculating the 95% bootstrap confidence intervals for the mean and the standard deviation of the cumulative gaussians (they are in the data frame psy.fits):

library("ggplot2") ; library("plyr"); library("Hmisc")

data<-read.table("datLopezMoliner.txt",header=TRUE)

data$speed<-as.factor(data$speed)

x.limits<-c(min(data$phase),max(data$phase))

split.variables<-c("subject","cond","speed")


psy.fit<-function(df) {

guess<-c(0.5,0.5)

number.simulations<-200

gaussian.op<-function(p,x,y){ ypred<- pnorm(x,p[1],p[2]); sum((y-ypred)^2) }

results<-ddply(df,.(phase),numcolwise(mean))

model<-optim(guess,gaussian.op,x=results$phase,y=results$response)$par

 

mean.sim<-c(); sd.sim<-c()

for (i in 1:number.simulations){

  sampling<-function(df) data.frame(response=mean(sample(df$response,replace=T)))

  results.sim<-ddply(df,.(phase),sampling)

  model.sim<-optim(guess,gaussian.op,x=results.sim$phase,y=results.sim$response)$par

  mean.sim[i]<-model.sim[1]; sd.sim[i]<-model.sim[2]

}

mean.ci<-quantile(mean.sim,c(0.025,0.975)); sd.ci<-quantile(sd.sim,c(0.025,0.975))

data.frame(mean=model[1],mean.inf=mean.ci[1],mean.sup=mean.ci[2],

                  sd=model[2], sd.inf=sd.ci[1],sd.sup=sd.ci[2])          

}

psy.fits<-ddply(data,split.variables,psy.fit)


psychometric<-function(df) {

grid <- data.frame(phase = seq(x.limits[1],x.limits[2], length = 200))

grid$response <- pnorm(grid$phase,mean=df$mean,sd=df$sd)

grid

}

psychometrics<-ddply(psy.fits,split.variables,psychometric)


p<-ggplot(data, aes(x=phase, y=response, colour=cond)) + facet_grid(speed~subject)

p<-p+stat_summary(fun.data = "mean_cl_boot",position=position_jitter(width=0.04))

p<-p+geom_line(data=psychometrics)

p

Plotting logistic functions using the library PsychoFun:

library("ggplot2")

library('PsychoFun')


dat<-read.table("datLopezMoliner.txt",header=TRUE)

x.limits<-c(min(dat$phase),max(dat$phase))

split.variables<-c("subject","cond","speed")


fit<-function(df) {

         dat.fit<-ddply(df,.(phase),summarize,resp=mean(response))

n<-length(df$response)/length(unique(df$phase))

dat.fit<-cbind(dat.fit,n)

set.up<-PsychoSetUp(dat.fit, plPrior="flat", t1Prior="flat",

t2Prior="flat", type = "logistic",chance = 0)

phase = seq(x.limits[1],x.limits[2],, length = 500)

MLEparameters <- MLEstimation(set.up)

r <- PsiFun(phase,MLEparameters,set.up)

grid<-data.frame(response=r,phase)

grid         

}

fits<-ddply(dat,split.variables,fit)


p<-ggplot(data=dat,aes(x=phase,y=response,color=cond))

p<-p+facet_grid(speed~subject)

p<-p+stat_summary(fun.data = "mean_cl_boot",position=position_jitter(width=0.04))

p<-p+geom_line(data=fits)

p

This library allows upper asymptotes smaller than 1 and lower asymptotes different from zero (the chance parameter should be changed).