A cortical semantic space integrating fractions and integers

1. ---
2. title: Analyses ligne numérique CSEN - Expérience 1
3. author:
4. - Maxime Cauté
5. lang: fr
6. output: rmdformats::readthedown
7. fig-dpi: 160
8. ---
10. ```{r settings, echo = FALSE, message = FALSE}
11. knitr::opts_chunk$set(echo = FALSE, warning = FALSE, message = FALSE)
13. # Package manager
14. if(!require("pacman")) install.packages("pacman")
15. # GGplot-friendly library for IRT.
16. if(!require("ggmirt")) devtools::install_github("masurp/ggmirt")
17. pacman::p_load(tidyr, dplyr, tibble, magrittr, # Tidy scripts
18. apa, # Report stats in APA format
19. mirt, ggmirt, # IRT analysis
20. egg, # Very neat ggarrange
21. ggplot2,
22. devtools,
23. ggprism, grid, ggpubr, # Plot legends
24. purrr,
25. officer, rvg, # Outputing pptx
26. pandoc, # HTML output
27. stringr, #String manipulations
28. #cooccur,
29. reshape, # Matrix manipulation and stuff
30. install = T)
32. ```
34. ```{r}
35. # All helpers in a dedicated file
36. source("all_helpers.R")
38. # Parameters for printing
39. PT_TO_PIXEL = 16/12
40. FAMILY = "NimbusSan"
42. # Add a theme that specifies the size
43. theme_text = function(size = 16, title.size = 18, family = FAMILY, scaling_factor = 1){
44. theme(text = element_text(family = family, size = size * scaling_factor),
45. axis.text = element_text(family = family, size = size * scaling_factor),
46. legend.text = element_text(family = family, size = size * scaling_factor),
47. axis.title = element_text(family = family, size = title.size * scaling_factor),
48. legend.title = element_text(family = family, size = title.size * scaling_factor),
49. strip.text = element_text(family = family, size = title.size * scaling_factor),
50. plot.subtitle = element_text(family = family, size = title.size * scaling_factor),
51. plot.title = element_text(family = family, size = title.size * scaling_factor)) %>%
52. return
53. }
55. # Use PT_TO_PIXEL to save the text as the correct size in pt in svg (svg saves raw text.size as px)
56. # If you want to export as knit, you may want to change this.
57. theme_set(theme_get() %+replace% theme_text(scaling_factor = PT_TO_PIXEL))
58. ```
60. # Data loading
62. Load the data, only drop the column index.
63. ```{r dataLoading}
64. d = read.csv("./data/cleaned_all_years.csv") %>% select(-X)
65. ```
67. ### Performance on integers and decimals
69. Only keep data relevant to the experiment: data from 6th graders with no tutorial.
70. ```{r dataFiltering}
71. d %<>%
72. # Remove the tutorial trial
73. filter(!is.tutorial) %>%
74. # Only keep 6th graders for Exp. 1
75. filter(class == "6ème") %>%
76. # Only keep data from the second line, on integers (6), decimals (7), and fractions (9)
77. filter(stimtype %in% c(6, 7, 9))
78. ```
80. Check that we have what we expect: average error rates with examples for each type of stimulus
81. ```{r performance on integers and decimals}
82. d %>%
83. # Compute the error rate for each type of stimulus, for each participant. Take the first stimulus as an example.
84. group_by(subj, stimtype) %>% summarise(
85. error_rate = mean(!is.correct), sample_stimulus = stimulus %>% first, .groups = "drop") %>%
86. # Average across targets
87. group_by(stimtype) %>% summarise(
88. error_rate.se = se(error_rate), error_rate = mean(error_rate),
89. sample_stimulus = sample_stimulus %>% first) %>%
90. # Round error rates at 3 decimals, and error rate at 4. Display both as %.
91. mutate(error_rate = error_rate %>% (scales::percent_format(accuracy = 0.1)),
92. error_rate.se = error_rate.se %>% (scales::percent_format(accuracy = 0.01))) %>%
93. # Only display relevent informations
94. select(stimtype, sample_stimulus, error_rate, error_rate.se)
95. ```
96. As we can see, performance on integers is correct, somewhat so for decimals. Fractions really stand out!
98. ### Fractions
100. ```{r fractionFilter}
101. # Now that we looked at error rates across all targets, drop everything but fractions.
102. d %<>% filter(stimtype == 9) %>% select(-stimtype)
103. ```
105. ```{r dataFormattingHelpers}
106. # Create pretty names for our different categories of fractions
107. fraction_categories = c("Equivalent to 1/2", "Integer result", "n / 1", "Tenths")
109. # Create a function that categorizes fractions...
110. categorize_fractions = function(numerator, denominator){
111. factor(
112. # ... if the denominator is ten.
113. if_else(denominator == 10, fraction_categories[4],
114. # ... if the denominator is 1
115. if_else(denominator == 1, fraction_categories[3],
116. # ... if the fraction is equivalent to 1/2
117. if_else(round(numerator / denominator, 1) == 0.5, fraction_categories[1],
118. # ... or else, then the fraction is equivalent to an integer
119. fraction_categories[2]))),
120. # Make this a factor with the desired order.
121. levels = fraction_categories) %>%
122. return
123. }
125. # Define a helper function that counts responses with an entry for all possible responses,
126. # including those that were never answered
127. # This is not done by count().
128. count_possible_responses = function(d, ...){
129. d %>%
130. # Count responses per stimulus and position and whatever additional argument is passed.
131. count(stimulus, actual_response, ...) %>%
132. # Right join with...
133. right_join(
134. # ... all the target stimuli (and whatever is also grouping)...
135. d %>% select(stimulus, ...) %>% unique %>%
136. # ... and all 51 positions on the number line
137. merge(tibble(position = seq(0, 5, 0.1) %>% round(1))) %>%
138. # And now order that to make it tidy
139. arrange(stimulus),
140. # Join by the specified columns, mapping actual response positions to the possible positions
141. by = d %>% select(stimulus, ...) %>% colnames %>% c("actual_response" = "position")
142. ) %>%
143. # If some stimulus was added with the right join, then it has no responses.
144. mutate(n = n %>% replace_na(0)) %>%
145. # Make everything easier to read.
146. arrange(stimulus, actual_response) %>%
147. return
148. }
150. # This loads a function that can compute_fractions_predictions, which is shared across analysis files.
151. source("src/fractions_predictions.R")
153. # Compute the predictions for our actual stimuli.
154. fractions_predictions.d = d %>%
155. # Select the relevant stimuli and make them unique
156. select(stimulus) %>% unique %>%
157. # Compute the predictions.
158. compute_fractions_predictions
159. ```
162. ```{r dataFormatting}
163. # Format the data to help later computations
164. d %<>%
165. # Get the exact actual response to avoid floating-point errors.
166. mutate(actual_response = actual.response %>% round(1), .keep = "unused") %>%
167. # Extract the components of the stimuli
168. mutate(numerator = stimulus %>% str_split_i(" / ", 1) %>% as.numeric,
169. denominator = stimulus %>% str_split_i(" / ", 1) %>% as.numeric) %>%
170. # Compute the fraction category from the components
171. mutate(category = categorize_fractions(numerator, denominator)) %>%
172. # Adopt a snake_case convention.
173. dplyr::rename(is_correct = is.correct)
174. ```
176. # Demographics
178. Number of girls
179. ```{r}
180. # Read demographics data
181. demographics.d = read.csv("./demographics/demographics_all_years.csv")
183. # Compute the number of each sex.
184. demographics.d %>%
185. # Only keep the participants whose data we care about
186. filter(subj %in% d$subj) %>%
187. # Count by sex
188. count(sex)
189. ```
190. # Poor and varied performance
192. ## Global error rate
194. How many errors do children do on average?
196. ```{r global error rate}
197. # Compute the error rate averaged across subjects
198. d %>%
199. # Compute each subject's error rate
200. group_by(subj) %>% summarise(error_rate = mean(!is_correct), .groups = "drop") %>%
201. # And now compute the average error rate across subjects, as well as standard error.
202. summarise(global_error_rate.se = se(error_rate), global_error_rate = mean(error_rate), .groups = "drop") %>%
203. (function(d){
204. cat("The average error rate in 6ème is ",
205. d$global_error_rate[[1]] %>% (scales::percent_format(accuracy = 0.1)),
206. " (se=",
207. d$global_error_rate.se[[1]] %>% (scales::percent_format(accuracy = 0.01)), ").", sep = "")
208. })
209. ```
211. ## Error rates per stimulus
213. Are some stimuli easier than others?
214. ```{r error rate per stimulus, fig.width = 8.27/2, fig.height = 11.69 / 2}
215. # Compute error rates per stimulus
216. error_rates.d = d %>%
217. # Group by stimulus and category...
218. group_by(stimulus, category) %>%
219. # ...and compute error rate and half binomial confidence interval (since we work with binary data)
220. summarise(error_rate = mean(!is_correct), error_rate.half_ci = half_binomial_ci(error_rate, n()), .groups = "drop") %>%
221. # Order by error rate to make everything easier to read.
222. arrange(-error_rate)
223. ```
225. Compute the maximum, minimum and 1/2's error rates
226. ```{r remarkableStimulusErrorRates}
227. error_rates.d %>%
228. select(stimulus, error_rate) %>%
229. filter(error_rate == max(error_rate)) %>%
230. (function(d){
231. cat("Maximum error rate:", (scales::percent_format)(accuracy = 2)(d$error_rate[[1]]),
232. "for fraction", d$stimulus[[1]])
233. cat(".\n")
234. })
236. error_rates.d %>%
237. select(stimulus, error_rate) %>%
238. filter(error_rate == min(error_rate)) %>%
239. (function(d){
240. cat("Minimum error rate:", (scales::percent_format)(accuracy = 2)(d$error_rate[[1]]),
241. "for fraction", d$stimulus[[1]])
242. cat(".\n")
243. })
245. error_rates.d %>%
246. select(stimulus, error_rate) %>%
247. filter(stimulus == "1 / 2") %>%
248. (function(d){
249. cat("Error rate for 1 / 2:", (scales::percent_format)(accuracy = 2)(d$error_rate[[1]]))
250. cat(".\n")
251. })
252. ```
254. Significance of differences
255. ```{r error rate per stimulus chi square test and anova}
256. d %>%
257. # Select each response as an entry with the stimulus and response correctness
258. select(stimulus, is_correct) %>%
259. # Turn this into a table...
260. table %>%
261. # On which we perform a chi-square test
262. chisq.test %>%
263. # Print it according to APA norms
264. apa(print = T)
265. ```
267. ## Error rates per category
269. ```{r error rates per category}
270. # Compute the average error rates for each category.
271. d %>%
272. # Compute average error rate per target
273. # (no need to average by participant since a target was only seen once at most by a given participant)
274. group_by(stimulus, category) %>% summarise(error_rate = mean(!is_correct), .groups = "drop") %>%
275. # Now average across stimuli
276. group_by(category) %>% summarise(error_rate = mean(error_rate), .groups = "drop") %>%
277. # Order by error rate to make it easier to read.
278. arrange(-error_rate)
279. ```
281. Are differences significant ? We use an anova on error rate by stimulus here, because we would have dependency in the data otherwise.
283. Note that theoretically we have similar dependency by subject (children saw an unbalanced amount of stimuli per category), but we assume that this will not be a problem because of a large sample size.
285. ```{r errorRatesPerCategory}
286. # Not happy with this since participants have not seen the same amount of stimuli in each category
287. # You can't group by subject either because you then have a dependency on stimulus
288. d %>%
289. # For each target, compute error rate
290. group_by(stimulus, category) %>% summarise(error_rate = mean(!is_correct), .groups = "drop") %>%
291. # Run an anova by category
292. rstatix::anova_test(formula = error_rate ~ category) %>%
293. # Print the results in APA-friendly format
294. rstatix_anova_to_APA()
295. ```
297. # Variations in difficulties - IRT
299. ## IRT difficulty
301. ```{r irtDifficulty}
302. irt_2PL_model.plci = d %>%
303. # Change correctness to numeric value, which is required for IRT analysis
304. mutate(is_correct = as.numeric(is_correct)) %>%
305. # Only keep the relevant stuff for IRT: subject, target, and outcome of the question
306. select(subj, stimulus, is_correct) %>%
307. # Format everything with one column per target to fit IRT format requirement.
308. pivot_wider(
309. names_from = stimulus,
310. values_from = is_correct
311. ) %>%
312. # Drop the subject model, as IRT does not care about identifiers.
313. select(-subj) %>%
314. # Fit IRT model
315. mirt(itemtype = "2PL", technical = list(NCYCLES = 10000)) %>%
316. # Get 95% CI - this may be long
317. PLCI.mirt %>%
318. # The reported d is easiness, so change it since we want difficulty in order to match error rates
319. # https://groups.google.com/g/mirt-package/c/RwGmO_fGi2c
320. mutate(parnam = if_else(parnam == "d", "difficulty", parnam),
321. value = if_else(parnam == "difficulty", -value, value),
322. lower_2.5 = if_else(parnam == "difficulty", -lower_2.5, lower_2.5),
323. upper_97.5 = if_else(parnam == "difficulty", -upper_97.5, upper_97.5))
324. ```
326. What is the link between difficulty and error rate
327. ```{r irtDifficultyAndErrorRateLM}
328. irt_2PL_model.plci %>%
329. # Only care about difficulty, not discrimination
330. filter(parnam == "difficulty") %>%
331. # Only look at the target and the result
332. select(Item, value) %>% dplyr::rename(stimulus = Item, difficulty = value) %>%
333. # Add the average error rate for each target
334. inner_join(d %>% group_by(stimulus) %>% summarise(error_rate = mean(!is_correct), .groups = "drop"), by = c("stimulus")) %>%
335. # Compute the linear regression.
336. lm(formula = error_rate ~ difficulty) %>%
337. # Print the result according to APA norms
338. summary %>% regression_summary_to_APA
339. ```
341. ### Error rate & IRT figure
343. ### Figure 2
344. ```{r errorRateAndDifficultyByStimulus, fig.width = 8.27, fig.height = 11.69/2}
345. # Make the plot for error rates
346. error_rates.pl = error_rates.d %>%
347. # Plot error rate by stimulus (ordered by error rate) and colored by category
348. ggplot(aes(x = error_rate, y = reorder(stimulus, error_rate), fill = category,
349. xmin = error_rate, xmax = error_rate + error_rate.half_ci)) +
350. # Plot the error rate as bars...
351. geom_bar(stat = "identity") +
352. # ... with standard error
353. geom_errorbar(width = 0) +
354. # Ensure that similar ylims and expansion are used on the two plots!
355. coord_cartesian(xlim = c(0.6, 1.001), ylim = c(0.5, d$stimulus %>% unique %>% length + 0.5), expand = F) +
356. # Plot error rate as a percentage
357. scale_x_continuous(labels = scales::percent_format()) +
358. # Define axis labels
359. labs(y = "Target", x = "Error rate", fill = "Category") +
360. # Make a clean plot
361. theme(axis.title.y = element_blank(), panel.background = element_blank(), panel.grid = element_blank(),
362. # Add some lines to delimit everything
363. axis.line = element_line(colour = "black"), panel.grid.major.x = element_line(colour = "gray", linetype = 1),
364. # Use the same margin on the two plots!
365. plot.margin = margin(0, 1,0,0.1, unit = "cm"),
366. # No use for a legend in this plot.
367. legend.position = "none")
369. # Make the difficulty plot
370. difficulty.pl = irt_2PL_model.plci %>%
371. # We only care about difficulty
372. filter(parnam == "difficulty") %>%
373. # Keep relevant parameters (target, value and 95% CI)
374. select(Item, value, lower_2.5, upper_97.5) %>%
375. # Add error rates so that target order matches the error rate plot.
376. inner_join(d %>% group_by(stimulus, category) %>% summarise(error_rate = mean(!is_correct), .groups = "drop"), by = c("Item" = "stimulus")) %>%
377. # Plot difficulty by item, with items ordered by error rates just like the others.
378. ggplot(aes(y = reorder(Item, error_rate), x = value, fill = category,
379. xmin = lower_2.5, xmax = upper_97.5)) +
380. # Add lines at hard and very hard boundaries.
381. geom_vline(xintercept = c(0.5, 2), colour = "red", linetype = 2) +
382. # Add errorbars
383. geom_errorbar(width = 0.5) +
384. # Add points, with a triangle shape like usual IRT plots.
385. geom_point(size = 5, shape = 22) +
386. # Make breaks for difficulty every 0.5
387. scale_x_continuous(breaks = c(0.5, seq(2,8, 2)), labels = function(x){signif(x, 0)}) +
388. # Ensure that similar ylims and expansion are used on the two plots!
389. coord_cartesian(xlim = c(0, 8), ylim = c(0.5, d$stimulus %>% unique %>% length + 0.5), expand = F) +
390. # Add axis labels
391. labs(x = "Difficulty", y = "Target", fill = "Category") +
392. # Clean the plot
393. theme(axis.title.y = element_blank(), panel.background = element_blank(), panel.grid = element_blank(),
394. # Add some lines to help readability
395. axis.line = element_line(colour = "black"), panel.grid.major.x = element_line(colour = "gray"),
396. # Use some margin to make axes legible - use the same margin on the two plots!
397. plot.margin = margin(t = 0, r = 1, b = 0, l = 0.1, unit = "cm"),
398. # Place the legend as inset
399. legend.position = c(0.7, 0.15), legend.background = element_rect(colour = "black"),
400. # Remove legend grey background
401. legend.key = element_blank(),
402. # Just set usual text (this is only meant for tweaking)
403. legend.text = element_text(hjust = 0),
404. # Black frame
405. legend.margin = margin(t = 0.3, r = 1.1, b = 0.3, l = 0.3, "cm"),
406. # Space between entries
407. legend.key.height = unit(1, "cm")
408. ) +
409. # Specifiy the size of fill indications in the legend.
410. guides(fill = guide_legend(override.aes = list(size = 7)))
412. # Create a two-plots figure and save it
413. ggarrange(error_rates.pl, difficulty.pl, ncol = 2) %>%
414. save_on_disk("Figure_2", path = "./figures/main/", a4_fraction = c(1, 1/2))
415. ```
418. # Non-random errors
419. ## Precision
420. ```{r precision}
421. # Create a function that finds an print the proportion of responses fallling within a given threshold of the correct answer.
422. (function(threshold){
423. d %>%
424. # Drop non-responses and wrong responses????
425. filter(!actual_response %>% is.na, !is_correct) %>%
426. # Only keep relevant data
427. select(stimulus, response, actual_response) %>%
428. # Compute the distance to the response: if it is below the threshold, call it a close response.
429. mutate(response_distance = abs(actual_response - response),
430. is_close_response = response_distance <= threshold) %>%
431. # For each target, compute the frequency of close response.
432. group_by(stimulus) %>% summarise(close_response_rate = mean(is_close_response),.groups = "drop") %>%
433. # Compute the average rate across targets.
434. summarise(close_response_rate.se = se(close_response_rate), close_response_rate = mean(close_response_rate)) %>%
435. # Print the result
436. (function(d){cat("The proportion of responses falling within ", threshold, " of the correct answer is ",
437. d$close_response_rate[[1]] %>% (scales::percent_format(accuracy=0.1)), "%",
438. " (se = ", d$close_response_rate.se[[1]] %>% (scales::percent_format(accuracy = 0.01)), "%).",
439. sep = "")
440. })
441. }) (0.25)
442. ```
444. ## Imbalanced distributions
445. Are there stimulus where the response distribution is not uniform?
446. ```{r distribution of errors per stimulus}
447. d %>%
448. # Count possible responses, it happens that each response that was never given is incorrect
449. count_possible_responses(is_correct) %>% mutate(is_correct = is_correct %>% replace_na(F)) %>%
450. # Only consider erroneous responses: drop non-responses and correct responses
451. filter(!actual_response %>% is.na, !is_correct) %>%
452. # Only care about the number of responses for each position and target -- the target itself is not relevant
453. select(stimulus, n) %>%
454. # Perform a chi-square test on each target and put the results in adequate columns.
455. group_by(stimulus) %>% summarise(test = chisq.test(n) %>% (broom::tidy), .groups = "drop") %>% unnest %>%
456. # Rename columns (and drop useless method report)
457. dplyr::rename(X2 = statistic, df = parameter, p = p.value) %>% select(-method) %>%
458. # FDR-correction for p-values, and put them in nice APA format
459. mutate(p = p %>% p.adjust("fdr") %>% p_value_to_APA()) %>%
460. # Order target from most significative test to the least significative one.
461. arrange(X2) %>%
462. # Round as per APA format
463. mutate(X2 = X2 %>% flexible_round(2)) %>%
464. # Report results in a human-friendly format.
465. mutate(stats = str_c("\\Chi2(", df, ") = ", X2, ", ", p), .keep = "unused")
466. ```
468. Can we explain the imbalance in responses based on integers?
469. ```{r error distribution without integer}
470. d %>%
471. # Count possible responses, it happens that each response that was never given is incorrect
472. count_possible_responses(is_correct) %>% mutate(is_correct = is_correct %>% replace_na(F)) %>%
473. # Identify integer responses
474. mutate(is_integer = actual_response %>% round == actual_response) %>%
475. # Only consider erroneous responses: drop non-responses, correct responses, and integer responses
476. filter(!actual_response %>% is.na, !is_correct, !is_integer) %>%
477. # For each target, perform a chi-square test and put the results in adequate columns.
478. group_by(stimulus) %>% summarise(test = chisq.test(n) %>% (broom::tidy), .groups = "drop") %>% unnest %>%
479. # Rename columns (and drop useless method report)
480. dplyr::rename(X2 = statistic, df = parameter, p = p.value) %>% select(-method) %>%
481. # FDR-correction for p-values, and put them in nice APA format
482. mutate(p = p %>% p.adjust("fdr") %>% p_value_to_APA()) %>%
483. # Order target from most significative test to the least significative one.
484. arrange(X2) %>%
485. # Round as per APA format
486. mutate(X2 = X2 %>% flexible_round(2)) %>%
487. # Report results in a human-friendly format.
488. mutate(stats = str_c("\\Chi2(", df, ") = ", X2, ", ", p), .keep = "unused")
489. ```
491. All the targets display a significant imbalance of error positions. The df varies between 43 and 44 across tests because sometimes the correct response is an integer (so we remove only 6 response positions for the test), sometimes it isn't (so we remove 7 response positions).
493. To what extent are error positions falling on a restricted set of responses?
494. ```{r error position stacking}
495. d %>%
496. # Drop non-responses and correct responses
497. filter(!actual_response %>% is.na, !is_correct) %>%
498. # Count how many times each response was given for each target.
499. count(stimulus, actual_response) %>%
500. # Perform this analysis per target
501. group_by(stimulus) %>%
502. # Turn the count into a proportion per target.
503. mutate(proportion = n / sum(n)) %>%
504. # Arrange responses by decreasing proportion: most given responses are on top...
505. arrange(stimulus, -proportion) %>%
506. # And now make a cumulated sum of the proportions: each sum covers the most given responses since we ordered them above!
507. mutate(cumsummed_proportion = cumsum(proportion),
508. # Also count the number of positions included in the cumsum.
509. actual_response_index = 1:n()) %>%
510. # Only keep the cumsums after which the 50% of responses threshold was hit.
511. filter(cumsummed_proportion > .50) %>%
512. # And keep the one with the least responses included: this is how many positions are required to reach the threshold (again because we ordered them by decreasing order of proportion of responses before)
513. filter(actual_response_index == min(actual_response_index)) %>%
514. # Ungroup, we now consider results across targets again.
515. ungroup %>%
516. # Only keep relevant info: targets, and the number of positions. required to reach 50% of responses.
517. select(stimulus, actual_response_index) %>%
518. # Compute the average number of responses required, and print the result in a nice
519. summarise(mean_index = mean(actual_response_index)) %>%
520. (function(d){cat("On average, half the responses fall on", flexible_round(d$mean_index, 2), "positions.")})
521. ```
523. # Inversion
524. ## Four targets
526. First, let's consider four remarkable targets: 1/2 (a basic fraction, if not the most basic one); 3/6 (an equivalent to 1/2); 2/1 (the inverse of 1/2); and 6/3 (the inverse of 3/6)
528. ```{r four targets example, fig.width = 8.27, fig.height = 11.69/2}
529. (d %>%
530. # Only keep our targets of interest
531. filter(stimulus %in% c("1 / 2", "3 / 6", "2 / 1", "6 / 3"),
532. # Drop non-responses.
533. !actual_response %>% is.na) %>%
534. # Count the number of responses per target (also keep some info about the target like components and whether the answer is correct or not)
535. group_by(stimulus, numerator, denominator, response) %>% count(actual_response, is_correct) %>%
536. # Turn the number of responses into a proportion per target, and count the total number of responses per target.
537. group_by(stimulus) %>% mutate(proportion = n / sum(n), sample_size_label = str_c("n = ", sum(n))) %>%
538. # Add a some infos about the target as label for later plotting.
539. # We only add a single label per target to avoid cluttering the plot, on the correct response for convenience
540. # (Note that here we have at least one correct responses per target, but this line may fail if there is none.)
541. mutate(sample_size_label = if_else(is_correct, sample_size_label, NA),
542. numerator_label = if_else(is_correct, numerator, NA),
543. denominator_label = if_else(is_correct, denominator, NA),
544. response_label_position = if_else(is_correct, response, NA)) %>%
545. ungroup %>%
546. # Match responses with predictions from our patterns to mark the responses on the plot
547. left_join(
548. # If several patterns yield the same prediction, only keep the first one as it is just illustrative
549. fractions_predictions.d %>% group_by(stimulus, prediction) %>% filter(1:n() == 1) %>%
550. # Separate the 'normative' pattern into correct (without inversion3) and inverted
551. mutate(pattern = if_else(pattern == "Normative", if_else(inverted, "Inverted", "Correct"), pattern)),
552. # Join the predictions with actual responses
553. by = c("stimulus", "actual_response" = "prediction")) %>%
554. # Order the stimuli neatly for the facetted display
555. mutate(stimulus = stimulus %>% factor(c("1 / 2", "3 / 6", "2 / 1", "6 / 3"))) %>%
556. # Order the patterns neatly for the legend
557. mutate(pattern = pattern %>% factor(levels =
558. patterns.colours[patterns.colours %>% names != "Normative"] %>% names %>%
559. # Add blanks to correctly shape the columns - go in decreasing order to avoid shifting indices
560. append(c("BLANK3", "BLANK2"), after = 5) %>% append("BLANK1", after = 3)
561. )) %>%
562. # Plot the proportion of responses per position, filling the bars matching specific patterns
563. ggplot(aes(x = actual_response, y = proportion, fill = pattern)) +
564. # Add bars to indicate the proportion of resopnses.
565. geom_bar(stat = "identity") +
566. # Divide by size in mm by .pt to turn it into point, and then multiply by PT_TO_PIXEL to correct SVG saving sizes in pixels.
567. # Write the sample size in the top right corner
568. geom_text(aes(label = sample_size_label), x = 4.5, y = .36, size = 15 * PT_TO_PIXEL / .pt, family = FAMILY) +
569. # Write the target fraction above the response (the correct one, since all others have NAs)
570. geom_text(aes(label = numerator_label, x = response), y = .38, size = 18 * PT_TO_PIXEL / .pt, family = FAMILY) +
571. geom_text(aes(label = denominator_label, x = response), y = .34, size = 18 * PT_TO_PIXEL / .pt, family = FAMILY) +
572. # Write the vinculum of the fraction
573. geom_segment(aes(x = response_label_position - 0.1, xend = response_label_position + 0.1),
574. y = .36, yend = .36, linewidth = 0.5) +
575. # Indicate the correct response with an arrow poiting to it.
576. ## Body
577. geom_segment(aes(x = response_label_position, xend = response_label_position,
578. y = proportion + 0.02),
579. yend = .31) +
580. ## Head
581. geom_segment(aes(x = response_label_position, xend = response_label_position - 0.1,
582. y = proportion + 0.02, yend = proportion + 0.04)) +
583. geom_segment(aes(x = response_label_position, xend = response_label_position + 0.1,
584. y = proportion + 0.02, yend = proportion + 0.04)) +
585. # Create one subplot for each of the four targets, in a 2x2 display
586. facet_wrap(~stimulus, nrow = 2, scales = "free_x") +
587. # Indicate axes labels
588. labs(x = "Responded position", y = "% of responses", fill = "Response type") +
589. # Set the limits of the y axis (determined by hand) and display % of responses as actual %s.
590. scale_y_continuous(limits = c(0, 0.4), labels = scales::percent_format()) +
591. # Clean the panel to make it easier to read.
592. theme(panel.background = element_blank(), panel.grid = element_blank(), panel.spacing = unit(0.5, "cm"),
593. # Add some margins to space everything.
594. plot.margin = margin(t = 0.6, r = 0.5, b = 0.5, l = 0.5, "cm"),
595. # Remove the strip from facet plot
596. strip.background = element_blank(), strip.text = element_blank(),
597. # Make the titles
598. #axis.title = element_text(),
599. axis.title.x = element_text(hjust = 0.14, vjust = 0),
600. axis.line = element_line(colour = "black"), panel.grid.major.y = element_line(colour = "gray", linetype = 1),
601. legend.position = "bottom",
602. # Remove grey background from keys
603. legend.key = element_blank(), legend.title = element_text(vjust = 0.75)) +
604. # legend.position = "none") +
605. # Add a horizontal y axis legend - small hacky stuff with subtitle
606. theme(axis.title.y = element_blank(),
607. plot.title = element_blank(), plot.subtitle = element_text(hjust = -0.08, vjust = 4)) +
608. ggtitle("", subtitle = "% of responses") +
609. scale_x_decimal_number_line(expand = F) +
610. guides(fill = guide_legend(ncol = 3), override.aes = aes(width = 1)) +
611. # This is a very hacky way to get the shape of the legend right
612. scale_fill_manual(values = c(patterns.colours, "BLANK1" = "white", "BLANK2" = "white", "BLANK3" = "white"),
613. labels = Vectorize(function(p){
614. if(p == "BLANK1") return("to the line")
615. else if(!p %in% names(pretty_pattern_names)) return("")
616. else pretty_pattern_names[[p]] %>% str_remove(" to the line")}),
617. drop = FALSE
618. )
619. ) %>%
620. save_on_disk("Figure_3", path = "./figures/main/SVG/", device = "svg", a4_fraction = c(1, .60), indicate_format = F)
621. ```
623. ## Pairwise confusion
625. Compute the similarity between distributions for pairs of inverted fractions.
627. ```{r fig.width= 8.27/2, fig.height = 11.69/2}
628. # Compute the list of pairs of inverse fractions, ordered to make plots using it easier to read.
629. ordered_inverse_targets_pairs = d %>%
630. # Only consider the targets and their components
631. select(stimulus, numerator, denominator) %>% unique %>%
632. # Sort the components so that pairs of inverse fractions have the same components in the same order
633. mutate(smaller_component = pmin(numerator, denominator),
634. greater_component = pmax(numerator, denominator)) %>%
635. # Only keep fractions where you can build a pair with same components (i.e., inverted fractions)
636. group_by(smaller_component, greater_component) %>% filter(n() == 2) %>%
637. # Reorder everything pairwise to have an order for future plots
638. arrange(greater_component, smaller_component, numerator) %$%
639. # Keep the list of targets
640. stimulus
642. # Compute the correlation between distributions -- to show that it is greater for inverted fractions
643. inverted_pairs_distributions_correlations = d %>%
644. # Count the number of responses per position for each target, including possible positions without any response, since we care about distributions.
645. count_possible_responses(numerator, denominator) %>%
646. # Only keep stimuli for which an inverse fraction was presented
647. filter(stimulus %in% ordered_inverse_targets_pairs) %>%
648. # Keep relevant infos: target, response position, number of responses.
649. select(stimulus, actual_response, n) %>%
650. # For each position, count the number of responses that fell on it per target
651. group_by(stimulus) %>% pivot_wider(names_from = stimulus, values_from = n) %>%
652. # Now drop the exact positions, since only care about the distributions
653. select(-actual_response) %>%
654. # Compute the correlation between each line of the matrix...
655. cor %>%
656. # ...and report it with one line per pair.
657. melt %>%
658. # Rename the columns and make the targets into factor for a consistent order.
659. mutate(target1 = X1 %>% factor(levels = ordered_inverse_targets_pairs),
660. target2 = X2 %>% factor(levels = ordered_inverse_targets_pairs),
661. pearson_r = value,
662. .keep = "unused")
664. # Make a plot
665. (inverted_pairs_distributions_correlations %>%
666. # Identify target 1 as x_label
667. dplyr::rename(x_label = target1) %>%
668. # Reinject numerators and denominators to make computations with them on x label.
669. inner_join(d %>% select(stimulus, numerator, denominator) %>% unique,
670. by = c("x_label" = "stimulus")) %>%
671. # Place fractions in a pair on alternating lines, by adding a line return for those > 1
672. mutate(x_label = str_c(if_else(numerator > denominator, "\n", ""), x_label %>% str_remove_all(" "))) %>%
673. # Plot the correlation between the first target (x-axis) and second target (y-axis)
674. ggplot(aes(x = reorder(x_label, target1 %>% as.numeric), y = target2, fill = pearson_r)) +
675. # Tile with pearson correlation between distributions
676. geom_tile(color = NA) +
677. # Ensure that the matrix is a square
678. coord_fixed() +
679. # Set a nicer scale for filling the tiles
680. scale_fill_viridis_c(option = "magma") +
681. # Set labels
682. labs(x = "Target 1", y = "Target 2", fill = "Pearson's r") +
683. # Clean the plot
684. theme(panel.background = element_blank(), panel.grid = element_blank(),
685. # Drop the axes ticks and title
686. axis.ticks = element_blank(), axis.title = element_blank(),
687. # Set the legend
688. legend.position = "top", legend.key.width = unit(1.5, "cm"), legend.title = element_text(vjust = 1))
689. ) %>%
690. # Save
691. save_on_disk("Figure_4", path = "./figures/main/SVG/", a4_fraction = c(1/2, 1/2), indicate_format = F)
692. ```
694. Is the correlation stronger within pairs (big squares on the diagonal) than across?
695. ```{r correlationForInvertedPairsSummary}
696. # Compute the correlation for each non-ordered pair
697. non_ordered_pairs_distributions_correlations = inverted_pairs_distributions_correlations %>%
698. # Add back the numerator and denominator of target 1
699. inner_join(d %>% select(stimulus, numerator, denominator) %>% unique %>% dplyr::rename(target1 = stimulus, numerator1 = numerator, denominator1 = denominator),
700. by = c("target1")) %>%
701. # Add back the numerator and denominator of target 1
702. inner_join(d %>% select(stimulus, numerator, denominator) %>% unique %>% dplyr::rename(target2 = stimulus, numerator2 = numerator, denominator2 = denominator),
703. by = c("target2")) %>%
704. # We only want to consider non-ordered pairs, so only keep those following the lexicographic order
705. filter(numerator1 < numerator2
706. | numerator1 == numerator2 & denominator1 < denominator2) %>%
707. # Flag inverted pairs
708. mutate(is_inverted_pair = numerator1 == denominator2 & denominator1 == numerator2, .keep = "unused")
710. non_ordered_pairs_distributions_correlations %>%
711. # Compute the average correlation within and across inverted pairs
712. group_by(is_inverted_pair) %>% summarise(
713. pearson_r.se = se(pearson_r), pearson_r = mean(pearson_r), .groups = "drop") %>%
714. # Report the results
715. select(is_inverted_pair, pearson_r, pearson_r.se)
716. ```
718. ```{r correlationForInvertedPairsStats}
719. # Look how the inverted pairs rank in terms of correlations compared to others
720. non_ordered_pairs_distributions_correlations %>% arrange(-pearson_r)
722. ```
724. There is only one non-inverted pair for which the correlation is better than a single inverted pair
726. Since we have 6 relevant pairs (those of inverted pairs) and 66 pairs in total, there are 6C66 possible orderings (~90,858,768)
728. The only best ordering would have been one where the 6 inverted pairs had a better correlation than all others, so we have the 2nd to best ordering(with 6 targets, that is inverse pairs): this means that in exact permutation test, p = 2/(6C66) = 2/~90,858,768 < .001
730. # Global error patterns
732. Now evaluate the patterns across all targets
733. ```{r total prediction}
734. d %>%
735. # Remove non-responses
736. filter(!actual_response %>% is.na, !is_correct) %>%
737. # Count reponses per stimulus and position
738. count(stimulus, actual_response) %>%
739. # Convert the count into a proportion to allow for comparison across stimulus
740. group_by(stimulus) %>% mutate(proportion = n/sum(n)) %>% ungroup %>%
741. # Only keep stimulus-response pairs that are also in the predictions df.
742. # Drop the patterns themselves as we don't care about them
743. inner_join(fractions_predictions.d %>% select(stimulus, prediction) %>% unique,
744. by = c("stimulus", "actual_response" = "prediction")) %>%
745. # Since we have disjoint entries in the set, we can compute the proportion explained per stimulus as the sum.
746. group_by(stimulus) %>% summarise(proportion_explained = sum(proportion), .groups = "drop") %>%
747. # Now compute the average proportion across all stimuli
748. mutate(mean_proportion_explained = mean(proportion_explained)) %>%
749. # And order them from the least explained.
750. arrange(proportion_explained)
751. ```
753. Compute the proportion of non-ambiguously explained responses
754. ```{r unambiguously predicted}
755. d %>%
756. # Ignore non-responses
757. filter(!actual_response %>% is.na) %>%
758. # Count each response on each stimulus
759. count_possible_responses() %>%
760. # Turn count into proportion to allow for comparison across stimuli
761. group_by(stimulus) %>% mutate(proportion = n/sum(n)) %>% ungroup %>%
762. # Match responses to pattern prediction, but only when not ambiguous!
763. inner_join(fractions_predictions.d %>% filter(!ambiguous),
764. by = c("stimulus", "actual_response" = "prediction")) %>%
765. # Compute the proportion of responses explained by a pattern across stimuli
766. group_by(pattern, inverted) %>% summarise(proportion.se = se(proportion), proportion = mean(proportion), .groups = "drop") %>%
767. # Reorder the columns...
768. select(pattern, inverted, proportion, proportion.se) %>%
769. # ...and round to ease readability
770. mutate(proportion = scales::label_percent(accuracy = 0.1)(proportion),
771. proportion.se = scales::label_percent(accuracy = 0.1)(proportion.se))
772. ```
774. ## Overrepresentation analyses across all predictions
776. Use all targets
777. ```{r global overrepresentation with all targets}
778. d %>%
779. # Remove non-responses
780. filter(!actual_response %>% is.na) %>%
781. # Count the number of responses on all positions
782. count_possible_responses(stimulus, actual_response) %>%
783. # Just order by stim and response for easier readability
784. arrange(stimulus, actual_response) %>%
785. # Turn count into proportion
786. group_by(stimulus) %>% mutate(proportion = n / sum(n)) %>%
787. # Compute the median proportion of responses on each position
788. group_by(actual_response) %>% mutate(median_proportion_on_position = median(proportion)) %>%
789. # Match responses to predictions -- only keep non ambiguous cases to avoid inflating our results from other patterns
790. inner_join(fractions_predictions.d %>% filter(!ambiguous),
791. by = c("stimulus", "actual_response" = "prediction")) %>%
792. # For each pattern (with/without inversion) do a paired wilcoxon test on whether the proportion is consistently greater than the median. Compute also the mean proportion and number of targets as supplementary info
793. group_by(pattern, inverted) %>%
794. summarise(number_of_unambiguous_targets = n(),
795. test = wilcox.test(proportion, median_proportion_on_position, paired = T, alternative = "greater") %>% (broom::tidy),
796. mean_proportion = mean(proportion),
797. .groups = "drop") %>%
798. # Unpack test data
799. unnest(cols = test) %>%
800. # Correct for multiple testing
801. mutate(p.value = p.value %>% p.adjust("fdr")) %>%
802. # Order by pattern to ease readability2
803. mutate(pattern = factor(pattern, levels = patterns.colours %>% names)) %>% arrange(pattern) %>%
804. # Clean the report -- merge pattern name and inversion
805. mutate(pattern = str_c(pattern, if_else(inverted, " (inverted)", "")), .keep = "unused") %>%
806. # Turn statistics into APA format
807. mutate(p.value = p.value %>% p_value_to_APA, mean_proportion = scales::percent_format(accuracy = 0.1)(mean_proportion)) %>%
808. # Reorder and only keep relevant columns
809. select(pattern, p.value, mean_proportion, number_of_unambiguous_targets, statistic)
810. ```
812. Exclude the four examples
813. ```{r global overrepresentation with naive targets}
814. d %>%
815. # Remove non-responses
816. filter(!actual_response %>% is.na) %>%
817. # Ignore the worked out examples
818. filter(!stimulus %in% c("1 / 2", "2 / 1", "3 / 6", "6 / 3")) %>%
819. # Count the number of responses on all positions
820. count_possible_responses(stimulus, actual_response) %>%
821. # Just order by stim and response for easier readability
822. arrange(stimulus, actual_response) %>%
823. # Turn count into proportion
824. group_by(stimulus) %>% mutate(proportion = n / sum(n)) %>%
825. # Compute the median proportion of responses on each position
826. group_by(actual_response) %>% mutate(median_proportion_on_position = median(proportion)) %>%
827. # Match responses to predictions -- only keep non ambiguous cases to avoid inflating our results from other patterns
828. inner_join(fractions_predictions.d %>% filter(!ambiguous),
829. by = c("stimulus", "actual_response" = "prediction")) %>%
830. # For each pattern (with/without inversion) do a paired wilcoxon test on whether the proportion is consistently greater than the median. Compute also the mean proportion and number of targets as supplementary info
831. group_by(pattern, inverted) %>%
832. summarise(number_of_unambiguous_targets = n(),
833. test = wilcox.test(proportion, median_proportion_on_position, paired = T, alternative = "greater") %>% (broom::tidy),
834. mean_proportion = mean(proportion),
835. .groups = "drop") %>%
836. # Unpack test data
837. unnest(cols = test) %>%
838. # Correct for multiple testing
839. mutate(p.value = p.value %>% p.adjust("fdr")) %>%
840. # Order by pattern to ease readability2
841. mutate(pattern = factor(pattern, levels = patterns.colours %>% names)) %>% arrange(pattern) %>%
842. # Clean the report -- merge pattern name and inversion
843. mutate(pattern = str_c(pattern, if_else(inverted, " (inverted)", "")), .keep = "unused") %>%
844. # Turn statistics into APA format
845. mutate(p.value = p.value %>% p_value_to_APA, mean_proportion = scales::percent_format(accuracy = 0.1)(mean_proportion)) %>%
846. # Reorder and only keep relevant columns
847. select(pattern, p.value, mean_proportion, number_of_unambiguous_targets, statistic)
848. ```
850. Restrict some patterns
851. ```{r global overrepresentation with restricted targets}
852. d %>%
853. # Remove non-responses
854. filter(!actual_response %>% is.na) %>%
855. # Count the number of responses on all positions
856. count_possible_responses(stimulus, actual_response) %>%
857. # Just order by stim and response for easier readability
858. arrange(stimulus, actual_response) %>%
859. # Turn count into proportion
860. group_by(stimulus) %>% mutate(proportion = n / sum(n)) %>%
861. # Compute the median proportion of responses on each position
862. group_by(actual_response) %>% mutate(median_proportion_on_position = median(proportion)) %>%
863. # Match responses to predictions -- only keep non ambiguous cases to avoid inflating our results from other patterns
864. inner_join(fractions_predictions.d %>%
865. filter(!ambiguous,
866. !pattern %in% c("Relative", "Mixed")
867. | pattern == "Relative" & prediction == 2.5
868. | pattern == "Mixed" & stimulus %>% str_count(" / 2") == 1
869. ),
870. by = c("stimulus", "actual_response" = "prediction")) %>%
871. # For each pattern (with/without inversion) do a paired wilcoxon test on whether the proportion is consistently greater than the median. Compute also the mean proportion and number of targets as supplementary info
872. group_by(pattern, inverted) %>%
873. summarise(number_of_unambiguous_targets = n(),
874. test = wilcox.test(proportion, median_proportion_on_position, paired = T, alternative = "greater") %>% (broom::tidy),
875. mean_proportion = mean(proportion),
876. .groups = "drop") %>%
877. # Unpack test data
878. unnest(cols = test) %>%
879. # Correct for multiple testing
880. mutate(p.value = p.value %>% p.adjust("fdr")) %>%
881. # Order by pattern to ease readability2
882. mutate(pattern = factor(pattern, levels = patterns.colours %>% names)) %>% arrange(pattern) %>%
883. # Clean the report -- merge pattern name and inversion
884. mutate(pattern = str_c(pattern, if_else(inverted, " (inverted)", "")), .keep = "unused") %>%
885. # Turn statistics into APA format
886. mutate(p.value = p.value %>% p_value_to_APA, mean_proportion = scales::percent_format(accuracy = 0.1)(mean_proportion)) %>%
887. # Reorder and only keep relevant columns
888. select(pattern, p.value, mean_proportion, number_of_unambiguous_targets, statistic)
889. ```
891. ## Overrepresentation of individual predictions
893. ```{r overrepresentation of individual predictions pvalues}
894. overrepresentation.d = d %>%
895. # Drop non-responses
896. filter(!actual_response %>% is.na) %>%
897. # Count answers
898. count_possible_responses(stimulus, actual_response) %>%
899. # Plug back some information
900. # - Category, numerator, denominator for ordering
901. inner_join(d %>% select(stimulus, category, numerator, denominator, response) %>% unique,
902. by = c("stimulus")) %>%
903. # Turn count into proportion
904. group_by(stimulus) %>% mutate(total_n = sum(n), proportion = n / sum(n), ) %>%
905. # Compute the median proportion of responses on each position
906. group_by(actual_response) %>% mutate(median_proportion_on_position = median(proportion)) %>% ungroup %>%
907. # Match the responses with patterns
908. # Do not care about ambiguity since we are concerned about individual response positions
909. inner_join(fractions_predictions.d %>% select(stimulus, prediction, pattern, inverted),
910. by = c("stimulus", "actual_response" = "prediction")) %>%
911. # Order by pattern for readability
912. arrange(pattern, inverted) %>%
913. # Do a one-proportion t-test, for each entry, against the median proportion
914. mutate(test = Vectorize(function(x, n, p){prop.test(x, n, p = p, alternative = "greater") %>% (broom::tidy) %>% list})(n, total_n, median_proportion_on_position)) %>% unnest(cols = c(test)) %>%
915. # Only keep some relevant values and reorder them
916. select(pattern, inverted, stimulus, actual_response, proportion, statistic, parameter, p.value, median_proportion_on_position, n, total_n) %>%
917. # Correct p-values
918. mutate(p.value = p.value %>% p.adjust("fdr"))
920. overrepresentation.d %>%
921. # Format statistics
922. mutate(p_value = p.value %>% p_value_to_APA,
923. proportion = proportion %>% (scales::percent_format(accuracy = 1)),
924. z = statistic %>% flexible_round(digits = 2),
925. .keep = "unused") %>%
926. # Report statistical test as per APA standards
927. mutate(statistical_test = str_c("z(", parameter, ", ", total_n, ") = ", z), .keep = "unused") %>%
928. # Merge pattern name and inversion
929. mutate(pattern = str_c(pattern, if_else(inverted, "_inverted", "")), .keep = "unused") %>%
930. # Reorder relevant columns
931. select(pattern, stimulus, proportion, statistical_test, p_value)
933. ```
935. #### Figure 5
937. Plot the distribution of responses as well as some focus on specific patterns
938. ```{r global response distribution, fig.width = 8.27}
939. full_response_patterns.pl = d %>%
940. # Drop non-responses
941. filter(!actual_response %>% is.na) %>%
942. # Count the number of responses per target and position, also keep some useful data
943. count(stimulus, category, numerator, denominator, is_correct, actual_response) %>%
944. # Turn count into proportion per target
945. group_by(stimulus) %>% mutate(proportion = n / sum(n)) %>% ungroup() %>%
946. # Order the categories
947. mutate(category = category %>% factor(c("n / 1", "Tenths", "Integer result", "Equivalent to 1/2"))) %>%
948. # Plot the distribution of response, with targets ordered
949. # - by category
950. # - then by denominator,
951. # - then by numerator
952. ggplot(aes(x = actual_response,
953. y = reorder(stimulus, 10000 * as.numeric(category) + 100 * denominator + numerator))) +
954. # Fill each tile with the proportion of correct responses
955. geom_tile(aes(fill = proportion)) +
956. # Add an outline to correct responses
957. geom_tile(aes(colour = is_correct, linewidth = is_correct), alpha = 0) +
958. # Add a number line for the x-axis
959. scale_x_decimal_number_line(expand = F, major.tick.size = 7) +
960. # Make the y-axis neater
961. scale_y_discrete(guide = "prism_offset") +
962. # Color the correct responses with the color of the 'Correct' pattern for consistency
963. scale_colour_manual(values = c("TRUE" = patterns.colours[["Normative"]]), breaks = c("TRUE"), labels = "Correct response") +
964. # Set the size of the outline of the correct responses
965. scale_linewidth_manual(values = c("TRUE" = 1.25), breaks = c()) +
966. # Set the fill scale for the proportion of responses as a grey-scale gradient
967. scale_fill_gradient(low = "white", high = "black", labels = scales::percent_format(), limits = c(0, .4), breaks = seq(0, .3, .1)) +
968. # Set the linewith for the 'correct' response legend.
969. guides(colour = guide_legend(override.aes = list(linewidth = 1.25), title = NULL)) +
970. # Set nicer labels for axes
971. labs(x = "Response position", y = "Target fraction", fill = "% of responses", colour = "Correct responses") +
972. # Set y label in a horizontal position
973. horizontal_ylab("Target\nfraction", vjust = 1.15) +
974. # Set the title to be a tiny bit greater.
975. theme_text(size = 14, title.size = 16, scaling_factor = PT_TO_PIXEL) +
976. # Clean the plot
977. theme(panel.background = element_blank(), panel.grid = element_blank(),
978. panel.grid.major.y = element_blank(),
979. # Set a line on the y-axis
980. axis.line = element_line(colour = "black"),
981. # Make the legend prettier
982. legend.key = element_blank(), legend.position = "bottom",
983. legend.text = element_text(vjust = 0, hjust = 0.5),
984. legend.box.just = "right",
985. legend.title = element_text(vjust = 0.9), legend.key.width = unit(1.6, "cm"),
986. legend.spacing.y = unit(0, "mm"),
987. # Adapt the title of the plot
988. title = element_text(margin = margin(b = 0.5, unit = "cm"))) +
989. # Set the title
990. ggtitle("Full distribution of responses")
992. # Load overrepresentation helper
993. source("./src/overrepresentation_plot.R")
995. # Make a 4-plot solution with the distribution in the top left corner and some overrepresentation lots in the others
996. # GGarrange from egg package handles legends magnificiently
997. egg::ggarrange(
998. # Full distributions
999. full_response_patterns.pl + theme(legend.box = "vertical"),
1000. # Overrepresentation plot for 'Relative' error pattern
1001. overrepresentation.d %>% overrepresentation_plot("Relative", plot_title = "Proportion of errors of type\nCorrect relative to the line"),
1002. # Overrepresentation plot for 'Count of units' error pattern
1003. overrepresentation.d %>% overrepresentation_plot("Counting_units", plot_title = "Proportion of errors of type\nCounting the numerator"),
1004. # Overrepresentation plot for 'Decimal' error pattern
1005. overrepresentation.d %>% overrepresentation_plot("Decimal", plot_title = "Proportion of errors of type\nDecimal reading"),
1006. # Make it 2x2
1007. ncol = 2, nrow = 2
1008. ) %>%
1009. # Save
1010. save_on_disk("Figure_5", path = "./figures/main/SVG/", device = "svg", a4_fraction = c(1, 3/4), indicate_format = F)
1011. ```
1013. ### Inversion
1015. Plot the frequency of answers applying the pattern on the fraction read in the conventional order or in the opposite one.
1017. ```{r inversion, fig.width = 8.27, fig.height = 11.69 / 4}
1018. # Save the data
1019. inversion.d = d %>%
1020. # Count responses -- we want all responses because we want to include non-realised predictions of the patterns
1021. count_possible_responses() %>%
1022. # Match responses with predictions from patterns...
1023. inner_join(
1024. # ...but only keep predictions that are unambiguous and have an unambiguous inverted counterpart
1025. fractions_predictions.d %>% filter(!ambiguous) %>% group_by(stimulus, pattern) %>% filter(n() == 2),
1026. by = c("stimulus", "actual_response" = "prediction")
1027. ) %>%
1028. # We no longer care about the actual value of the positions
1029. select(-actual_response) %>%
1030. # Compute proportions of responses on target for each position
1031. group_by(stimulus, pattern) %>% mutate(proportion = n / sum(n), proportion_half_ci = half_binomial_ci(proportion, sum(n))) %>% ungroup() %>%
1032. # Turn inverted into a label to make a cleaner pivot_wider
1033. mutate(label = if_else(inverted, "_inverted", "_direct"), .keep = "unused") %>%
1034. # Pivot so that every pattern-stimulus pair is matched with its _direct and _inverted proportions
1035. pivot_wider(names_from = label, values_from = c(proportion, n, proportion_half_ci), names_sep = "") %>%
1036. # The confidence interval is the same for direct and inversion matches, so remove the duplicate
1037. select(-proportion_half_ci_inverted) %>% dplyr::rename(proportion_half_ci = proportion_half_ci_direct) %>%
1038. # Retrieve useful data about the stimuli
1039. inner_join(d %>% select(stimulus, response, numerator, denominator) %>% unique, by = c("stimulus")) %>%
1040. # Reorder stimuli for comparison on graphs...
1041. arrange(pattern, !response %in% c(0.5, 2),
1042. pmax(denominator, numerator), pmin(denominator, numerator), numerator, denominator) %>%
1043. # And save this order into a dedicated column
1044. group_by(pattern) %>% mutate(order = 1:n()) %>% ungroup
1046. ```
1048. Are there differences in inversion across stimuli for some patterns?
1050. ```{r inversion difference between stimuli}
1051. inversion.d %>%
1052. # Chi-square proportion test over responses matching the direct pattern vs. those matching any of the two
1053. group_by(pattern) %>% summarise(distribution_test = prop.test(n_direct, n_direct + n_inverted) %>% (broom::tidy), .groups = "drop") %>% unnest(cols = c(distribution_test)) %>%
1054. # Correct and format statistics nicely
1055. mutate(z_test = str_c("X2(", parameter, ") = ", statistic %>% flexible_round(2)),
1056. # Fdr-correction for p-values
1057. p_value = p.value %>% p.adjust("fdr") %>% p_value_to_APA,
1058. alternative = alternative,
1059. pattern = pattern, .keep = "none")
1060. ```
1062. ```{r inversions plot, fig.width = 8.27, fig.height = 11.69 * 2/3}
1063. # Create a helper function for inversion plot, displaying the bias towards applying the pattern on the target fraction or its inverse.
1064. inversion_plot = function(data, pattern_name, never_inverted_name = "Never\ninverted", always_inverted_name = "Always\ninverted", use_quincunx_for_x_axis = F){
1065. data %>%
1066. # Place stimuli on a W shape -- only if asked for
1067. mutate(stimulus = if_else(order %% 2 == 0 & use_quincunx_for_x_axis, "\n", "") %>% str_c(stimulus)) %>%
1068. # Plot proportion of inversion per stimulus
1069. ggplot(aes(x = reorder(stimulus, order), y = proportion_direct - 0.5, fill = pattern,
1070. ymin = proportion_direct - 0.5 - proportion_half_ci, ymax = proportion_direct - 0.5 + proportion_half_ci)) +
1071. # Proportions are bars
1072. geom_bar(stat = "identity") +
1073. # Add half binominal confidence interval as error bars
1074. geom_errorbar(width = 0.1) +
1075. # Mark the 50-50 line
1076. geom_hline(yintercept = 0, linewidth = 1.5, linetype = 2) +
1077. # Do not stick the bars on the axis as this makes ugly overlaps
1078. coord_cartesian(expand = TRUE) +
1079. # Use color to indicate the pattern
1080. scale_fill_manual(values = patterns.colours) +
1081. # Set the scale to indicate relative preferences
1082. scale_y_continuous(
1083. limits = c(-0.5,0.5), breaks = c(-0.5, 0, 0.5),
1084. labels = c(always_inverted_name, "50-50", never_inverted_name)) +
1085. # Set labels
1086. labs(x = "Target", y = "% of responses without inversion", fill = "Pattern") +
1087. # Clean the theme
1088. theme(
1089. # Remove background elements
1090. panel.background = element_blank(), panel.grid = element_blank(),
1091. # Remove everything but the stimuli from the x-axis
1092. axis.line.x = element_blank(), axis.ticks.x = element_blank(), axis.title.x = element_blank(),
1093. # Remove the ticks and title but restore the line of the y-axis
1094. axis.line.y = element_line(colour = "black"), axis.ticks.y = element_blank(), axis.title.y = element_blank(),
1095. # No need for a legend, fill colors are only evocative
1096. legend.position = "none",
1097. # Add some margins so that everything can be read
1098. plot.margin = unit(c(1, 1, 0.1, 0.1), "cm")
1099. ) +
1100. ggtitle(label = pattern_name)
1101. }
1103. # Make a triple plot with inversion plots on each line
1104. egg::ggarrange(
1105. # One for correct vs.inverted errors
1106. inversion.d %>% filter(pattern == "Normative") %>% inversion_plot("Correct or inverted responses", "Correct", "Inverted") + horizontal_ylab("Bias towards", vjust = 1.2) ,
1107. # One for decimal confusion errors
1108. inversion.d %>% filter(pattern == "Decimal") %>% inversion_plot("Decimal reading errors", "Numerator first", "Denominator first"),
1109. # One for count of tenths errors
1110. inversion.d %>% filter(pattern == "Counting_tenths") %>% inversion_plot("Count of tenths errors", "Count numerator", "Count denominator", use_quincunx_for_x_axis = T)
1111. ) %>%
1112. # Save
1113. save_on_disk("Figure_6", path = "./figures/main/SVG/",
1114. device = "svg", a4_fraction = c(1, 2/3), indicate_format = F)
1115. ```
1117. # Session ending details
1119. ```{r saveFigs}
1120. global_pptx %>% print(target = "./analysis_all_figs.pptx")
1121. sessionInfo()
1122. ```

paper_analysis.Rmd, no license · at the source