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