Lab 6

Author

Filippo Gambarota

devtools::load_all()
library(tidyverse)
library(ggeffects)

Overview

EDA

Loading data

data("psych")
dat <- psych

Given that we did not introduced random-effects models, we select a single subject to analyze.

cc <- seq(0, 1, 0.1)

dat$contrast_c <- cut(dat$contrast, cc, include.lowest = TRUE)

dat |> 
  filter(id %in% sample(unique(dat$id), 5)) |> 
  group_by(id, cond, contrast_c) |> 
  summarise(y = mean(y)) |> 
  ggplot(aes(x = contrast_c, y = y, color = cond, group = cond)) +
  geom_point() +
  geom_line() +
  theme(axis.text.x = element_text(angle = 90)) +
  facet_wrap(~id)

`summarise()` has grouped output by 'id', 'cond'. You can override using the
`.groups` argument.

dat <- filter(dat, id == 6)

We have several interesting stuff to estimate. Let’s start by fitting a simple model:

fit1 <- glm(y ~ contrast, data = dat, family = binomial(link = "logit"))
summary(fit1)


Call:
glm(formula = y ~ contrast, family = binomial(link = "logit"), 
    data = dat)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.8826  -0.6003   0.2343   0.5974   2.1113  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -2.11482    0.09237  -22.90   <2e-16 ***
contrast     6.41786    0.22165   28.95   <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: 4038.1  on 2999  degrees of freedom
Residual deviance: 2500.5  on 2998  degrees of freedom
AIC: 2504.5

Number of Fisher Scoring iterations: 5

The parameters (Intercept) and contrast are respectively the probability of saying yes for stimuli with 0 contrast. Seems odd but in Psychophysics this is a very interesting information. We can call it the false alarm rate. We usually expect this rate to be low, ideally 0.

plogis(coef(fit1)[1])

(Intercept) 
   0.107665

The contrast is the slope i.e. the increase in the log odds of saying yes for a unit increase in contrast. In this case this parameter is hard to intepret, let’s change the scale of the contrast:

dat$contrast10 <- dat$contrast * 10
fit1 <- glm(y ~ contrast10, data = dat, family = binomial(link = "logit"))
summary(fit1)


Call:
glm(formula = y ~ contrast10, family = binomial(link = "logit"), 
    data = dat)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.8826  -0.6003   0.2343   0.5974   2.1113  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -2.11482    0.09237  -22.90   <2e-16 ***
contrast10   0.64179    0.02217   28.95   <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: 4038.1  on 2999  degrees of freedom
Residual deviance: 2500.5  on 2998  degrees of freedom
AIC: 2504.5

Number of Fisher Scoring iterations: 5

Now the contrast10 is the increase in the log odds of saying yes for an increase of 10% contrast. Using the divide-by-4 rule we obtain an maximal increase of 0.1604465 of probability of saying yes.

Another interesting parameter is the threshold. In psychophysics the threshold is the required \(x\) level (in this case contrast) to obtain a certain proportions of \(y\) response.

For a logistic distribution (see Knoblauch and Maloney 2012) the 50% threshold can be estimated as \(-\frac{\beta_0}{\beta_1}\) thus:

-(coef(fit1)[1]/coef(fit1)[2])

(Intercept) 
   3.295206

Then the slope is simply the inverse of the regression slope and represent the increase in performance/visibility for a unit increase in \(x\):

1/coef(fit1)[2]

contrast10 
  1.558152

In fact, we can use these parameters to plot a logistic distribution:

curve(plogis(x, -(coef(fit1)[1]/coef(fit1)[2]), 1/coef(fit1)[2]),
      0, 10)

That is very similar to the effects estimated by our model:

plot(ggeffect(fit1))

$contrast10

References

Knoblauch, Kenneth, and Laurence T Maloney. 2012. Modeling Psychophysical Data in r. New York, NY: Springer New York. https://doi.org/10.1007/978-1-4614-4475-6.