1 Data Management

Input & display data

dta <- read.table("C:/Users/ASUS/Desktop/data/brightness.txt", h=T)
head(dta)
##   Stimulus Intensity Trial Yes
## 1        1        10    49   2
## 2        2        20    48   3
## 3        3        30    48  13
## 4        4        40    50  31
## 5        5        50    47  38
## 6        6        60    49  48

Add a new variable Prop defining as the % of yeas in Trial

# observed proportions
dta$Prop <-  dta$Yes/dta$Trial

2 Analysis

2.1 m0

Model Specification
Define m0 as a model using Intensity to predict y^ [y^= combining yes & No(Trial-yes)],
the dependent variable is now binary rather than continuous
Interpretation: The output showed that with the increase in stimulus intensity, it is more likely to increase correct discrimination

# use glm to fit a psychophysical function
summary(m0 <- glm(cbind(Yes, Trial - Yes) ~ Intensity, data = dta, 
                  family = binomial(link = "probit")))
## 
## Call:
## glm(formula = cbind(Yes, Trial - Yes) ~ Intensity, family = binomial(link = "probit"), 
##     data = dta)
## 
## Deviance Residuals: 
##       1        2        3        4        5        6  
##  1.0388  -0.7613  -0.3124   0.4486  -0.5979   0.7050  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -2.874139   0.286787  -10.02   <2e-16 ***
## Intensity    0.077472   0.007367   10.52   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 184.5899  on 5  degrees of freedom
## Residual deviance:   2.8119  on 4  degrees of freedom
## AIC: 26.549
## 
## Number of Fisher Scoring iterations: 4

Add a new variable (newxval) derived from Intensity, definig its score range as (5,70) with 0.5 as the unit of change.

# new values in the intensity range 
newxval <- with(dta, data.frame(Intensity = seq(5,70, by = 0.5)))

Using predict() function to generate y^ values (phat) with x=newxval, phat is defined
as the correct discrimination % (y^) predicted by x =newxal (intensity range between 5:70)

# predicted probabilities
phat <- predict(m0, newdata = newxval, type = "response" )

3 Visdualization

Plot the relationship between newxval(as Intensity) and phat (the predicted y^ values
with x=newxval)

library(ggplot2)
ot <- theme_set(theme_bw())
# fitted psychophysical curve
ggplot() +
 geom_line(aes(x=newxval[,1], y = phat)) +
 geom_point(aes(x = dta$Intensity, y = dta$Prop)) +
 geom_vline(xintercept = 35, lty = 2, col = "gray") +
 geom_hline(yintercept = .5, lty = 2, col = "gray") +
 labs(x = "Intensity", y = "Probability of correnct responses")

# The end