RatingAnalysis.R 8.6 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169
  1. # _____ _ _ _ ______
  2. # / ____| | | (_) | | | ____|
  3. # | (___ _ __ ___ ___ | | ___ _ __ __ _ __ _ _ __ __| | | |__ ___ __ _ _ __
  4. # \___ \| '_ ` _ \ / _ \| |/ / | '_ \ / _` | / _` | '_ \ / _` | | __/ _ \/ _` | '__|
  5. # ____) | | | | | | (_) | <| | | | | (_| | | (_| | | | | (_| | | | | __/ (_| | |
  6. # |_____/|_| |_| |_|\___/|_|\_\_|_| |_|\__, | \__,_|_| |_|\__,_| |_| \___|\__,_|_|
  7. # __/ |
  8. # |___/
  9. #Analysis script Fear ratings and US expectancy ratings (day1 and day 2)
  10. #written by Madeleine Mueller, February 2024
  11. library("lme4")
  12. library("car")
  13. library("effects")
  14. library("emmeans")
  15. # __ _ _
  16. # / _| | | (_)
  17. # | |_ ___ __ _ _ __ _ __ __ _| |_ _ _ __ __ _ ___
  18. # | _/ _ \/ _` | '__| | '__/ _` | __| | '_ \ / _` / __|
  19. # | || __/ (_| | | | | | (_| | |_| | | | | (_| \__ \
  20. # |_| \___|\__,_|_| |_| \__,_|\__|_|_| |_|\__, |___/
  21. # _ _ _ _ __/ |
  22. # | | | | | | | | |___/
  23. # | |__ ___ | |_| |__ __| | __ _ _ _ ___
  24. # | '_ \ / _ \| __| '_ \ / _` |/ _` | | | / __|
  25. # | |_) | (_) | |_| | | | | (_| | (_| | |_| \__ \
  26. # |_.__/ \___/ \__|_| |_| \__,_|\__,_|\__, |___/
  27. # __/ |
  28. # |___/
  29. #read in data of fear rating for both days
  30. onstud_fear<-read.csv(file='/Madeleine/Online Studie/analysis/R/FearRatingsOnstudforR.csv',head=TRUE,sep=",")
  31. #all four groups (1=non-smoker, 2=smoker with 6h break after ACQ, 3=smoker with cigarette after ACQ, 4=smoker with cigarette before ACQ)
  32. onstud_fear$group<-as.factor(onstud_fear$group)
  33. #all 4 groups, see above, but smokers from group 2 that did smoke during the 6h smoke break were regrouped into group 3
  34. onstud_fear$regroup<-as.factor(onstud_fear$regroup)
  35. #stimulus: CSP=CS+, CSM=CS-
  36. onstud_fear$stim<-as.factor(onstud_fear$stim)
  37. #time: pre ACQ/pre GEN (depending on day)=1; post ACQ/ post GEN = 2
  38. onstud_fear$time<-as.factor(onstud_fear$prepost)
  39. #smoker(group2,3 and 4)= 1, non-smoker(group1) = 0
  40. onstud_fear$smoker<-as.factor(onstud_fear$smoker)
  41. #subject number
  42. onstud_fear$sub<-onstud_fear$ID
  43. #fear ratings day 1 models with different groups
  44. both.f1s<-(lmer(ratingd1~(1|sub)+stim*time*smoker, data=onstud_fear,control = lmerControl(optimizer = "bobyqa")))
  45. both.f1g<-(lmer(ratingd1~(1|sub)+stim*time*regroup, data=onstud_fear,control = lmerControl(optimizer = "bobyqa")))
  46. #fear ratings day 1 models with different groups
  47. both.f2s<-(lmer(ratingd2~(1|sub)+stim*time*smoker, data=onstud_fear,control = lmerControl(optimizer = "bobyqa")))
  48. both.f2g<-(lmer(ratingd2~(1|sub)+stim*time*regroup, data=onstud_fear,control = lmerControl(optimizer = "bobyqa")))
  49. #anova based on models, include different model name for all results
  50. Anova(both.f1s, type="3", test="F" )
  51. Anova(both.f2s, type="3", test="F" )
  52. #post-hoc tests, include effects in which you are interested
  53. emmeans(both.f2s, list(pairwise ~ stim:time:smoker), adjust = "none")
  54. #plot effects, change effect name and model name as needed
  55. plot(effect("stim*smoker",both.f1s))
  56. plot(effect("stim*time*smoker",both.f2s))
  57. # _ _ _____ _
  58. # | | | |/ ____| | |
  59. # | | | | (___ _____ ___ __ ___ ___| |_ __ _ _ __ ___ _ _
  60. # | | | |\___ \ / _ \ \/ / '_ \ / _ \/ __| __/ _` | '_ \ / __| | | |
  61. # | |__| |____) | | __/> <| |_) | __/ (__| || (_| | | | | (__| |_| |
  62. # \____/|_____/ \___/_/\_\ .__/ \___|\___|\__\__,_|_| |_|\___|\__, |
  63. # | | /_ | | | __/ |
  64. # __| | __ _ _ _ | | |_| |___/
  65. # / _` |/ _` | | | | | |
  66. # | (_| | (_| | |_| | | |
  67. # \__,_|\__,_|\__, | |_|
  68. # __/ |
  69. # |___/
  70. #read in data
  71. data_ex<-read.csv(file='/Madeleine/Online Studie/analysis/R/D1USexponstudR.csv',head=TRUE,sep=",")
  72. #all four groups (1=non-smoker, 2=smoker with 6h break after ACQ, 3=smoker with cigarette after ACQ, 4=smoker with cigarette before ACQ)
  73. data_ex$group<-as.factor(data_ex$group)
  74. #all 4 groups, see above, but smokers from group 2 that did smoke during the 6h smoke break were regrouped into group 3
  75. data_ex$regroup<-as.factor(data_ex$regroup)
  76. #stimulus: CSP=CS+, CSM=CS-
  77. data_ex$stim<-as.factor(data_ex$stim)
  78. # trial 1-12 per stimulus, so in total 24 trials
  79. data_ex$trial<-as.factor(data_ex$trial)
  80. #block 1-3, one block consists of 4 trials per stimulus, so 8 trials in total
  81. data_ex$block<-as.factor(data_ex$block)
  82. #difference rated between pleasantness US and nUS (mean US - mean nUS)
  83. data_ex$usdiff<-as.factor(data_ex$usdiff)
  84. # avoidance of screen when US is presented, higher rating means more avoidance
  85. data_ex$avoid<-as.factor(data_ex$avoid)
  86. # CS-US contingency, 1=not aware of correct CS-US contingency, 2=aware of correct CS-US contingency
  87. data_ex$aware<-as.factor(data_ex$aware)
  88. #smoker(group2,3 and 4)= 1, non-smoker(group1) = 0
  89. data_ex$smoker<-as.factor(data_ex$smoker)
  90. #subject number
  91. data_ex$sub<-data_ex$sub
  92. #########################################
  93. #models
  94. #########################################
  95. USex1a<-lmer(rating~(1|sub)+stim*block*regroup,data = data_ex, control = lmerControl(optimizer = "bobyqa"))
  96. Anova(USex1a, type="3", test="F" )
  97. USex1sm<-lmer(rating~(1|sub)+stim*block*smoker,data = data_ex,control = lmerControl(optimizer = "bobyqa"))
  98. Anova(USex1sm, type="3", test="F" )
  99. emmeans(USex1a, list(pairwise ~ stim), adjust = "none")
  100. plot(effect("block*regroup",USex1sa))
  101. plot(effect("stim*block*smoker",USex1sm))
  102. plot(effect("regroup",USex1a))
  103. # _ _ _____ _
  104. # | | | |/ ____| | |
  105. # | | | | (___ _____ ___ __ ___ ___| |_ __ _ _ __ ___ _ _
  106. # | | | |\___ \ / _ \ \/ / '_ \ / _ \/ __| __/ _` | '_ \ / __| | | |
  107. # | |__| |____) | | __/> <| |_) | __/ (__| || (_| | | | | (__| |_| |
  108. # \____/|_____/ \___/_/\_\ .__/ \___|\___|\__\__,_|_| |_|\___|\__, |
  109. # | | |__ \ | | __/ |
  110. # __| | __ _ _ _ ) ||_| |___/
  111. # / _` |/ _` | | | | / /
  112. # | (_| | (_| | |_| | / /_
  113. # \__,_|\__,_|\__, | |____|
  114. # __/ |
  115. # |___/
  116. d_ex<-read.csv(file='/Madeleine/Online Studie/analysis/R/D2USexponstudR.csv',head=TRUE, sep=",")
  117. #all four groups (1=non-smoker, 2=smoker with 6h break after ACQ, 3=smoker with cigarette after ACQ, 4=smoker with cigarette before ACQ)
  118. d_ex$group<-as.factor(d_ex$group)
  119. #all 4 groups, see above, but smokers from group 2 that did smoke during the 6h smoke break were regrouped into group 3
  120. d_ex$regroup<-as.factor(d_ex$regroup)
  121. #stimulus: CSP=CS+, GS2-GS9=generalized stimuli, CSM=CS-
  122. d_ex$stim<-as.factor(d_ex$stim)
  123. #block 1 and 2, one block consists of 1 trial per stimulus, so 10 trials in total
  124. d_ex$block<-as.factor(d_ex$block)
  125. # CS-US contingency, 1=not aware of correct CS-US contingency, 2=aware of correct CS-US contingency
  126. d_ex$aware<-as.factor(d_ex$aware)
  127. #smoker(group2,3 and 4)= 1, non-smoker(group1) = 0
  128. d_ex$smoker<-as.factor(d_ex$smoker)
  129. #subject number
  130. d_ex$sub<-d_ex$sub
  131. #########################################
  132. #models
  133. #########################################
  134. #models for either smoker vs non-smoker or regroup
  135. D2onstud<-lmer(rating~(1|sub)+stim*block*smoker,data = d_ex,control = lmerControl(optimizer = "bobyqa"))
  136. D2bonstud<-lmer(rating~(1|sub)+stim*block*regroup,data = d_ex,control = lmerControl(optimizer = "bobyqa"))
  137. Anova(D2onstud, type="3", test="F" )
  138. plot(effect("stim",D2onstud))
  139. plot(effect("regroup",D2bonstud))
  140. plot(effect("block",D2onstud))
  141. plot(effect("stim*regroup",D2bonstud))
  142. plot(effect("stim*block*smoker",D2onstud))
  143. plot(effect("stim*smoker",D2onstud))
  144. #posthoc
  145. emmeans(D2onstud, list(pairwise ~ stim:block:smoker), adjust = "none")