library("ggplot2") #plotting

library("plyr")  # to split and work with dataframes

library("Hmisc") # to use the function mean_cl_boot

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) { # funtion to fit cumulative gaussians

guess<-c(0.5,0.5) # guess mean and slope

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    

data.frame(mean=model[1],sd=model[2])           

}

psy.fits<-ddply(data,split.variables,psy.fit) # plyr (fitting by split variable)


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)


theme_set(theme_bw()) # plotting; the error bars are bootstrapping conf int

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)) #avoid overplotting

p<-p+geom_line(data=psychometrics) # plot the psychometric functions

p

Below, some examples of how to calculate and plot psychometric functions using “plyr” and “ggplot2” packages. I used some data (datLopezMoliner.txt) from the paper López-Moliner & Linares 2006.

1) Plotting cumulative gaussians:

2) The same graph but fitting logistic functions can be generated with the code below:

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))


psy.fit<-function(df) {

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

model<-glm(response~phase,data=results,family=binomial)

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

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

grid         

}

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


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=psy.fits)

p

3) The code below generates the first graph but also calculates 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

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))


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,.(subject,cond),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$response.inf <- pnorm(grid$phase,mean=df$mean.inf,sd=df$sd)

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

grid

}

psychometrics<-ddply(psy.fits,.(subject,cond),psychometric)


p<-ggplot() + facet_grid(.~subject)

p<-p+geom_ribbon(data=psychometrics,aes(x=phase,ymin=response.inf, ymax=response.sup,fill=cond),

                 alpha=0.3)

p<-p+geom_line(data=psychometrics,aes(x=phase, y=response, colour=cond))

p<-p+labs(x="Lag (deg)", y="Proportion",colour="Condition",fill="Condition")

p<-p+coord_cartesian(xlim = c(-0.2, 1.5))

p

4) Collapsing data across speed (the confidence intervals are plotted):