Data processing

Before beginning analyses, we will load the raw data and do some light processing, including: filtering out bad trials, computing summary statistics and sanity checking.

  1. Read the raw trial data from a CSV file. The dataset contains nearly 5,000 trials completed by three subjects.
# setwd("~/R/Schmidt-Boehm-Tuten-Roorda_2019")
library(tidyverse)
library(ggthemes)
trial_data <- read_csv('all_trial_data.csv')
# The column isPair is meant to be a logical.
trial_data <- mutate(trial_data, isPair=as.logical(isPair))
summary(select(trial_data, yb, gr, delivery_error, lConeNeighbors, isPair))
       yb                gr          delivery_error    lConeNeighbors     isPair       
 Min.   :-0.8333   Min.   :-1.0000   Min.   : 0.0000   Min.   :0.1667   Mode :logical  
 1st Qu.: 0.0000   1st Qu.:-0.2000   1st Qu.: 0.1618   1st Qu.:0.4545   FALSE:2487     
 Median : 0.0000   Median : 0.0000   Median : 0.1848   Median :0.5455   TRUE :2481     
 Mean   :-0.0120   Mean   :-0.0626   Mean   : 0.2310   Mean   :0.5723                  
 3rd Qu.: 0.0000   3rd Qu.: 0.2000   3rd Qu.: 0.2120   3rd Qu.:0.6667                  
 Max.   : 0.6000   Max.   : 1.0000   Max.   :10.1016   Max.   :1.0000                  
 NA's   :730       NA's   :730       NA's   :41        NA's   :1224                    
  1. Trials with delivery errors greater than 0.35 or less than 0.01 arcmin (values below 0.01 do not occur naturally) were considered bad deliveries. In those trials, we cannot be confident that the correct cone was targeted. After removing bad trials (3.6%), 4,788 trials remained for further analysis. Below, we see the distribution of delivery errors was roughly normal with a mean of 0.19 arcmin and a standard deviation of 0.036 arcmin.
trial_data <- trial_data %>%
  filter(., delivery_error < 0.35) %>%
  filter(., delivery_error > 0.01)

summarize(trial_data, mean=mean(delivery_error), 
          std=sd(delivery_error), N=length(delivery_error))

ggplot(trial_data, aes(delivery_error)) + 
  geom_histogram(bins = 30)
  1. Frequency of seeing (FoS) was computed for S- and L/M-cones and serves as a sanity check. FoS should fall close to zero when an S-cone was targeted because they are insensitive to the stimulus. When an L/M-cone was targeted FoS should have been close to 0.85 since the intensity of each flash was scaled to achieve 85% FoS via a detection task. Below we see that these expectations were approximately borne out. Trials that either targeted an S-cone or were not detected were not analyzed further.
Scones <- filter(trial_data, type1 == 1 | type2 == 1)
Scones <- filter(Scones, isPair == FALSE)
print(paste(c('S-cone FoS: ', 
              round(sum(!is.na(Scones$yb)) / length(Scones$delivery_error) * 100, 1)), 
            collapse = ' '))
[1] "S-cone FoS:  25.7"
print(paste(c('L/M-cone FoS: ', 
              round(sum(!is.na(trial_data$yb)) / length(trial_data$delivery_error) * 100, 1)), 
            collapse = ' '))
[1] "L/M-cone FoS:  85.5"
# Now filter out bad trials
trial_data <- filter(trial_data, type1 != 1 & type2 != 1)
trial_data <- filter(trial_data, !is.na(yb))
  1. The remaining dataset contained trials in which individual or pairs of L- and M-cones were stimulated (N=4,057). For some analyses, we wanted to know both the type tested and whether it was a pair or a single cone. To facilitate these analyses, a type ID was created in which unknown-cone=0, unknown-pair=1, 1M=2, 1L=3, 2M=4, L+M=5 and 2L=6. Note: S-cones were already removed from the data and the mosaic of S20075 has not been classified.
trial_data <- trial_data %>%
  mutate(., typeID=if_else(isPair, type1 + type2, type1)) %>%
  mutate(., typeID=if_else(typeID == 0, 
                           if_else(isPair, as.integer(1), as.integer(0)), as.integer(typeID))) %>%
  mutate(., typeIDstr=factor(typeID, labels = 
                               c('unknown-cone', 'unknown-pair', 'M-cone', 
                                 'L-cone', 'M+M-pair', 'L+M-pair', 'L+L-pair')))
tally(trial_data)
  1. Finally, saturation was computed from a sum of the absolute values of the green-red (gr) and yellow-blue (yb) dimensions (\(\mid yb \mid + \mid gr \mid\)). Below we see a plot of saturation as a function of type ID. Saturation was similar across all cone types, with exception of the unidentified cones in S20075. This subject reported more saturated sensations than the other two subjects.
trial_data <- mutate(trial_data, saturation = abs(yb) + abs(gr))
ggplot(trial_data, aes(x=saturation)) + 
  geom_histogram(binwidth=0.1) + facet_grid(. ~ typeIDstr)

Before saving the data, we visualize each column of the dataset with a histogram to get a better sense of the data. These plots illustrate the distribution of each variable and help make sure there are no outliers or unexpected values. As we can see some of the variables are categorical and some are highly non-gaussian. For instance, distance_btwn_cones_arcmin has a large peak at 0 due to single cone trials.

trial_data %>%
  select(., -x1, -x2, -y1, -y2, -session, -trial, -type1, -type2, -masterID1, -masterID2) %>%
  keep(is.numeric) %>% 
  gather() %>% 
  ggplot(aes(value)) +
    facet_wrap(~ key, scales = "free") +
    geom_histogram()

Finally, we save the filtered dataset for subsequent analysis.

# Save the data for analyses.
write_csv(trial_data, 'filtered_data.csv')
LS0tCnRpdGxlOiAnIFNwYXRpYWwgc3VtbWF0aW9uIG9mIGluZGl2aWR1YWwgY29uZXMgaW4gaHVtYW4gY29sb3IgdmlzaW9uJwphdXRob3I6ICdbQnJpYW4gUC4gU2NobWlkdF0oaHR0cHM6Ly9icHMxMC5naXRodWIuaW8pLCBBbGV4YW5kcmEgRS4gQm9laG0sIFdpbGxpYW0gUy4gVHV0ZW4sIEF1c3RpbiBSb29yZGEnCmRhdGU6ICJgciBTeXMuRGF0ZSgpYCIKb3V0cHV0OgogIGh0bWxfbm90ZWJvb2s6IGRlZmF1bHQKICBodG1sX2RvY3VtZW50OgogICAgZGZfcHJpbnQ6IHBhZ2VkCiAgcGRmX2RvY3VtZW50OiBkZWZhdWx0Ci0tLQoKIyBEYXRhIHByb2Nlc3NpbmcKCkJlZm9yZSBiZWdpbm5pbmcgYW5hbHlzZXMsIHdlIHdpbGwgbG9hZCB0aGUgcmF3IGRhdGEgYW5kIGRvIHNvbWUgbGlnaHQgcHJvY2Vzc2luZywgaW5jbHVkaW5nOiBmaWx0ZXJpbmcgb3V0IGJhZCB0cmlhbHMsIGNvbXB1dGluZyBzdW1tYXJ5IHN0YXRpc3RpY3MgYW5kIHNhbml0eSBjaGVja2luZy4gCgoxLiBSZWFkIHRoZSByYXcgdHJpYWwgZGF0YSBmcm9tIGEgQ1NWIGZpbGUuIFRoZSBkYXRhc2V0IGNvbnRhaW5zIG5lYXJseSA1LDAwMCB0cmlhbHMgY29tcGxldGVkIGJ5IHRocmVlIHN1YmplY3RzLgoKYGBge3IgcmVzdWx0cz0iaGlkZSIsIG1lc3NhZ2U9RkFMU0V9CiMgc2V0d2QoIn4vUi9TY2htaWR0LUJvZWhtLVR1dGVuLVJvb3JkYV8yMDE5IikKbGlicmFyeSh0aWR5dmVyc2UpCmxpYnJhcnkoZ2d0aGVtZXMpCmBgYAoKYGBge3IsIG1lc3NhZ2U9RkFMU0V9CnRyaWFsX2RhdGEgPC0gcmVhZF9jc3YoJ2FsbF90cmlhbF9kYXRhLmNzdicpCiMgVGhlIGNvbHVtbiBpc1BhaXIgaXMgbWVhbnQgdG8gYmUgYSBsb2dpY2FsLgp0cmlhbF9kYXRhIDwtIG11dGF0ZSh0cmlhbF9kYXRhLCBpc1BhaXI9YXMubG9naWNhbChpc1BhaXIpKQpzdW1tYXJ5KHNlbGVjdCh0cmlhbF9kYXRhLCB5YiwgZ3IsIGRlbGl2ZXJ5X2Vycm9yLCBsQ29uZU5laWdoYm9ycywgaXNQYWlyKSkKYGBgCgoyLiBUcmlhbHMgd2l0aCBkZWxpdmVyeSBlcnJvcnMgZ3JlYXRlciB0aGFuIDAuMzUgb3IgbGVzcyB0aGFuIDAuMDEgYXJjbWluICh2YWx1ZXMgYmVsb3cgMC4wMSBkbyBub3Qgb2NjdXIgbmF0dXJhbGx5KSB3ZXJlIGNvbnNpZGVyZWQgYmFkIGRlbGl2ZXJpZXMuIEluIHRob3NlIHRyaWFscywgd2UgY2Fubm90IGJlIGNvbmZpZGVudCB0aGF0IHRoZSBjb3JyZWN0IGNvbmUgd2FzIHRhcmdldGVkLiBBZnRlciByZW1vdmluZyBiYWQgdHJpYWxzICgzLjYlKSwgNCw3ODggdHJpYWxzIHJlbWFpbmVkIGZvciBmdXJ0aGVyIGFuYWx5c2lzLiBCZWxvdywgd2Ugc2VlIHRoZSBkaXN0cmlidXRpb24gb2YgZGVsaXZlcnkgZXJyb3JzIHdhcyByb3VnaGx5IG5vcm1hbCB3aXRoIGEgbWVhbiBvZiAwLjE5IGFyY21pbiBhbmQgYSBzdGFuZGFyZCBkZXZpYXRpb24gb2YgMC4wMzYgYXJjbWluLiAKCmBgYHtyIGRlbGl2ZXJ5LCAsIGZpZy5oZWlnaHQ9MiwgZmlnLndpZHRoPTR9CnRyaWFsX2RhdGEgPC0gdHJpYWxfZGF0YSAlPiUKICBmaWx0ZXIoLiwgZGVsaXZlcnlfZXJyb3IgPCAwLjM1KSAlPiUKICBmaWx0ZXIoLiwgZGVsaXZlcnlfZXJyb3IgPiAwLjAxKQoKc3VtbWFyaXplKHRyaWFsX2RhdGEsIG1lYW49bWVhbihkZWxpdmVyeV9lcnJvciksIAogICAgICAgICAgc3RkPXNkKGRlbGl2ZXJ5X2Vycm9yKSwgTj1sZW5ndGgoZGVsaXZlcnlfZXJyb3IpKQoKZ2dwbG90KHRyaWFsX2RhdGEsIGFlcyhkZWxpdmVyeV9lcnJvcikpICsgCiAgZ2VvbV9oaXN0b2dyYW0oYmlucyA9IDMwKQpgYGAKCjMuIEZyZXF1ZW5jeSBvZiBzZWVpbmcgKEZvUykgd2FzIGNvbXB1dGVkIGZvciBTLSBhbmQgTC9NLWNvbmVzIGFuZCBzZXJ2ZXMgYXMgYSBzYW5pdHkgY2hlY2suIEZvUyBzaG91bGQgZmFsbCBjbG9zZSB0byB6ZXJvIHdoZW4gYW4gUy1jb25lIHdhcyB0YXJnZXRlZCBiZWNhdXNlIHRoZXkgYXJlIGluc2Vuc2l0aXZlIHRvIHRoZSBzdGltdWx1cy4gV2hlbiBhbiBML00tY29uZSB3YXMgdGFyZ2V0ZWQgRm9TIHNob3VsZCBoYXZlIGJlZW4gY2xvc2UgdG8gMC44NSBzaW5jZSB0aGUgaW50ZW5zaXR5IG9mIGVhY2ggZmxhc2ggd2FzIHNjYWxlZCB0byBhY2hpZXZlIDg1JSBGb1MgdmlhIGEgZGV0ZWN0aW9uIHRhc2suIEJlbG93IHdlIHNlZSB0aGF0IHRoZXNlIGV4cGVjdGF0aW9ucyB3ZXJlIGFwcHJveGltYXRlbHkgYm9ybmUgb3V0LiBUcmlhbHMgdGhhdCBlaXRoZXIgdGFyZ2V0ZWQgYW4gUy1jb25lIG9yIHdlcmUgbm90IGRldGVjdGVkIHdlcmUgbm90IGFuYWx5emVkIGZ1cnRoZXIuIAoKYGBge3J9ClNjb25lcyA8LSBmaWx0ZXIodHJpYWxfZGF0YSwgdHlwZTEgPT0gMSB8IHR5cGUyID09IDEpClNjb25lcyA8LSBmaWx0ZXIoU2NvbmVzLCBpc1BhaXIgPT0gRkFMU0UpCgpwcmludChwYXN0ZShjKCdTLWNvbmUgRm9TOiAnLCAKICAgICAgICAgICAgICByb3VuZChzdW0oIWlzLm5hKFNjb25lcyR5YikpIC8gbGVuZ3RoKFNjb25lcyRkZWxpdmVyeV9lcnJvcikgKiAxMDAsIDEpKSwgCiAgICAgICAgICAgIGNvbGxhcHNlID0gJyAnKSkKcHJpbnQocGFzdGUoYygnTC9NLWNvbmUgRm9TOiAnLCAKICAgICAgICAgICAgICByb3VuZChzdW0oIWlzLm5hKHRyaWFsX2RhdGEkeWIpKSAvIGxlbmd0aCh0cmlhbF9kYXRhJGRlbGl2ZXJ5X2Vycm9yKSAqIDEwMCwgMSkpLCAKICAgICAgICAgICAgY29sbGFwc2UgPSAnICcpKQoKIyBOb3cgZmlsdGVyIG91dCBiYWQgdHJpYWxzCnRyaWFsX2RhdGEgPC0gZmlsdGVyKHRyaWFsX2RhdGEsIHR5cGUxICE9IDEgJiB0eXBlMiAhPSAxKQp0cmlhbF9kYXRhIDwtIGZpbHRlcih0cmlhbF9kYXRhLCAhaXMubmEoeWIpKQpgYGAKCjQuIFRoZSByZW1haW5pbmcgZGF0YXNldCBjb250YWluZWQgdHJpYWxzIGluIHdoaWNoIGluZGl2aWR1YWwgb3IgcGFpcnMgb2YgTC0gYW5kIE0tY29uZXMgd2VyZSBzdGltdWxhdGVkIChOPTQsMDU3KS4gRm9yIHNvbWUgYW5hbHlzZXMsIHdlIHdhbnRlZCB0byBrbm93IGJvdGggdGhlIHR5cGUgdGVzdGVkIGFuZCB3aGV0aGVyIGl0IHdhcyBhIHBhaXIgb3IgYSBzaW5nbGUgY29uZS4gVG8gZmFjaWxpdGF0ZSB0aGVzZSBhbmFseXNlcywgYSB0eXBlIElEIHdhcyBjcmVhdGVkIGluIHdoaWNoIHVua25vd24tY29uZT0wLCB1bmtub3duLXBhaXI9MSwgMU09MiwgMUw9MywgMk09NCwgTCtNPTUgYW5kIDJMPTYuIE5vdGU6IFMtY29uZXMgd2VyZSBhbHJlYWR5IHJlbW92ZWQgZnJvbSB0aGUgZGF0YSBhbmQgdGhlIG1vc2FpYyBvZiBTMjAwNzUgaGFzIG5vdCBiZWVuIGNsYXNzaWZpZWQuCgpgYGB7cn0KdHJpYWxfZGF0YSA8LSB0cmlhbF9kYXRhICU+JQogIG11dGF0ZSguLCB0eXBlSUQ9aWZfZWxzZShpc1BhaXIsIHR5cGUxICsgdHlwZTIsIHR5cGUxKSkgJT4lCiAgbXV0YXRlKC4sIHR5cGVJRD1pZl9lbHNlKHR5cGVJRCA9PSAwLCAKICAgICAgICAgICAgICAgICAgICAgICAgICAgaWZfZWxzZShpc1BhaXIsIGFzLmludGVnZXIoMSksIGFzLmludGVnZXIoMCkpLCBhcy5pbnRlZ2VyKHR5cGVJRCkpKSAlPiUKICBtdXRhdGUoLiwgdHlwZUlEc3RyPWZhY3Rvcih0eXBlSUQsIGxhYmVscyA9IAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgYygndW5rbm93bi1jb25lJywgJ3Vua25vd24tcGFpcicsICdNLWNvbmUnLCAKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgJ0wtY29uZScsICdNK00tcGFpcicsICdMK00tcGFpcicsICdMK0wtcGFpcicpKSkKCnRhbGx5KHRyaWFsX2RhdGEpCmBgYAoKNS4gRmluYWxseSwgc2F0dXJhdGlvbiB3YXMgY29tcHV0ZWQgZnJvbSBhIHN1bSBvZiB0aGUgYWJzb2x1dGUgdmFsdWVzIG9mIHRoZSBncmVlbi1yZWQgKGdyKSBhbmQgeWVsbG93LWJsdWUgKHliKSBkaW1lbnNpb25zICgkXG1pZCB5YiBcbWlkICsgXG1pZCBnciBcbWlkJCkuIEJlbG93IHdlIHNlZSBhIHBsb3Qgb2Ygc2F0dXJhdGlvbiBhcyBhIGZ1bmN0aW9uIG9mIHR5cGUgSUQuIFNhdHVyYXRpb24gd2FzIHNpbWlsYXIgYWNyb3NzIGFsbCBjb25lIHR5cGVzLCB3aXRoIGV4Y2VwdGlvbiBvZiB0aGUgdW5pZGVudGlmaWVkIGNvbmVzIGluIFMyMDA3NS4gVGhpcyBzdWJqZWN0IHJlcG9ydGVkIG1vcmUgc2F0dXJhdGVkIHNlbnNhdGlvbnMgdGhhbiB0aGUgb3RoZXIgdHdvIHN1YmplY3RzLgoKYGBge3IsIGZpZy5oZWlnaHQ9MywgZmlnLndpZHRoPTksIG1lc3NhZ2U9RkFMU0V9CnRyaWFsX2RhdGEgPC0gbXV0YXRlKHRyaWFsX2RhdGEsIHNhdHVyYXRpb24gPSBhYnMoeWIpICsgYWJzKGdyKSkKCmdncGxvdCh0cmlhbF9kYXRhLCBhZXMoeD1zYXR1cmF0aW9uKSkgKyAKICBnZW9tX2hpc3RvZ3JhbShiaW53aWR0aD0wLjEpICsgZmFjZXRfZ3JpZCguIH4gdHlwZUlEc3RyKQoKYGBgCgpCZWZvcmUgc2F2aW5nIHRoZSBkYXRhLCB3ZSB2aXN1YWxpemUgZWFjaCBjb2x1bW4gb2YgdGhlIGRhdGFzZXQgd2l0aCBhIGhpc3RvZ3JhbSB0byBnZXQgYSBiZXR0ZXIgc2Vuc2Ugb2YgdGhlIGRhdGEuIFRoZXNlIHBsb3RzIGlsbHVzdHJhdGUgdGhlIGRpc3RyaWJ1dGlvbiBvZiBlYWNoIHZhcmlhYmxlIGFuZCBoZWxwIG1ha2Ugc3VyZSB0aGVyZSBhcmUgbm8gb3V0bGllcnMgb3IgdW5leHBlY3RlZCB2YWx1ZXMuIEFzIHdlIGNhbiBzZWUgc29tZSBvZiB0aGUgdmFyaWFibGVzIGFyZSBjYXRlZ29yaWNhbCBhbmQgc29tZSBhcmUgaGlnaGx5IG5vbi1nYXVzc2lhbi4gRm9yIGluc3RhbmNlLCBgZGlzdGFuY2VfYnR3bl9jb25lc19hcmNtaW5gIGhhcyBhIGxhcmdlIHBlYWsgYXQgMCBkdWUgdG8gc2luZ2xlIGNvbmUgdHJpYWxzLgoKYGBge3IsIHdhcm5pbmc9RkFMU0UsbWVzc2FnZT1GQUxTRX0KCnRyaWFsX2RhdGEgJT4lCiAgc2VsZWN0KC4sIC14MSwgLXgyLCAteTEsIC15MiwgLXNlc3Npb24sIC10cmlhbCwgLXR5cGUxLCAtdHlwZTIsIC1tYXN0ZXJJRDEsIC1tYXN0ZXJJRDIpICU+JQogIGtlZXAoaXMubnVtZXJpYykgJT4lIAogIGdhdGhlcigpICU+JSAKICBnZ3Bsb3QoYWVzKHZhbHVlKSkgKwogICAgZmFjZXRfd3JhcCh+IGtleSwgc2NhbGVzID0gImZyZWUiKSArCiAgICBnZW9tX2hpc3RvZ3JhbSgpCgpgYGAKCkZpbmFsbHksIHdlIHNhdmUgdGhlIGZpbHRlcmVkIGRhdGFzZXQgZm9yIHN1YnNlcXVlbnQgYW5hbHlzaXMuCgpgYGB7cn0KIyBTYXZlIHRoZSBkYXRhIGZvciBhbmFseXNlcy4Kd3JpdGVfY3N2KHRyaWFsX2RhdGEsICdmaWx0ZXJlZF9kYXRhLmNzdicpCgpgYGA=