This notebook shows how to analyze data obtained in an MLCM experiment using the R package MLCM. It is similar in structure as the accompaning notebook on MLDS. The notebook is divided into four parts:

  1. Loading and parsing data
  2. Estimating the perceptual scale
  3. Calculating confidence intervals
  4. Evaluating goodness of fit.

To start, we need to have the MLDS packaged installed and loading in our R session.

# uncomment to install package if you don't have it.
# you only need to do this once.

#install.packages('MLCM')

# loading package into the session
library(MLCM)

1. Loading and parsing data

Loading data

The first step in the analysis is to load the data from the experiment, and parse them to the format the functions in the MLDS package understands.

The following code loads a CSV file containing data from an MLCM experiment on brightness perception of the White’s illusion. This was an experiment programmed with Psychopy and run online with Pavlovia.

The following file contains data of the MLCM experiment run six times, for a total of 572 trials. I added it here already aggregated so that you don’t have to run so many trials on your own.

# Loading a single file
fname <- 'data/mlcm/GA.csv'

df <- read.csv(fname)
df

Now we only keep the columns relevant for our analysis. These are:

  • lum1, lum2: the luminance values presented in each trial, for the left and right stimuli
  • the “context” on which each stimulus was presented, 0 for on a black bar, 1 for on a white bar.
  • resp.keys: what button the observer pressed.
df <- df[c('resp.keys', 'lum1', 'lum2', 'c1', 'c2')]
df

Parsing data to format needed by MLCM

We need to map the key presses to 0’s or 1’s. We do this with the following call and save the responses in a new column called resp.

df$resp <- as.numeric(df$resp.keys == 'left')
df

Now we remove NaNs entered automatically by Psychopy

df <- na.omit(df)
df

Then we need to put the data in the format that the MLCM package likes. Each row is a trial, and we will have 5 columns:

  • Resp: containing the binary response (0 or 1)
  • L1, L2: containing the stimulus indices along the first stimulus dimension, in this case luminance.
  • C1, C2: containing the stimulus indices along the second stimulus dimension, in this case the “context” (on a black or white bar)

The following code rearranges the data to that format, saving it in variable data

lum_vector <- c(0.25, 0.33, 0.42, 0.5, 0.58, 0.67,0.75)

lum_index <- seq(length(lum_vector))

for (i in 1:length(lum_vector)) {
  df$lum1 <- replace(df$lum1, df$lum1==lum_vector[i], lum_index[i])
  df$lum2 <- replace(df$lum2, df$lum2==lum_vector[i], lum_index[i])
}

# we will map context 0 to index 2, and context 1 to index 1. 
# this last one doesn't change, so we only need to change from 0 --> 2
df$c1 <- replace(df$c1, df$c1==0, 2)
df$c2 <- replace(df$c2, df$c2==0, 2)

data <- df[c('resp', 'lum1', 'lum2', 'c1', 'c2')]

# renaming columns
colnames(data) <- c("Resp", "L1", "L2", "C1", "C2")

data

This last data frame, data, is in the correct format ready for MLDS.

2. Estimating the perceptual scales

Scales can be fitted using one of three different decisions models in MLCM. The most general one is the ‘full’ model, in which MLCM finds a scale value to each stimulus.

Then there is the ‘additive’ model, which assumes that the effect of the stimulus dimension is additive. It finds a scale for one dimension (luminance), and for the other dimension (context) only finds an additive factor.

Finally the most constrained model is the ‘independent’ model, where it is assumed that only one stimulus dimension determines appearance.

Here we fit the three models, then compared them to determine which one is the most appropriate one, and then plot the scales for that model only. Finally we calculate confidence intervals for them.

‘Full’ model

# full model
scales_full <- mlcm(data, model='full')
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
scales_full
## 
## Maximum Likelihood Conjoint Measurement
## 
## Model:    Full
## Perceptual Scale:
##    C1    C2   
## L1  0.00  7.65
## L2  7.87  9.60
## L3  8.83 10.71
## L4 10.38 12.02
## L5 10.50 13.21
## L6 12.00 14.20
## L7 13.25 21.74
## 
## sigma:
## [1] 1

‘Additive’ and ‘independent’ model

# additive
scales_add <- mlcm(data, model='add')
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
scales_add
## 
## Maximum Likelihood Conjoint Measurement
## 
## Model:    Additive
## Perceptual Scale:
##      L    C   
## Lev1 0.00 0.00
## Lev2 2.20 2.06
## Lev3 3.32   NA
## Lev4 4.80   NA
## Lev5 5.40   NA
## Lev6 6.59   NA
## Lev7 7.81   NA
## 
## sigma:
## [1] 1
# independent
scales_ind <- mlcm(data, model='ind', whichdim = 1)
scales_ind
## 
## Maximum Likelihood Conjoint Measurement
## 
## Model:    Independence
## Perceptual Scale:
##    [,1]
## L1 0.00
## L2 1.20
## L3 1.77
## L4 2.48
## L5 2.83
## L6 3.45
## L7 4.02
## 
## sigma:
## [1] 1

Comparing models

The likelihood ratio test is used for model comparison. For this statistical test the statistic follows the Chi-squared distribution.

# comparing models
anova(scales_ind, scales_add, scales_full, test = "Chisq")

For these data the model comparison tells us that the additive model is signifficantly better than the full model in explaining the data (see asterisks next to p-value). For the rest of the notebook we then only keep analysing the additive model.

# change here to 'full' if in the comparison above the additive model is not 
# sig. better than the full
modeltype = 'add' # 'full'
if (modeltype=='add'){
  scales <- scales_add
  
} else if(modeltype=='full'){
  scales <- scales_full
}

3. Calculating confidence intervals

Calculating confidence intervals follows the same logic as for MLDS. Refer to the MLDS tutorial for details.

scalesboot <- boot.mlcm(scales, nsim=1000)

## here  we manually calculate the percentile confidence intervals from the bootstrap samples
samples <- scalesboot$boot.samp

## if the model is additive, we need to rearrange the bootstrap samples a bit to get a full matrix of samples
if (modeltype=='add'){
  additiveshift <- samples[nrow(samples),]
  
  context1 <- samples[1:nrow(samples)-1,]
  context2 <- rbind(rep(0, ncol(samples)), samples[1:nrow(samples)-1,]) + additiveshift
  
  samples <- rbind(context1, context2)
}

mns <- apply(samples, 1, mean)  # bootstrap mean
#ci95 <- qnorm(0.975) * apply(samples, 1, sd) 

# low and high bounds of confidence interval
low <- apply(samples, 1, quantile, probs = 0.025)
high <- apply(samples, 1, quantile, probs = 0.975)

Plotting scales with CIs

par(pty="s")
# plot one scale for first context
plot(lum_vector, scales$pscale[,1], 
     xlab = "Luminance [relative]", ylab = "Perceived brightness", 
     col="darkgray",
     ylim=c(0, max(mns)*1.2), 
     #xlim=c(0, 1),
     )

# and for the second context
if (modeltype=='add'){
  # manually add the additive factor to the first scale
  points(lum_vector, scales$pscale[,1] + scales$pscale[2,2] , col="black")
} else {
  if (modeltype=='full'){
    # not necessary to add the factor manually
    points(lum_vector,scales$pscale[,2], col="black")
  }
}

# CIs
ns <- length(lum_vector)
segments(lum_vector[2:ns], low[1:ns-1], lum_vector[2:ns], high[1:ns-1], col="darkgray")
segments(lum_vector[1:ns], low[(ns+1):(2*ns)-1], lum_vector[1:ns], high[(ns+1):(2*ns)-1], col="black")

4. Evaluating goodness of fit

The type of evaluation is identical as in MLDS. Check out the MLDS tutorial (tutorial_mlds.Rmd) and the slides for details.

scales.gof <- binom.diagnostics(scales, nsim=1000)
plot(scales.gof)

print(scales.gof$p)
## [1] 0.171