OSCR

Exposure to false cardiac feedback alters pain perception and anticipatory cardiac frequency.

Code ↔ Paper

The paper beside its authors' code: matches between them have not been computed for this paper yet.

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The authors' code

R · 947 lines · 36 KB · MIT

  1. # Clear variables
  2. rm(list = ls())
  4. library(readxl)
  5. library(lme4)
  6. library(optimx)
  7. library(lmerTest)
  8. library(effects)
  9. library(emmeans) #
  10. library(brms)
  12. # load data
  13. setwd("C:/Users/Eleonora/OneDrive/MyExperiments/Interoception_exp/2_pain perception/")
  15. data <- read_excel("Data_v3 - 2504.xlsx")
  17. # recode factors
  18. data$Exp <- factor(data$Exp, levels = c("Intero", "Extero"))
  19. data$Feedback <- factor(data$Feedback, levels = c("Congruent", "Slower", "Faster", "No Feedback"))
  21. # remove stim int levels
  22. data = data[data$StimInt != "5", ]
  23. data = data[data$StimInt != "1", ]
  25. # Center Stim intensity (2 = -.5; 3 = 0; 4 = .5)
  26. data$StimInt <- (data$StimInt-3)/2
  28. #################
  29. ### LMM Btw ####
  30. #################
  32. ################### HR BETWEEN ###################
  35. MBtw.HR_N0 = lmer(HR ~ Exp*Feedback*Trial + (1|SSID),data = data, REML = TRUE,
  36. control = lmerControl(optimizer = "optimx", calc.derivs = FALSE,
  37. optCtrl = list(method = "nlminb", starttests = FALSE, kkt = FALSE)))
  39. MBtw.HR_N = lmer(HR ~ Exp*Feedback*Trial + (Trial|SSID),data = data, REML = TRUE,
  40. control = lmerControl(optimizer = "optimx", calc.derivs = FALSE,
  41. optCtrl = list(method = "nlminb", starttests = FALSE, kkt = FALSE)))
  43. model_comparison <- anova(MBtw.HR_N0, MBtw.HR_N)
  44. model_comparison
  46. MBtw.HR_N = lmer(HR ~ Exp*Feedback*Trial + (Trial|SSID),data = data, REML = TRUE,
  47. control = lmerControl(optimizer = "optimx", calc.derivs = FALSE,
  48. optCtrl = list(method = "nlminb", starttests = FALSE, kkt = FALSE)))
  50. # remove outliers (|standardized residuals| > 2 SDs)
  51. MBtw.HRt_N = lmer(HR ~ Exp*Feedback*Trial + (Trial|SSID),data = data,REML = TRUE,
  52. subset = abs(scale(resid(MBtw.HR_N)))<2,
  53. control = lmerControl(optimizer = "optimx", calc.derivs = FALSE,
  54. optCtrl = list(method = "nlminb", starttests = FALSE, kkt = FALSE)))
  56. # residual analysis
  57. qqnorm(resid(MBtw.HRt_N))
  58. qqline(resid(MBtw.HRt_N))
  60. # show model results
  61. anova(MBtw.HRt_N)
  62. summary(MBtw.HRt_N)
  64. # save model results
  65. F_Btw.HR_N = anova(MBtw.HRt_N)
  66. T_Btw.HR_N = data.frame(summary(MBtw.HRt_N)$coefficients)
  68. posthoc_results <- emmeans(MBtw.HRt_N, pairwise ~ Feedback, adjust="fdr")
  70. posthoc_results
  72. # forziamo il download della versione binaria senza compilare
  73. install.packages("DHARMa", type = "binary", dependencies = TRUE, repos = "https://cran.rstudio.com")
  74. library(DHARMa)library(DHARMa)
  75. sim <- simulateResiduals(fittedModel = MBtw.HRt_N, n = 1000)
  76. plotQQunif(sim)
  77. testNormality(sim)
  78. testDispersion(sim)
  82. ################### HR BETWEEN INTERACTIONS ###################
  83. ### FEEDBACK * EXP
  84. library(emmeans)
  85. emmeans_FE <- emmeans(MBtw.HRt_N, ~ Feedback * Exp)
  86. summary(emmeans_FE) #summary(emmeans_FE, infer = c(TRUE, TRUE)) # HA SENSO?
  87. #Contrasti tra feedback entro ciascun esperimento
  88. contrast(emmeans_FE, method = "pairwise", by = "Exp", adjust = "fdr")
  89. #Contrasti di interazione (cioè: la differenza delle differenze)
  90. contrast(emmeans_FE, interaction = "pairwise")
  91. #Contrasti tra feedback entro ciascun esperimento
  92. contrast(emmeans_FE, method = "pairwise", by = "Feedback", adjust = "fdr")
  94. ###FEEDBACK*TRIAL
  95. library(emmeans)
  96. # Calcola la slope del Trial nei due livelli di Exp
  97. emtrends(MBtw.HRt_N, var = "Trial", specs = "Exp")
  98. emtrends(MBtw.HRt_N, var = "Trial", specs = "Exp", infer = c(TRUE, TRUE), adjust = "none")
  99. contrast(emtrends(MBtw.HRt_N, var = "Trial", specs = "Exp"), method = "pairwise")
  104. ## FEEDBACK * TRIAL * EXPERIMENT
  105. anova(MBtw.HRt_N)
  106. summary(MBtw.HRt_N)
  107. library(emmeans)
  108. # Stima degli slopes per ogni condizione Feedback × Exp
  109. # Estrai le slope di HR su Trial
  110. slopes_HR <- emtrends(MBtw.HRt_N, ~ Feedback * Exp, var = "Trial")
  111. # slope significative (aumento/decremento diverso da zero per ogni cond (exp e feedback))
  112. summary(slopes_HR, infer = c(TRUE, TRUE)) # per avere anche t, p e CI, usa questo, non contrast(slopes_HR)
  113. #dentro a ogni esperimento
  114. contrast(slopes_HR, method = "pairwise", by = "Exp", adjust = "fdr")
  115. # differenze delle differenze TRA esperimenti
  116. interaction_contrast <- contrast(slopes_HR, interaction = "pairwise")
  117. summary(interaction_contrast)
  118. #Tutti i confronti possibili tra condizioni (es. Congruent Intero vs Faster Extero)
  119. contrast(slopes_HR, method = "pairwise", adjust = "fdr")
  122. trial_points <- c(-0.5, 0, 0.5)
  124. # Estimated marginal means at specific trial values
  125. em_trial <- emmeans(MBtw.HRt_N, ~ Feedback * Exp | Trial, at = list(Trial = trial_points))
  127. # View predictions
  128. summary(em_trial)
  129. contrast(em_trial, method = "pairwise", by = c("Trial", "Feedback"), adjust = "fdr")
  131. library(emmeans)
  132. library(ggplot2)
  134. # 1. Definisci i tre livelli di interesse per Trial (già centrato)
  135. trial_levels <- c(-0.5, 0, 0.5)
  137. # 2. Calcola gli emmeans a ciascun livello di Trial
  138. em_trial <- emmeans(MBtw.HRt_N, ~ Feedback | Exp * Trial, at = list(Trial = trial_levels))
  140. # 3. Contrasti tra condizioni di Feedback all’interno di ogni Exp × Trial
  141. contrasts_by_trial <- contrast(em_trial, method = "pairwise", by = c("Exp", "Trial"), adjust = "fdr")
  143. # 4. Visualizza i risultati dei contrasti
  144. summary(contrasts_by_trial)
  151. # 1. Cambiamento nel tempo (slopes) per ogni Feedback: confronto Intero vs Extero
  152. slopes_HR <- emtrends(MBtw.HRt_N, ~ Feedback * Exp, var = "Trial")
  154. # Contrasti di interazione: slope_Intero - slope_Extero per ciascun Feedback
  155. slope_diff <- contrast(slopes_HR, interaction = "pairwise")
  156. summary(slope_diff)
  158. # → Se estimate ≠ 0 e p < .05, significa che quel Feedback cambia con Trial in modo diverso tra Exp
  160. # 2. Faster vs Slower a Trial = -0.5, 0, 0.5 dentro ciascun Exp
  161. trial_levels <- c(-0.5, 0, 0.5)
  163. # Calcola gli emmeans per Feedback in ciascun Exp ai trial prescelti
  164. em_trial <- emmeans(
  165. MBtw.HRt_N,
  166. ~ Feedback | Exp * Trial,
  167. at = list(Trial = trial_levels)
  168. )
  170. # Contrasto custom: Faster vs Slower (vettore su livelli di Feedback: Congruent, Slower, Faster, No Feedback)
  171. fs_contrasts <- contrast(
  172. em_trial,
  173. method = list("Faster vs Slower" = c(0, -1, 1, 0)),
  174. by = c("Exp", "Trial"),
  175. adjust = "fdr"
  176. )
  178. # Risultati
  179. summary(fs_contrasts)
  180. # → Per ogni Exp × Trial ottieni stima, t-ratio e p-value di Faster–Slower
  187. ################### LIKERT BETWEEEN ################### qui ok feedback
  190. MBtw.LIK_N0 = lmer(LIK ~ Exp*Feedback*StimInt*Trial + (1|SSID),data = data,
  191. control = lmerControl(optimizer = "optimx", calc.derivs = FALSE,
  192. optCtrl = list(method = "nlminb", starttests = FALSE, kkt = FALSE)))
  194. MBtw.LIK_N = lmer(LIK ~ Exp*Feedback*StimInt*Trial + (StimInt+Trial|SSID),data = data,
  195. control = lmerControl(optimizer = "optimx", calc.derivs = FALSE,
  196. optCtrl = list(method = "nlminb", starttests = FALSE, kkt = FALSE)))
  198. MBtw.LIK_N1 = lmer(LIK ~ Exp*Feedback*StimInt*Trial + (StimInt+Trial+Feedback|SSID),data = data,
  199. control = lmerControl(optimizer = "optimx", calc.derivs = FALSE,
  200. optCtrl = list(method = "nlminb", starttests = FALSE, kkt = FALSE)))
  202. model_comparison <- anova(MBtw.LIK_N, MBtw.LIK_N1)
  203. model_comparison
  208. MBtw.LIK_N = lmer(LIK ~ Exp*Feedback*StimInt*Trial + (StimInt+Trial+Feedback|SSID),data = data,
  209. control = lmerControl(optimizer = "optimx", calc.derivs = FALSE,
  210. optCtrl = list(method = "nlminb", starttests = FALSE, kkt = FALSE)))
  212. # remove outliers (|standardized residuals| > 2 SDs)
  213. MBtw.LIKt_N = lmer(LIK ~ Exp*Feedback*StimInt*Trial + (StimInt+Trial+Feedback|SSID),data = data,
  214. subset = abs(scale(resid(MBtw.LIK_N)))<2,
  215. control = lmerControl(optimizer = "optimx", calc.derivs = FALSE,
  216. optCtrl = list(method = "nlminb", starttests = FALSE, kkt = FALSE)))
  218. # residual analysis
  219. qqnorm(resid(MBtw.LIKt_N))
  220. qqline(resid(MBtw.LIKt_N))
  222. # show model results
  223. anova(MBtw.LIKt_N)
  224. summary(MBtw.LIKt_N)
  226. # save model results
  227. F_Btw.LIK_N = anova(MBtw.LIKt_N)
  228. T_Btw.LIK_N = data.frame(summary(MBtw.LIKt_N)$coefficients)
  230. posthoc_results <- emmeans(MBtw.LIKt_N, pairwise ~ Feedback, adjust="fdr")
  232. posthoc_results
  234. ################### LIKERT BETWEEN INTERACTIONS ###################
  235. # FEEDBACK * EXP
  236. library(emmeans)
  237. emmeans_FE <- emmeans(MBtw.LIKt_N, ~ Feedback * Exp)
  238. summary(emmeans_FE)
  239. #Contrasti tra feedback entro ciascun esperimento
  240. contrast(emmeans_FE, method = "pairwise", by = "Exp", adjust = "fdr")
  241. #Contrasti di interazione (cioè: la differenza delle differenze)
  242. contrast(emmeans_FE, interaction = "pairwise")
  244. # STIMINT*TRIAL
  245. slopes_stimint_at <- emtrends(MBtw.LIKt_N, ~ StimInt, var = "Trial", at = list(StimInt = c(-0.5, 0, 0.5)))
  246. summary(slopes_stimint_at)
  247. test(slopes_stimint_at)
  248. contrast(slopes_stimint_at, method = "pairwise", adjust = "fdr")
  251. #FEEDBACK * STIMINT
  252. em_stimint <- emmeans(MBtw.LIKt_N, ~ Feedback | StimInt, at = list(StimInt = c(-0.5, 0, 0.5)))
  253. summary(em_stimint, infer = c(TRUE, TRUE))
  254. # Pairwise comparisons tra livelli di feedback, per ogni StimInt
  255. contrast(em_stimint, method = "pairwise", adjust = "fdr")
  256. # 1. Calcola le slope (trend) di StimInt per ciascun livello di Feedback
  257. slopes_feedback <- emtrends(MBtw.LIKt_N, ~ Feedback, var = "StimInt")
  258. # 2. Riassunto delle stime: slope, SE, t, p, CI
  259. summary(slopes_feedback)
  260. test(slopes_feedback)
  261. # Contrasti sulle slope di StimInt tra condizioni di feedback
  262. contrast(slopes_feedback, method = "pairwise", adjust = "fdr")
  265. ###FEEDBACK*TRIAL
  266. # Calcolo delle slope dell’effetto di Trial all’interno di ogni livello di Feedback
  267. slopes_trial_by_feedback <- emtrends(MBtw.LIKt_N, ~ Feedback, var = "Trial")
  268. summary(slopes_trial_by_feedback, infer = c(TRUE, TRUE)) # per avere anche t, p e CI
  269. pairs(slopes_trial_by_feedback, adjust = "fdr")
  270. # Contrasti tra le pendenze di Trial tra i diversi livelli di Feedback
  271. trial_slope_contrasts <- pairs(slopes_trial_by_feedback, adjust = "fdr")
  272. summary(trial_slope_contrasts)
  273. plot(slopes_trial_by_feedback)
  274. ## what happens at each trial level
  275. feedback_by_trial_levels <- emmeans(MBtw.LIKt_N, ~ Feedback | Trial,
  276. at = list(Trial = c(-0.5, 0, 0.5)))
  277. pairs(feedback_by_trial_levels, adjust = "fdr")
  279. ################### NPS BETWEEN ###################
  281. MBtw.VAS_NO = lmer(VAS ~ Exp*Feedback*StimInt*Trial +(1|SSID),data = data,
  282. control = lmerControl(optimizer = "optimx", calc.derivs = FALSE,
  283. optCtrl = list(method = "nlminb", starttests = FALSE, kkt = FALSE)))
  284. MBtw.VAS_N = lmer(VAS ~ Exp*Feedback*StimInt*Trial +(StimInt+Trial|SSID),data = data,
  285. control = lmerControl(optimizer = "optimx", calc.derivs = FALSE,
  286. optCtrl = list(method = "nlminb", starttests = FALSE, kkt = FALSE)))
  288. model_comparison <- anova(MBtw.VAS_NO, MBtw.VAS_N)
  289. model_comparison
  291. MBtw.VAS_N = lmer(VAS ~ Exp*Feedback*StimInt*Trial +(StimInt+Trial|SSID),data = data,
  292. control = lmerControl(optimizer = "optimx", calc.derivs = FALSE,
  293. optCtrl = list(method = "nlminb", starttests = FALSE, kkt = FALSE)))
  294. # remove outliers (|standardized residuals| > 2 SDs)
  295. MBtw.VASt_N = lmer(VAS ~ Exp*Feedback*StimInt*Trial +(StimInt+Trial|SSID),data = data,
  296. subset = abs(scale(resid(MBtw.VAS_N)))<2,
  297. control = lmerControl(optimizer = "optimx", calc.derivs = FALSE,
  298. optCtrl = list(method = "nlminb", starttests = FALSE, kkt = FALSE)))
  299. # residual analysis
  300. qqnorm(resid(MBtw.VASt_N))
  301. qqline(resid(MBtw.VASt_N))
  303. # show model results
  304. anova(MBtw.VASt_N)
  305. summary(MBtw.VASt_N)
  307. # save model results
  308. F_Btw.VAS_N = anova(MBtw.VASt_N)
  309. T_Btw.VAS_N = data.frame(summary(MBtw.VASt_N)$coefficients)
  311. posthoc_results <- emmeans(MBtw.VASt_N, pairwise ~ Feedback, adjust="fdr")
  313. posthoc_results
  315. ################### NPS BETWEEN INTERACTIONS ###################
  317. # FEEDBACK * EXP
  318. library(emmeans)
  319. emmeans_FE <- emmeans(MBtw.VASt_N, ~ Feedback * Exp)
  320. summary(emmeans_FE)
  321. #Contrasti tra feedback entro ciascun esperimento
  322. contrast(emmeans_FE, method = "pairwise", by = "Exp", adjust = "fdr")
  323. #Contrasti di interazione (cioè: la differenza delle differenze)
  324. contrast(emmeans_FE, interaction = "pairwise")
  326. # STIMINT*TRIAL
  327. slopes_stimint_at <- emtrends(MBtw.VASt_N, ~ StimInt, var = "Trial", at = list(StimInt = c(-0.5, 0, 0.5)))
  328. summary(slopes_stimint_at)
  329. test(slopes_stimint_at)
  330. contrast(slopes_stimint_at, method = "pairwise", adjust = "fdr")
  333. #FEEDBACK * STIMINT
  334. em_stimint <- emmeans(MBtw.VASt_N, ~ Feedback | StimInt, at = list(StimInt = c(-0.5, 0, 0.5)))
  335. summary(em_stimint, infer = c(TRUE, TRUE))
  336. # Pairwise comparisons tra livelli di feedback, per ogni StimInt
  337. contrast(em_stimint, method = "pairwise", adjust = "fdr")
  338. # 1. Calcola le slope (trend) di StimInt per ciascun livello di Feedback
  339. slopes_feedback <- emtrends(MBtw.VASt_N, ~ Feedback, var = "StimInt")
  340. # 2. Riassunto delle stime: slope, SE, t, p, CI
  341. summary(slopes_feedback)
  342. test(slopes_feedback)
  343. # Contrasti sulle slope di StimInt tra condizioni di feedback
  344. contrast(slopes_feedback, method = "pairwise", adjust = "fdr")
  347. ###FEEDBACK*TRIAL
  348. # Calcolo delle slope dell’effetto di Trial all’interno di ogni livello di Feedback
  349. slopes_trial_by_feedback <- emtrends(MBtw.VASt_N, ~ Feedback, var = "Trial")
  350. summary(slopes_trial_by_feedback, infer = c(TRUE, TRUE)) # per avere anche t, p e CI
  351. pairs(slopes_trial_by_feedback, adjust = "fdr")
  352. # Contrasti tra le pendenze di Trial tra i diversi livelli di Feedback
  353. trial_slope_contrasts <- pairs(slopes_trial_by_feedback, adjust = "fdr")
  354. summary(trial_slope_contrasts)
  355. plot(slopes_trial_by_feedback)
  356. ## what happens at each trial level
  357. feedback_by_trial_levels <- emmeans(MBtw.VASt_N, ~ Feedback | Trial,
  358. at = list(Trial = c(-0.5, 0, 0.5)))
  359. pairs(feedback_by_trial_levels, adjust = "fdr")
  362. ###################
  363. ### LMM Intero ####
  364. ###################
  366. data1 = data[data$Exp=="Intero", ]
  369. ################### HR Interoceptive ###################
  371. MInt.HRN0 = lmer(HR ~ Feedback*Trial + (1|SSID),data = data1, REML = TRUE,
  372. control = lmerControl(optimizer = "optimx", calc.derivs = FALSE,
  373. optCtrl = list(method = "nlminb", starttests = FALSE, kkt = FALSE)))
  375. MInt.HRN = lmer(HR ~ Feedback*Trial + (Trial|SSID),data = data1, REML = TRUE,
  376. control = lmerControl(optimizer = "optimx", calc.derivs = FALSE,
  377. optCtrl = list(method = "nlminb", starttests = FALSE, kkt = FALSE)))
  379. model_comparison <- anova(MInt.HRN0, MInt.HRN)
  380. model_comparison
  382. # remove outliers (|standardized residuals| > 2 SDs)
  383. MInt.HRt = lmer(HR ~ Trial*Feedback + (Trial|SSID),data = data1, REML = TRUE,
  384. subset = abs(scale(resid(MInt.HRN)))<2,
  385. control = lmerControl(optimizer = "optimx", calc.derivs = FALSE,
  386. optCtrl = list(method = "nlminb", starttests = FALSE, kkt = FALSE)))
  390. # residual analysis
  391. qqnorm(resid(MInt.HRt))
  392. qqline(resid(MInt.HRt))
  394. # show model results
  395. anova(MInt.HRt)
  396. summary(MInt.HRt)
  398. # save model results
  399. F_Int.HR = anova(MInt.HRt)
  400. T_Int.HR = data.frame(summary(MInt.HRt)$coefficients)
  402. posthoc_results <- emmeans(MInt.HRt, pairwise ~ Feedback, adjust="fdr")
  404. posthoc_results
  406. ################### HR Interoceptive INTERACTIONS ###################
  407. ## FEEDBACK * TRIAL
  408. # Calcolo delle slope dell’effetto di Trial all’interno di ogni livello di Feedback
  409. slopes_trial_by_feedback <- emtrends(MInt.HRt, ~ Feedback, var = "Trial")
  410. summary(slopes_trial_by_feedback, infer = c(TRUE, TRUE)) # per avere anche t, p e CI
  411. pairs(slopes_trial_by_feedback, adjust = "fdr")
  414. # Contrasti tra le pendenze di Trial tra i diversi livelli di Feedback
  415. trial_slope_contrasts <- pairs(slopes_trial_by_feedback, adjust = "fdr")
  416. summary(trial_slope_contrasts)
  418. plot(slopes_trial_by_feedback)
  420. ## what happens at each trial level
  421. feedback_by_trial_levels <- emmeans(MInt.HRt, ~ Feedback | Trial,
  422. at = list(Trial = c(-0.5, 0, 0.5)))
  423. feedback_by_trial_levels
  424. pairs(feedback_by_trial_levels, adjust = "fdr")
  427. ################### LIKERT Interoceptive ################### qui ok feedback
  430. MInt.LIK1 = lmer(LIK ~ Trial*Feedback*StimInt + (1|SSID),data = data1, REML = TRUE,
  431. control = lmerControl(optimizer = "optimx", calc.derivs = FALSE,
  432. optCtrl = list(method = "nlminb", starttests = FALSE, kkt = FALSE)))
  434. MInt.LIK = lmer(LIK ~ Trial*Feedback*StimInt + (StimInt+Trial|SSID),data = data1, REML = TRUE,
  435. control = lmerControl(optimizer = "optimx", calc.derivs = FALSE,
  436. optCtrl = list(method = "nlminb", starttests = FALSE, kkt = FALSE)))
  438. model_comparison <- anova(MInt.LIK1, MInt.LIK)
  439. model_comparison
  441. # remove outliers (|standardized residuals| > 2 SDs)
  442. MInt.LIKt = lmer(LIK ~ Trial*Feedback*StimInt + (StimInt+Trial|SSID),data = data1, REML = TRUE,
  443. subset = abs(scale(resid(MInt.LIK)))<2,
  444. control = lmerControl(optimizer = "optimx", calc.derivs = FALSE,
  445. optCtrl = list(method = "nlminb", starttests = FALSE, kkt = FALSE)))
  447. # residual analysis
  448. qqnorm(resid(MInt.LIKt))
  449. qqline(resid(MInt.LIKt))
  451. # show model results
  452. anova(MInt.LIKt)
  453. summary(MInt.LIKt)
  455. # save model results
  456. F_Int.LIK = anova(MInt.LIKt)
  457. T_Int.LIK = data.frame(summary(MInt.LIKt)$coefficients)
  459. posthoc_results <- emmeans(MInt.LIKt, pairwise ~ Feedback, adjust="fdr")
  461. posthoc_results
  463. ################### LIKERT Interoceptive INTERACTIONS ###################
  465. #FEEDBACK * STIMINT
  466. em_stimint <- emmeans(MInt.LIKt, ~ Feedback | StimInt, at = list(StimInt = c(-0.5, 0, 0.5)))
  467. summary(em_stimint, infer = c(TRUE, TRUE))
  469. # Pairwise comparisons tra livelli di feedback, per ogni StimInt
  470. contrast(em_stimint, method = "pairwise", adjust = "fdr")
  472. # 1. Calcola le slope (trend) di StimInt per ciascun livello di Feedback
  473. slopes_feedback <- emtrends(MInt.LIKt, ~ Feedback, var = "StimInt")
  474. # 2. Riassunto delle stime: slope, SE, t, p, CI
  475. summary(slopes_feedback)
  476. test(slopes_feedback)
  477. # Contrasti sulle slope di StimInt tra condizioni di feedback
  478. contrast(slopes_feedback, method = "pairwise", adjust = "fdr")
  481. ## FEEDBACK*TRIAL
  482. # Calcolo delle slope dell’effetto di Trial all’interno di ogni livello di Feedback
  483. slopes_trial_by_feedback <- emtrends(MInt.LIKt, ~ Feedback, var = "Trial")
  484. summary(slopes_trial_by_feedback, infer = c(TRUE, TRUE)) # per avere anche t, p e CI
  485. pairs(slopes_trial_by_feedback, adjust = "fdr")
  486. # Contrasti tra le pendenze di Trial tra i diversi livelli di Feedback
  487. trial_slope_contrasts <- pairs(slopes_trial_by_feedback, adjust = "fdr")
  488. summary(trial_slope_contrasts)
  489. #plot(slopes_trial_by_feedback)
  490. ## what happens at each trial level
  491. feedback_by_trial_levels <- emmeans(MInt.LIKt, ~ Feedback | Trial,
  492. at = list(Trial = c(-0.5, 0, 0.5)))
  493. pairs(feedback_by_trial_levels, adjust = "fdr")
  494. # STIMINT * TRIAL
  495. slopes_stimint_at <- emtrends(MInt.LIKt, ~ StimInt, var = "Trial", at = list(StimInt = c(-0.5, 0, 0.5)))
  496. summary(slopes_stimint_at, infer = c(TRUE, TRUE)) # p values diverso da zero?
  497. #differenze delle differenze
  498. contrast(slopes_stimint_at, method = "pairwise", adjust = "fdr")
  501. ################### NPS Interoceptive ###################
  504. MInt.VAS1 = lmer(VAS ~ Trial*Feedback*StimInt + (1|SSID),data = data1, REML = TRUE,
  505. control = lmerControl(optimizer = "optimx", calc.derivs = FALSE,
  506. optCtrl = list(method = "nlminb", starttests = FALSE, kkt = FALSE)))
  509. MInt.VAS = lmer(VAS ~ Trial*Feedback*StimInt + (StimInt+Trial|SSID),data = data1, REML = TRUE,
  510. control = lmerControl(optimizer = "optimx", calc.derivs = FALSE,
  511. optCtrl = list(method = "nlminb", starttests = FALSE, kkt = FALSE)))
  513. model_comparison <- anova(MInt.VAS1, MInt.VAS)
  514. model_comparison
  516. # remove outliers (|standardized residuals| > 2 SDs)
  517. MInt.VASt = lmer(VAS ~ Trial*Feedback*StimInt + (StimInt+Trial|SSID),data = data1, REML = TRUE,
  518. subset = abs(scale(resid(MInt.VAS)))<2,
  519. control = lmerControl(optimizer = "optimx", calc.derivs = FALSE,
  520. optCtrl = list(method = "nlminb", starttests = FALSE, kkt = FALSE)))
  522. # residual analysis
  523. qqnorm(resid(MInt.VASt))
  524. qqline(resid(MInt.VASt))
  526. # show model results
  527. anova(MInt.VASt)
  528. summary(MInt.VASt)
  530. # save model results
  531. F_Int.VAS = anova(MInt.VASt)
  532. T_Int.VAS = data.frame(summary(MInt.VASt)$coefficients)
  534. posthoc_results <- emmeans(MInt.VASt, pairwise ~ Feedback, adjust="fdr")
  536. posthoc_results
  538. ################### NPS Interoceptive INTERACTIONS ###################
  540. # FEEDBACK*STIMINT
  541. # Estimated marginal means di VAS per ciascun livello di Feedback, ai tre livelli di StimInt
  542. em_stimint <- emmeans(MInt.VASt, ~ Feedback | StimInt, at = list(StimInt = c(-0.5, 0, 0.5)))
  543. summary(em_stimint, infer = c(TRUE, TRUE))
  544. # Pairwise comparisons dentro a ogni StimInt
  545. contrast(em_stimint, method = "pairwise", by = "StimInt", adjust = "fdr")
  546. slopes_feedback <- emtrends(MInt.VASt, ~ Feedback, var = "StimInt")
  547. # 2. Riassunto delle stime: slope, SE, t, p, CI
  548. summary(slopes_feedback, infer = c(TRUE, TRUE))
  549. # Contrasti sulle slope di StimInt tra condizioni di feedback
  550. contrast(slopes_feedback, method = "pairwise", adjust = "fdr")
  554. ###FEEDBACK * TRIAL
  555. # Calcolo delle slope dell’effetto di Trial all’interno di ogni livello di Feedback
  556. slopes_trial_by_feedback <- emtrends(MInt.VASt, ~ Feedback, var = "Trial")
  557. summary(slopes_trial_by_feedback, infer = c(TRUE, TRUE)) # per avere anche t, p e CI
  558. pairs(slopes_trial_by_feedback, adjust = "fdr")
  559. # Contrasti tra le pendenze di Trial tra i diversi livelli di Feedback
  560. trial_slope_contrasts <- pairs(slopes_trial_by_feedback, adjust = "fdr")
  561. summary(trial_slope_contrasts)
  562. #plot(slopes_trial_by_feedback)
  563. ## what happens at each trial level
  564. feedback_by_trial_levels <- emmeans(MInt.VASt, ~ Feedback | Trial,
  565. at = list(Trial = c(-0.5, 0, 0.5)))
  566. pairs(feedback_by_trial_levels, adjust = "fdr")
  568. ###################
  569. ### LMM Extero ####
  570. ###################
  572. data2 = data[data$Exp=="Extero", ]
  574. ################### HR exteroceptive ###################
  576. MExt.HR0 = lmer(HR ~ Trial+Feedback + (1|SSID),data = data2, REML = TRUE,
  577. control = lmerControl(optimizer = "optimx", calc.derivs = FALSE,
  578. optCtrl = list(method = "nlminb", starttests = FALSE, kkt = FALSE)))
  580. MExt.HR = lmer(HR ~ Trial*Feedback + (Trial|SSID),data = data2, REML = TRUE,
  581. control = lmerControl(optimizer = "optimx", calc.derivs = FALSE,
  582. optCtrl = list(method = "nlminb", starttests = FALSE, kkt = FALSE)))
  584. model_comparison <- anova(MExt.HR0, MExt.HR)
  585. model_comparison
  588. # remove outliers (|standardized residuals| > 2 SDs)
  589. MExt.HRt = lmer(HR ~ Trial*Feedback + (Trial|SSID),data = data2, REML = TRUE,
  590. subset = abs(scale(resid(MExt.HR)))<2,
  591. control = lmerControl(optimizer = "optimx", calc.derivs = FALSE,
  592. optCtrl = list(method = "nlminb", starttests = FALSE, kkt = FALSE)))
  594. # residual analysis
  595. qqnorm(resid(MExt.HRt))
  596. qqline(resid(MExt.HRt))
  598. # show model results
  599. anova(MExt.HRt)
  600. summary(MExt.HRt)
  602. # save model results
  603. F_Ext.HR = anova(MExt.HRt)
  604. T_Ext.HR = data.frame(summary(MExt.HRt)$coefficients)
  607. posthoc_results <- emmeans(MExt.HRt, pairwise ~ Feedback, adjust="fdr")
  609. posthoc_results
  611. ################### likert exteroceptive ###################
  612. MExt.LIK0 = lmer(LIK ~ Trial*Feedback*StimInt + (1|SSID),data = data2, REML = TRUE,
  613. control = lmerControl(optimizer = "optimx", calc.derivs = FALSE,
  614. optCtrl = list(method = "nlminb", starttests = FALSE, kkt = FALSE)))
  616. MExt.LIK = lmer(LIK ~ Trial*Feedback*StimInt + (StimInt+Trial|SSID),data = data2, REML = TRUE,
  617. control = lmerControl(optimizer = "optimx", calc.derivs = FALSE,
  618. optCtrl = list(method = "nlminb", starttests = FALSE, kkt = FALSE)))
  620. model_comparison <- anova(MExt.LIK0, MExt.LIK)
  621. model_comparison
  624. # remove outliers (|standardized residuals| > 2 SDs)
  625. MExt.LIKt = lmer(LIK ~ Trial*Feedback*StimInt + (StimInt+Trial|SSID),data = data2, REML = TRUE,
  626. subset = abs(scale(resid(MExt.LIK)))<2,
  627. control = lmerControl(optimizer = "optimx", calc.derivs = FALSE,
  628. optCtrl = list(method = "nlminb", starttests = FALSE, kkt = FALSE)))
  630. # residual analysis
  631. qqnorm(resid(MExt.LIKt))
  632. qqline(resid(MExt.LIKt))
  634. # show model results
  635. anova(MExt.LIKt)
  636. summary(MExt.LIKt)
  638. # save model results
  639. F_Ext.LIK = anova(MExt.LIKt)
  640. T_Ext.LIK = data.frame(summary(MExt.LIKt)$coefficients)
  642. posthoc_results <- emmeans(MExt.LIKt, pairwise ~ Feedback, adjust="fdr")
  644. posthoc_results
  647. ### EXTEROCEPTIVE likert INTERACTIONS
  648. #STIMULUSINTESITY * FEEDBACK
  649. em_stimint <- emmeans(MExt.LIKt, ~ Feedback | StimInt, at = list(StimInt = c(-0.5, 0, 0.5)))
  650. summary(em_stimint)
  651. # Pairwise comparisons tra livelli di feedback, per ogni StimInt
  652. contrast(em_stimint, method = "pairwise", adjust = "fdr")
  653. # 1. Calcola le slope (trend) di StimInt per ciascun livello di Feedback
  654. slopes_feedback <- emtrends(MExt.LIKt, ~ Feedback, var = "StimInt")
  655. # 2. Riassunto delle stime: slope, SE, t, p, CI
  656. summary(slopes_feedback)
  657. test(slopes_feedback)
  658. # Contrasti sulle slope di StimInt tra condizioni di feedback
  659. contrast(slopes_feedback, method = "pairwise", adjust = "fdr")
  663. ### FEEDBACK * TRIAL
  664. # Calcolo delle slope dell’effetto di Trial all’interno di ogni livello di Feedback
  665. slopes_trial_by_feedback <- emtrends(MExt.LIKt, ~ Feedback, var = "Trial")
  666. summary(slopes_trial_by_feedback, infer = c(TRUE, TRUE)) # per avere anche t, p e CI
  667. # Contrasti tra le pendenze di Trial tra i diversi livelli di Feedback
  668. pairs(slopes_trial_by_feedback, adjust = "fdr")
  669. #plot(slopes_trial_by_feedback)
  670. ## what happens at each trial level
  671. feedback_by_trial_levels <- emmeans(MExt.LIKt, ~ Feedback | Trial,
  672. at = list(Trial = c(-0.5, 0, 0.5)))
  673. feedback_by_trial_levels
  674. pairs(feedback_by_trial_levels, adjust = "fdr")
  676. # STIMINT * TRIAL
  677. slopes_stimint_at <- emtrends(MExt.LIKt, ~ StimInt, var = "Trial", at = list(StimInt = c(-0.5, 0, 0.5)))
  678. summary(slopes_stimint_at, infer = c(TRUE, TRUE)) # p values diverso da zero?
  679. #differenze delle differenze
  680. contrast(slopes_stimint_at, method = "pairwise", adjust = "fdr")
  684. ################### VAS exteroceptive ###################
  685. MExt.VAS0 = lmer(VAS ~ Trial*Feedback*StimInt + (1|SSID),data = data2, REML = TRUE,
  686. control = lmerControl(optimizer = "optimx", calc.derivs = FALSE,
  687. optCtrl = list(method = "nlminb", starttests = FALSE, kkt = FALSE)))
  689. ## VAS
  690. MExt.VAS = lmer(VAS ~ Trial*Feedback*StimInt + (Trial+StimInt|SSID),data = data2, REML = TRUE,
  691. control = lmerControl(optimizer = "optimx", calc.derivs = FALSE,
  692. optCtrl = list(method = "nlminb", starttests = FALSE, kkt = FALSE)))
  695. model_comparison <- anova(MExt.VAS0, MExt.VAS)
  696. model_comparison
  699. # remove outliers (|standardized residuals| > 2 SDs)
  700. MExt.VASt = lmer(VAS ~ Trial*Feedback*StimInt + (Trial+StimInt|SSID),data = data2, REML = TRUE,
  701. subset = abs(scale(resid(MExt.VAS)))<2,
  702. control = lmerControl(optimizer = "optimx", calc.derivs = FALSE,
  703. optCtrl = list(method = "nlminb", starttests = FALSE, kkt = FALSE)))
  705. # residual analysis
  706. qqnorm(resid(MExt.VASt))
  707. qqline(resid(MExt.VASt))
  709. # show model results
  710. anova(MExt.VASt)
  711. summary(MExt.VASt)
  713. # save model results
  714. F_Ext.VAS = anova(MExt.VASt)
  715. T_Ext.VAS = data.frame(summary(MExt.VASt)$coefficients)
  718. posthoc_results <- emmeans(MExt.VASt, pairwise ~ Feedback, adjust="fdr")
  720. posthoc_results
  723. ### VAS
  725. # FEEDBACK*STIMINT i
  726. # Estimated marginal means di VAS per ciascun livello di Feedback, ai tre livelli di StimInt
  727. em_stimint <- emmeans(MExt.VASt, ~ Feedback | StimInt, at = list(StimInt = c(-0.5, 0, 0.5)))
  728. summary(em_stimint, infer = c(TRUE, TRUE))
  729. # Pairwise comparisons dentro a ogni StimInt
  730. contrast(em_stimint, method = "pairwise", by = "StimInt", adjust = "fdr")
  731. ## per differenze con zero e contrasti tra feedback (quanto l'increase del main effect stimint e' diverso per ogni feed)
  732. # 1. Calcola le slope (trend) di StimInt per ciascun livello di Feedback
  733. slopes_feedback <- emtrends(MExt.VASt, ~ Feedback, var = "StimInt")
  734. # 2. Riassunto delle stime: slope, SE, t, p, CI
  735. summary(slopes_feedback, infer = c(TRUE, TRUE))
  736. # Contrasti sulle slope di StimInt tra condizioni di feedback
  737. contrast(slopes_feedback, method = "pairwise", adjust = "fdr")
  739. # STIMINT * TRIAL
  740. slopes_stimint_at <- emtrends(MExt.VASt, ~ StimInt, var = "Trial", at = list(StimInt = c(-0.5, 0, 0.5)))
  741. summary(slopes_stimint_at, infer = c(TRUE, TRUE)) # p values diverso da zero?
  742. #differenze delle differenze
  743. contrast(slopes_stimint_at, method = "pairwise", adjust = "fdr")
  745. ## FEEDBACK * TRIAL
  746. # Calcolo delle slope dell’effetto di Trial all’interno di ogni livello di Feedback
  747. slopes_trial_by_feedback <- emtrends(MExt.VASt, ~ Feedback, var = "Trial")
  748. summary(slopes_trial_by_feedback, infer = c(TRUE, TRUE)) # per avere anche t, p e CI
  749. pairs(slopes_trial_by_feedback, adjust = "fdr")
  750. # Contrasti tra le pendenze di Trial tra i diversi livelli di Feedback
  751. trial_slope_contrasts <- pairs(slopes_trial_by_feedback, adjust = "fdr")
  752. summary(trial_slope_contrasts)
  753. plot(slopes_trial_by_feedback)
  754. ## what happens at each trial level
  755. feedback_by_trial_levels <- emmeans(MExt.VASt, ~ Feedback | Trial,
  756. at = list(Trial = c(-0.5, 0, 0.5)))
  757. feedback_by_trial_levels
  758. pairs(feedback_by_trial_levels, adjust = "fdr")
  760. ######################## pain response
  763. model_comparison <- anova(painrespBTW1, painrespBTW)
  764. model_comparison
  766. painrespBTW2 = lmer(HRPR ~ Trial*Feedback*StimInt*Exp + (Trial+StimInt|SSID),data = data, REML = TRUE,
  767. control = lmerControl(optimizer = "optimx", calc.derivs = FALSE,
  768. optCtrl = list(method = "nlminb", starttests = FALSE, kkt = FALSE)))
  772. painrespBTW1 = lmer(HRPR ~ Trial*Feedback*StimInt*Exp + (Trial|SSID),data = data, REML = TRUE,
  773. control = lmerControl(optimizer = "optimx", calc.derivs = FALSE,
  774. optCtrl = list(method = "nlminb", starttests = FALSE, kkt = FALSE)))
  776. painrespBTW = lmer(HRPR ~ Trial*Feedback*StimInt*Exp + (1|SSID),data = data, REML = TRUE,
  777. control = lmerControl(optimizer = "optimx", calc.derivs = FALSE,
  778. optCtrl = list(method = "nlminb", starttests = FALSE, kkt = FALSE)))
  780. # remove outliers (|standardized residuals| > 2 SDs)
  781. painrespBTW_n = lmer(HRPR ~ Trial*Feedback*StimInt*Exp + (1|SSID),data = data, REML = TRUE,
  782. subset = abs(scale(resid(painrespBTW)))<2,
  783. control = lmerControl(optimizer = "optimx", calc.derivs = FALSE,
  784. optCtrl = list(method = "nlminb", starttests = FALSE, kkt = FALSE)))
  787. # residual analysis
  788. qqnorm(resid(painrespBTW_n))
  789. qqline(resid(painrespBTW_n))
  791. # show model results
  792. anova(painrespBTW_n)
  793. summary(painrespBTW_n)
  795. # save model results
  796. F_Ext.VAS = anova(painrespBTW_n)
  797. T_Ext.VAS = data.frame(summary(painrespBTW_n)$coefficients)
  799. posthoc_results <- emmeans(F_Ext.VAS, pairwise ~ Exp, adjust="fdr")
  801. posthoc_results
  803. # Calcola i simple slopes di StimInt per ciascun livello di Exp
  804. emtrends(painrespBTW_n, specs = ~Exp, var = "StimInt") %>%
  805. summary(infer = c(TRUE, TRUE))
  807. # Calcolo dei simple slopes dell’effetto di Trial a ciascun livello di StimInt e per ciascun Exp
  808. emtrends(painrespBTW_n, specs = ~Exp | StimInt, var = "Trial", at = list(StimInt = c(-0.5, 0, 0.5))) %>%
  809. summary(infer = c(TRUE, TRUE))
  811. emtrends(painrespBTW_n, specs = ~Exp | Feedback, var = "StimInt") %>%
  812. summary(infer = c(TRUE, TRUE))
  815. # plot effects
  816. # If not installed yet
  817. # install.packages("effects")
  819. library(effects)
  821. # Get the effect of the 3-way interaction
  822. eff <- effect("Feedback:StimInt:Exp", painrespBTW_n)
  824. # Plot
  825. plot(eff, multiline = TRUE, ci.style = "bands",
  826. main = "Feedback × StimInt × Exp Interaction (Model-Based)",
  827. xlab = "StimInt", ylab = "Predicted HRPR")
  830. ####
  831. # If not installed yet
  832. # install.packages("ggeffects")
  834. library(ggeffects)
  835. library(ggplot2)
  837. # Get model-based predictions for the 3-way interaction
  838. preds <- ggpredict(painrespBTW_n, terms = c("StimInt", "Feedback", "Exp"))
  840. # Plot
  841. ggplot(preds, aes(x = x, y = predicted, color = group)) +
  842. geom_line(size = 1) +
  843. geom_point(size = 2) +
  844. geom_ribbon(aes(ymin = conf.low, ymax = conf.high, fill = group), alpha = 0.2, color = NA) +
  845. facet_wrap(~ facet) +
  846. labs(title = "Model-Based Interaction Plot: Feedback × StimInt × Exp",
  847. x = "StimInt",
  848. y = "Predicted HRPR",
  849. color = "Feedback",
  850. fill = "Feedback") +
  851. theme_minimal(base_size = 14)
  854. ###############################################
  855. library(ggeffects)
  856. library(ggplot2)
  858. # Get model-based predictions for the interaction between Trial and Exp
  859. preds <- ggpredict(painrespBTW_n, terms = c("Trial", "Exp"))
  861. # Plot
  862. ggplot(preds, aes(x = x, y = predicted, color = group)) +
  863. geom_line(size = 1) +
  864. geom_point(size = 2) +
  865. geom_ribbon(aes(ymin = conf.low, ymax = conf.high, fill = group), alpha = 0.2, color = NA) +
  866. labs(title = "Model-Based Interaction Plot: Trial × Exp",
  867. x = "Trial",
  868. y = "Predicted HRPR",
  869. color = "Exp",
  870. fill = "Exp") +
  871. theme_minimal(base_size = 14)
  874. #####
  876. # FEEDBACK * EXP
  877. library(emmeans)
  878. emmeans_FE <- emmeans(MBtw.VASt_N, ~ Feedback * Exp)
  879. summary(emmeans_FE)
  880. #Contrasti tra feedback entro ciascun esperimento
  881. contrast(emmeans_FE, method = "pairwise", by = "Exp", adjust = "fdr")
  882. #Contrasti di interazione (cioè: la differenza delle differenze)
  883. contrast(emmeans_FE, interaction = "pairwise")
  887. #########################
  889. correl_model = lmer(VAS ~ HR*Exp*Feedback + (1|SSID),data = data, REML = TRUE,
  890. control = lmerControl(optimizer = "optimx", calc.derivs = FALSE,
  891. optCtrl = list(method = "nlminb", starttests = FALSE, kkt = FALSE)))
  893. # remove outliers (|standardized residuals| > 2 SDs)
  894. correl_model_n = lmer(VAS ~ HR*Exp*Feedback + (1|SSID),data = data, REML = TRUE,
  895. subset = abs(scale(resid(correl_model)))<2,
  896. control = lmerControl(optimizer = "optimx", calc.derivs = FALSE,
  897. optCtrl = list(method = "nlminb", starttests = FALSE, kkt = FALSE)))
  899. # residual analysis
  900. qqnorm(resid(correl_model_n))
  901. qqline(resid(correl_model_n))
  903. # show model results
  904. anova(correl_model_n)
  905. summary(correl_model_n)
  907. # save model results
  908. F_Ext.VAS = anova(correl_model_n)
  909. T_Ext.VAS = data.frame(summary(correl_model_n)$coefficients)
  915. ###FEEDBACK
  916. emmeans_feedback <- emmeans(correl_model_n, ~ Feedback)
  917. # Confronti post hoc a coppie tra le condizioni di Feedback (con correzione Tukey)
  918. pairs(emmeans_feedback, adjust = "fdr")
  921. library(emmeans)
  923. ###FEEDBACK*vas
  924. emtrends(correl_model_n, ~ Feedback, var = "VAS")
  925. pairs(emtrends(correl_model_n, ~ Feedback, var = "VAS"), adjust = "fdr")
  928. # Post hoc: slope di VAS per ciascun tipo di Feedback
  929. em_slopes <- emtrends(correl_model_n, ~ Feedback, var = "VAS")
  930. # Visualizza le slopes
  931. summary(em_slopes)
  932. # Confronti post hoc tra le slopes
  933. pairs(em_slopes, adjust = "fdr")
  940. # slope di VAS → HR per ciascuna combinazione Feedback × Exp
  941. em_triple <- emtrends(correl_model_n, ~ Feedback * Exp, var = "VAS")
  943. # mostra le slopes
  944. summary(em_triple)
  946. # confronti post hoc tra combinazioni
  947. pairs(em_triple, adjust = "fdr")

mixed_models_all.R at commit c1c8721, under MIT · at the source

Overview

Authors: Eleonora Parrotta1,2,3, Patric Bach2, Giovanni Pezzulo4, Andrea Zaccaro3, Mauro Gianni Perrucci3,5,6, Marcello Costantini5,7, Francesca Ferri3,5
ORCID iDs: Eleonora Parrotta, Patric Bach, Giovanni Pezzulo, Andrea Zaccaro
  1. Department of Psychology, Sapienza University of Rome, Rome, Italy
  2. School of Psychology, University of Aberdeen, Aberdeen, United Kingdom
  3. Department of Neuroscience, Imaging and Clinical Sciences, “G. d'Annunzio” University of Chieti-Pescara, Chieti, Italy
  4. Institute of Cognitive Sciences and Technologies, National Research Council, Rome, Italy
  5. Institute for Advanced Biomedical Technologies ‑ ITAB, “G. d'Annunzio” University of Chieti-Pescara, Chieti, Italy
  6. UdA-TechLab, Research Center, University “G. d’Annunzio” of Chieti-Pescara, Chieti, Italy
  7. Department of Psychology, “G. d’Annunzio” University of Chieti-Pescara, Chieti, Italy
Institutions: University of Aberdeen (United Kingdom); University of Chieti-Pescara (Italy); National Research Council (Italy); Institute of Cognitive Sciences and Technologies (Italy); Sapienza University of Rome (Italy); University of Plymouth (United Kingdom); National Academies of Sciences, Engineering, and Medicine (United States)
Journal: eLife, volume 12, article RP90013
Dates: published online 15 June 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.7554/elife.90013 · PMID 42294597 · PMCID PMC13268646 · OpenAlex W4387132705
Open access: gold, a free copy (OpenAlex)
Status: code verified
Categories: human (organism), pain (population), cognitive (subfield)
Methods: Statistics
Keywords: Human
MeSH: Anticipation, Psychological*, Heart Rate*, Pain Perception*, Adult, Female, Humans, Male, Young Adult (* major topic)
Topic: Psychosomatic Disorders and Their Treatments (Psychiatry and Mental health, Medicine), according to OpenAlex
Funding: European Union's Horizon 2020 (10.3030/945539, 10.3030/952215); European Research Council (10.3030/820213, 820213); the PNRR MUR (PE0000013-FAIR); Leverhulme Trust (RPG-2019-248)
Citations: cited by 1 paper (Europe PMC); 146 references in the paper

Abstract

The experience of pain, like other interoceptive processes, has recently been conceptualized in terms of predictive coding and free energy frameworks. In these views, the brain integrates sensory, proprioceptive, and interoceptive signals to generate probabilistic inferences about upcoming events, which shape both the state and the perception of our inner body. Here, we ask whether it is possible to induce pain expectations by providing false faster (vs. slower) acoustic cardiac feedback before administering electrical cutaneous shocks. We test whether these expectations will shape both the perception of pain and the body’s physiological state toward prior predictions. Results confirmed that faster cardiac feedback elicited pain expectations that affected both perceptual pain judgments and the body’s physiological response. Perceptual pain judgments were biased toward the expected level of pain, such that participants illusorily perceived identical noxious stimuli as more intense and unpleasant. Physiological changes mirrored the predicted level of pain, such that participants’ actual cardiac response in anticipation of pain stimuli showed a deceleration in heart rate, in line with the well-known orienting cardiac response in anticipation of threatening stimuli (Experiment 1). In a control experiment, such perceptual and cardiac modulations were dramatically reduced when the feedback reproduced an exteroceptive, instead of interoceptive, cardiac feedback (Experiment 2). These findings show that cardiac perception can be understood as interoceptive inference that modulates both our perception and the physiological state of the body, thereby actively generating the interoceptive and autonomic consequences that have been predicted.

Reproduced under the paper's license (CC BY), from the paper cited above.

Repository

Its files are read in the Code ↔ Paper reader above.

EP171993/painperception

License: MIT
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: c1c8721f0e3d00fefc602ed32b84e24c218c56d3, 3 June 2026
Languages: R (1)
Size: 5 files, 1 script
Software Heritage: not archived
Found in: “Data availability”
Holds: README, license file
Not found: CITATION.cff, environment file, tests, continuous integration, documentation
Tools: brms (1 file), emmeans (1 file), ggplot2 (1 file), lme4 (1 file), lmerTest (1 file)
Availability: 1 check, the latest on 27 September 2026: the link answers
  • 27 September 2026: the link answers
3 files
At the source: github.com/EP171993/painperception

The paper's code and data availability statement is in the Data section.

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

  • 1 repository of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 1 script, each with its path and the digest of its content;
  • no match between paragraphs and code yet;
  • neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

Datasets cited

Data availability

Anonymized raw ECG recordings underlying the results reported in this study are publicly available on OSF at: https://doi.org/10.17605/OSF.IO/5SW3M. Anonymized trial-level behavioural data, ECG-derived physiological measures, and analysis code are publicly available at: https://github.com/EP171993/painperception (copy archived at Parrotta, 2026).

The following dataset was generated:

Parrotta E. 2026. Exposure to false cardiac feedback alters pain perception and anticipatory cardiac frequency. Open Science Framework.

Reproduced under the paper's license (CC BY), from the paper cited above.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 1, 27 September 2026: the first record

Recorded: type, language, journal, volume, pages, dates, 7 authors, 1 keyword, 8 MeSH terms, 4 funders, 141 references.

Cite

This paper

Parrotta, E., Bach, P., Pezzulo, G., Zaccaro, A., Perrucci, M. G., Costantini, M., & Ferri, F. (2026). Exposure to false cardiac feedback alters pain perception and anticipatory cardiac frequency. eLife, 12, RP90013. https://doi.org/10.7554/elife.90013

BibTeX

@article{parrotta2026exposure,
author = {Parrotta, Eleonora and Bach, Patric and Pezzulo, Giovanni and Zaccaro, Andrea and Perrucci, Mauro Gianni and Costantini, Marcello and Ferri, Francesca},
title = {{Exposure to false cardiac feedback alters pain perception and anticipatory cardiac frequency}},
journal = {eLife},
year = {2026},
month = jun,
volume = {12},
pages = {RP90013},
publisher = {eLife Sciences Publications, Ltd},
issn = {2050-084X},
doi = {10.7554/elife.90013},
url = {https://doi.org/10.7554/elife.90013},
pmid = {42294597},
pmcid = {PMC13268646}
}

RIS

TY - JOUR
AU - Parrotta, Eleonora
AU - Bach, Patric
AU - Pezzulo, Giovanni
AU - Zaccaro, Andrea
AU - Perrucci, Mauro Gianni
AU - Costantini, Marcello
AU - Ferri, Francesca
TI - Exposure to false cardiac feedback alters pain perception and anticipatory cardiac frequency
T2 - eLife
J2 - Elife
PY - 2026
DA - 2026/06/15
VL - 12
SP - RP90013
SN - 2050-084X
PB - eLife Sciences Publications, Ltd
DO - 10.7554/elife.90013
UR - https://doi.org/10.7554/elife.90013
LA - en
ER -

CSL-JSON

{
"id": "10.7554/elife.90013",
"type": "article-journal",
"title": "Exposure to false cardiac feedback alters pain perception and anticipatory cardiac frequency",
"container-title": "eLife",
"author": [
{
"family": "Parrotta",
"given": "Eleonora"
},
{
"family": "Bach",
"given": "Patric"
},
{
"family": "Pezzulo",
"given": "Giovanni"
},
{
"family": "Zaccaro",
"given": "Andrea"
},
{
"family": "Perrucci",
"given": "Mauro Gianni"
},
{
"family": "Costantini",
"given": "Marcello"
},
{
"family": "Ferri",
"given": "Francesca"
}
],
"container-title-short": "Elife",
"volume": "12",
"page": "RP90013",
"DOI": "10.7554/elife.90013",
"PMID": "42294597",
"PMCID": "PMC13268646",
"ISSN": "2050-084X",
"publisher": "eLife Sciences Publications, Ltd",
"URL": "https://doi.org/10.7554/elife.90013",
"language": "en",
"issued": {
"date-parts": [
[
2026,
6,
15
]
]
}
}

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.1038/s41467-026-74743-0 [code]
Meta-analytic evidence for distinct neural correlates of conditioned versus verbally induced placebo analgesia.
Journal: Nature communications
In common: lmerTest, lme4, ggplot2, pain, 8 references
[2] doi:10.3389/fnhum.2026.1820376 [code]
The neural dynamics of political socio-pragmatic violations: an ERP study.
Journal: Frontiers in human neuroscience
In common: emmeans, lmerTest, lme4, 1 other tool, 3 references
[3] doi:10.1038/s41598-026-52930-9
Cardiac systole is associated with enhanced go responding in an orthogonalized go/nogo task.
Journal: Scientific reports
In common: cognitive, 6 references
[4] doi:10.1371/journal.pone.0355165 [code]
Pupillary dynamics during hands-off L2 driving and transitions of control under high cognitive load.
Journal: PloS one
In common: brms, emmeans, lmerTest, 2 other tools, cognitive
[5] doi:10.1016/j.isci.2026.116747 [code]
Age and loneliness relate to reduced trust learning and alterations in amygdala function.
Journal: iScience
In common: brms, emmeans, lmerTest, 2 other tools, cognitive
[6] doi:10.1038/s41467-026-73865-9 [code]
Histamine shapes the neurocomputational dynamics of human learning.
Journal: Nature communications
In common: brms, emmeans, lmerTest, 2 other tools, cognitive
[7] doi:10.1523/eneuro.0316-25.2026 [code]
Neural Mechanisms of Self-Generated Action Sequences.
Journal: eNeuro
In common: brms, emmeans, lmerTest, 2 other tools, cognitive
[8] doi:10.1093/nc/niag046 [code]
Awareness of being: a computational neurophenomenological model of mindfulness, mind-wandering, and meta-attentional control.
Journal: Neuroscience of consciousness
In common: 6 references
[9] doi:10.1093/braincomms/fcag120 [code]
The heartbeat evoked potential and the prediction of functional seizure semiology.
Journal: Brain communications
In common: 5 references
[10] doi:10.1038/s41467-026-72916-5 [code]
Precision fMRI reveals that the language network exhibits adult-like left-hemispheric lateralization by 4 years of age.
Journal: Nature communications
In common: brms, emmeans, lmerTest, 2 other tools

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.