# Figure caption: proportionortion of trials, out of 50, on which # light flashes were perceived by subject S.S. as a function of log10 # intensity, together with fits. Data from Hecht et al (first seris # of trials) are shown as circles. Dashed line is the fit obtained # by linear regression. Solid curve is the fit obtained by logistic # regression. expit <- function(x) { exp(x) / (1 + exp(x)) } responses <- c(0, 2, 9, 27, 47, 50) noresponses <- 50 - responses proportion <- responses / 50 percept <- cbind(responses, noresponses) intensity<- log(c(24.1, 37.6, 58.6, 91, 141.9, 221.3), base = 10) subject1.df <- data.frame(intensity,percept) subject1.lm <- lm(proportion ~ intensity) subject1.glm <- glm(percept ~ intensity, family = binomial, data = subject1.df) x <- seq(1, 2.5, by = .01) linearRegressionLine <- subject1.lm$coef[1] + subject1.lm$coef[2]*x logisticRegressionLine <- expit(subject1.glm$coef[1]+subject1.glm$coef[2]*x) png("../figures/hechtSS.png") par(cex.lab=1.75) par(cex.axis=1.5) par(las=1) plot(intensity, proportion, xlim = c(1, 2.5), ylim = c(0, 1), xlab = "log intensity", ylab = "proportion. perceived") lines(x, linearRegressionLine, lwd=3, lty=3) lines(x, logisticRegressionLine, lwd=3) dev.off() # The analysis below is not necessary for creating the figure. summary(subject1.lm) summary(subject1.glm) b0<-subject1.glm$coef[1] b1<-subject1.glm$coef[2] vmatr<-summary(subject1.glm)$cov.unscaled dg1<- -1/b1 dg2<- b0/b1^2 x50<- -b0/b1 se<-sqrt(crossprod(c(dg1,dg2), vmatr %*% c(dg1, dg2))) library(MASS) # needed for mvrnorm beta<- mvrnorm(10000, c(b0,b1), vmatr) x50vec<- -beta[,1] / beta[,2] sqrt(var(x50vec))