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.
- 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
- 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)
- 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))
- 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)
- 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=