I'm trying to estimate the confidence interval of a psychometric curve (binomial probit GLM), for a population (now only two subjects).
Suppose I've subject "a" and subject "b", which performs separately a task. The answers are recorded, and plotted against the stimulus level, giving us a sygmoidal curve.
Here is the data, and the code to plot the psychometric curve with its S.E.:
# data:
mydata <- structure(list(exp_dur = c(3250L, 2850L, 2250L, 2450L, 3450L,
3450L, 3050L, 3450L, 2450L, 2650L, 3250L, 2650L, 2850L, 3250L,
3050L, 2250L, 2450L, 2650L, 2250L, 2850L, 3050L, 2850L, 3050L,
2850L, 2250L, 2450L, 3450L, 3250L, 3450L, 2850L, 3050L, 2250L,
2450L, 3050L, 3450L, 2250L, 3050L, 3250L, 2450L, 2650L, 3250L,
2650L, 2850L, 2450L, 2650L, 3450L, 2250L, 3250L, 2650L, 3250L,
2450L, 2650L, 2650L, 3450L, 3450L, 3250L, 2450L, 2650L, 2850L,
3050L, 3050L, 3250L, 2250L, 3050L, 3450L, 2250L, 2450L, 2850L,
2850L, 2250L, 2250L, 3050L, 2650L, 2650L, 2450L, 3250L, 3050L,
2850L, 2850L, 3250L, 3250L, 3050L, 2450L, 3050L, 2450L, 3450L,
3250L, 2250L, 3450L, 2250L, 3450L, 2250L, 3450L, 2850L, 2650L,
2450L, 2850L, 2650L), response_key_exp_resp = structure(c(2L,
2L, 1L, 1L, 2L, 2L, 2L, 2L, 1L, 1L, 2L, 1L, 2L, 2L, 2L, 1L, 1L,
2L, 1L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 2L, 2L, 2L, 1L, 2L, 1L, 1L,
2L, 2L, 1L, 2L, 2L, 1L, 1L, 2L, 1L, 2L, 1L, 1L, 2L, 1L, 2L, 1L,
2L, 1L, 1L, 2L, 2L, 2L, 2L, 1L, 2L, 1L, 2L, 2L, 2L, 1L, 2L, 2L,
1L, 1L, 1L, 2L, 1L, 2L, 2L, 2L, 2L, 1L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 1L, 2L, 1L, 2L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 1L, 1L, 1L,
1L), .Label = c("c", "l"), class = "factor"), lcorr = c(1, 1,
0, 0, 1, 1, 1, 1, 0, 0, 1, 0, 1, 1, 1, 0, 0, 1, 0, 1, 1, 1, 1,
1, 0, 0, 1, 1, 1, 0, 1, 0, 0, 1, 1, 0, 1, 1, 0, 0, 1, 0, 1, 0,
0, 1, 0, 1, 0, 1, 0, 0, 1, 1, 1, 1, 0, 1, 0, 1, 1, 1, 0, 1, 1,
0, 0, 0, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 0, 1,
1, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 0), subj = structure(c(1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L
), .Label = c("a", "b"), class = "factor")), .Names = c("exp_dur",
"response_key_exp_resp", "lcorr", "subj"), row.names = c(NA,
98L), class = "data.frame")
For subject "a":
# subject a:
myfit.a <- glm(lcorr~exp_dur, binomial(probit), mydata, subset=subj=="a")
suba <- mydata[which(mydata$subj=="a"),]
sc1a <- with(suba, tapply(lcorr, exp_dur, mean))
plot(as.numeric(names(sc1a)), sc1a, log="x", xlab="Contrast", ylab="P correct", pch=21, cex=1.5)
cnt <- seq(2250, 3450, len=2000)
Long.Pred.a <- predict(myfit.a, newdata=data.frame(exp_dur=cnt), type="response", se.fit=T)
polygon(c(cnt, rev(cnt)), c(Long.Pred.a$fit + Long.Pred.a$se.fit, rev(Long.Pred.a$fit - Long.Pred.a$se.fit)), border="white", col="lightgrey")
lines(cnt, Long.Pred.a$fit, lwd=2)
points(as.numeric(names(sc1a)), sc1a, pch=21, cex=1.5)
For subject "b":
# subject b:
myfit.b <- glm(lcorr~exp_dur, binomial(probit), mydata, subset=subj=="b")
subb <- mydata[which(mydata$subj=="b"),]
sc1b <- with(subb, tapply(lcorr, exp_dur, mean))
plot(as.numeric(names(sc1b)), sc1b, log="x", xlab="Contrast", ylab="P correct", pch=19, cex=1.5)
cnt <- seq(2250, 3450, len=2000)
Long.Pred.b <- predict(myfit.b, newdata=data.frame(exp_dur=cnt), type="response", se.fit=T)
polygon(c(cnt, rev(cnt)), c(Long.Pred.b$fit + Long.Pred.b$se.fit, rev(Long.Pred.b$fit - Long.Pred.b$se.fit)), border="white", col="lightgrey")
lines(cnt, Long.Pred.b$fit, lwd=2, lty=2)
points(as.numeric(names(sc1b)), sc1b, pch=19, cex=1.5)
Now the psychometric curve of the group (subjects a and b) whith its plotted SE is:
# subject a+b
myfit.t <- glm(lcorr~exp_dur*subj, binomial(probit), mydata)
sc1t <- with(mydata, tapply(lcorr, exp_dur, mean))
plot(as.numeric(names(sc1t)), sc1t, log="x", xlab="Contrast", ylab="P correct", cex=1.5)
cnt <- seq(2250, 3450, len=2000)
Long.Pred.t <- predict(myfit.t, newdata=data.frame(exp_dur=cnt), type="response", se.fit=T)
polygon(c(cnt, rev(cnt)), c(Long.Pred.b$fit + Long.Pred.b$se.fit, rev(Long.Pred.b$fit - Long.Pred.b$se.fit)), border="white", col="lightgrey")
lines(cnt, Long.Pred.t$fit, lwd=2, lty=3)
points(as.numeric(names(sc1t)), sc1t, cex=1.5)
At a visual inspection, I don't think this is the confidence interval of the group. In fact, in the formula used, I did not specified that the data comes from 2 subject, and R interpret like this is a single subject.
So, what is the method to get the S.E. for the group in a GLM analysis?
I'm quite sure that the correct approach to this kind of problem is GLMM (lme4 package), but I'm not familiar with it, and it seems to be hard to plot the predict data for glmer models, which would look like this:
mod1 <- glmer(lcorr ~ exp_dur + (1 + exp_dur|subj), family = binomial(link = "probit"), data=mydata)
How can I plot the S.E. of the group' psychometric curve?
