TEXT 39
Code.txt Guest on 24th August 2020 05:50:48 AM
  1. ## R version information
  2. # R.Version()
  3. # $platform
  4. # [1] "x86_64-w64-mingw32"
  5. # $version.string
  6. # [1] "R version 3.2.1 (2015-06-18)"
  7. # $nickname
  8. # [1] "World-Famous Astronaut"
  9.  
  10. ## NOTE: Split_df_v4_2/ filelist[[2]] is the main experiment reported in the manuscript.
  11. ## Split_df_v4_1 and _3 / filelist[[1]] & [[3]] are the experiments reported in the appendix (experiment 2 and 3, respectively)
  12. ## THESE NAMES ARE NOW CHANGED to split1.txt ...
  13.  
  14. ## load packages
  15. require(plyr)
  16. require(ggplot2)
  17. require(grid)
  18. require(lme4)
  19. require(multcomp)
  20. require(reshape2)
  21.  
  22. ## workspace (after running the following commented code)
  23. load("/DDM_JDM_v2.RData")
  24.  
  25. #
  26. # ## read in raw data
  27. # setwd("/dataframes_v4/")
  28. # filenames <- list.files(pattern="split_df_v4*.*txt")
  29. # filelist <- as.list(lapply(filenames, read.table, header=T))
  30. # names(filelist) <- c("S1","S2","S3")
  31. ## note: S1 = reduced nutrition, S2=full info (main paper), S3= reduced, no percentages
  32. #
  33. # ## data type conversion
  34. # for (x in seq_along(filelist)) {
  35. #   filelist[[x]]$GDA_indi = as.factor(filelist[[x]]$GDA_indi) ## 1=GDA, 2=TL
  36. #   filelist[[x]]$RT_choice = as.numeric(as.character(filelist[[x]]$RT_choice))
  37. #   filelist[[x]]$rating_bins = as.factor(filelist[[x]]$rating_bins)
  38. #   filelist[[x]]$rating_diff = as.numeric(as.character(filelist[[x]]$rating_diff))
  39. #   }
  40. #
  41. # ## join dataframe for data cleaning
  42. # joined_df <- do.call("rbind", filelist)
  43.  
  44. # ######################################
  45. # ######## data cleaning ###############
  46. # ######################################
  47. #
  48. # ## first remove all RTs > 30 -- ztree errors, mean and SD not robust enough for really long outliers
  49. # joined_df<-joined_df[joined_df$RT_choice<=30,]
  50. #
  51. # ## subject - specific RT, remove 2 SD above subject-specific mean RT
  52. # temp<-split(joined_df, f=joined_df$ID)
  53. # ## create empty list of dataframes
  54. # mittel<-lapply(1:130,function(i) data.frame())
  55. # stabw<-lapply(1:130,function(i) data.frame())
  56. #
  57. # for (x in seq_along(temp)){
  58. #   mittel[[x]]<-mean(temp[[x]]$RT_choice)
  59. #   stabw[[x]]<-sd(temp[[x]]$RT_choice)
  60. #   ## then remove RT more than 2 SD above individual mean
  61. #   temp[[x]]<-temp[[x]][temp[[x]]$RT_choice<= mittel[[x]]+2*stabw[[x]],]
  62. # }
  63. #
  64. # ## how many participants excluded in main experiment (sbj- IDs 44-87 are from main experiment)
  65. # excl <- sapply(temp, nrow)
  66. # 350-mean(excl[44:87])
  67. # sd(excl[44:87])
  68. # range(excl[44:87])
  69. #
  70. # ## bind again
  71. # joined_df <- do.call("rbind", temp)
  72. # filelist <-split(joined_df, f=joined_df$session)
  73. #
  74. # ## save (for diffusion modeling/ drift as a function (here session 2 only))
  75. # write.table(filelist[[2]], "/one_drift_theta/split_df_2_v4.csv", sep=";", quote=FALSE, row.names=FALSE)
  76. # write.table(filelist[[2]][,c(8,10,11)], "/one_drift_theta/split_df_2_v4_fastdm.txt", sep="\t", quote=FALSE, row.names=FALSE)
  77. # write.table(filelist[[2]][,c(8,10,11,12)], "/one_drift_theta/split_df_2_v4_fastdm_ratingbins.txt", sep="\t", quote=FALSE, row.names=FALSE)
  78. #
  79. # ###########################################################################################
  80. # ###############################log-mixed model ###########################################
  81. # ###########################################################################################
  82.  
  83.  
  84. ## model "Label"
  85. logmixed1<-lapply(filelist, function(z) glmer(h_choice ~ GDA_indi  + (1 + GDA_indi|ID), data=z, family="binomial"))
  86. ## stats
  87. stats1<-lapply(logmixed1, function(z) summary(z))
  88.  
  89.  
  90. ## model Liking + Label
  91. ## interaction with ratings
  92. logmixed2<-lapply(filelist, function(z) glmer(h_choice ~ GDA_indi + rating_diff  + (1 + GDA_indi+rating_diff |ID), data=z, family="binomial"))
  93. ## stats
  94. stats2<-lapply(logmixed2, function(z) summary(z))
  95.  
  96.  
  97.  ## model Liking * Label
  98. logmixed3<-lapply(filelist, function(z) glmer(h_choice ~  rating_diff*GDA_indi + (1 + GDA_indi*rating_diff | ID),data=z ,family="binomial"))
  99. ## stats
  100. stats3<-lapply(logmixed3, function(z) summary(z))
  101.  
  102. ##############################
  103. ### Reaction time ############
  104. ##############################
  105. ## log transform RT due to skewness in binary choice paradigms
  106. for (b in seq_along(filelist)) {
  107.   filelist[[b]]$RT_choice_log =  log(filelist[[b]]$RT_choice)
  108. }
  109.  
  110. RT.model <- lapply(filelist, function (z) lmer(RT_choice_log ~ GDA_indi  + (1 + GDA_indi|ID), data =z, REML=FALSE))
  111. ## stats
  112.  
  113. ### summary RT
  114. RT_sum <-ddply(filelist[[2]], "GDA_indi", summarize, RT.mean=mean(RT_choice_log),
  115.            RT.sd=sd(RT_choice_log))
  116.  
  117. ############################################
  118. ######## plots #############################
  119. ############################################
  120.  
  121. ## percentage of healthy choices session 2 (main experiment)
  122. ### FIGURE 3
  123.  
  124. ## split into individual dataframes
  125. temp<-split(filelist[[2]], paste(filelist[[2]]$ID, filelist[[2]]$GDA_indi))
  126. ## empty vector for looping
  127. vec<-vector("list", length(temp))
  128.  
  129. ## sum healthy choices
  130. for (x in seq_along(temp)){
  131.   vec[[x]] <- sum(temp[[x]]$h_choice)/ nrow(temp[[x]])
  132. }
  133.  
  134. ## bind into dataframe
  135. temp_binded<-as.data.frame(do.call("rbind", vec))
  136. ## GDA_indi
  137. temp_binded$GDA_indi<-c(1,2)
  138. ## add column names
  139. colnames(temp_binded)<-c("Percent_healthy_choice", "GDA_indi")
  140. temp_binded$GDA_indi<-as.factor(temp_binded$GDA_indi)
  141. ### summary stats
  142. fig3_plo<-ddply(temp_binded, "GDA_indi", summarize, Percent.mean=mean(Percent_healthy_choice),
  143.            Percent.sd=sd(Percent_healthy_choice),
  144.            Lower=Percent.mean - Percent.sd/sqrt(length(unique(filelist[[2]]$ID))),
  145.            Upper=Percent.mean + Percent.sd/sqrt(length(unique(filelist[[2]]$ID))))
  146.  
  147. ## plot Figure 3
  148. fig3 <- ggplot(fig3_plo, aes(x=GDA_indi, y=Percent.mean, group=factor(GDA_indi),shape=factor(GDA_indi), colour=factor(GDA_indi)))+
  149.   geom_bar(fill=c("gray80", "gray20"),stat="identity",  colour="black", size=.6, width=0.5) +
  150.   geom_errorbar(aes(ymin=Lower, ymax=Upper), width=0.1, colour="black")+
  151.   theme_classic(base_size = 14, base_family ="serif")+theme(panel.grid = element_blank()) +
  152.   scale_colour_manual(values=c("gray80", "gray20"), guide=FALSE)+
  153.   coord_cartesian(ylim=c(0.0, 1))+  scale_x_discrete(name="Label", labels=c("GDA", "TL"))+
  154.   scale_y_continuous(name="Probability of healthy choices", expand=c(0,0))+
  155.   theme(legend.key.width =unit(1.5, "cm"), axis.ticks.x = element_blank()) +
  156.   geom_text(x = 1.5, y = 0.63, label = "**", size=10) ## add aterisks
  157.  
  158. fig3
  159. ggsave("E:/Figure3.tiff", width=6, height=8, dpi=300)
  160. ggsave("E:/Figure3.eps", family="serif", width=6, height=8, dpi=300)
  161.  
  162.  
  163. ### Figure 4
  164. ### Interaction with liking of taste, using empirical vs. predicted probabilities
  165. ## predict with random effect (ID) # IDs:  44-87
  166. ## note that liking ratings were binned for illustration purposes only
  167.  
  168. logmixed_ranks_fig4<-glmer(h_choice ~ rating_bins + GDA_indi + rating_bins*GDA_indi + (1|ID),data=filelist[[2]],family="binomial")
  169.  
  170. new1<-data.frame(GDA_indi=factor(c(1,1,1,1,1,1,1,2,2,2,2,2,2,2)), rating_bins=factor(1:7), ID=factor(rep(c(44:87), each=14)), h_choice=0)
  171. new1$prob<-predict(logmixed_ranks_fig4,newdata=new1, re.form=NULL, type="response")
  172. ## predicted data
  173. new1<-ddply(new1, .(GDA_indi, rating_bins), summarize, prob.mean=mean(prob))
  174.  
  175. ## empirical data
  176. h_choice<-ddply(filelist[[2]], .(rating_bins,GDA_indi), summarise, h_choice_sum = sum(h_choice))
  177. occurences<-data.frame(table(filelist[[2]]$GDA_indi, filelist[[2]]$rating_bins))
  178. colnames(occurences)<-(c("GDA_indi", "rating_bins", "Number_of_choices"))
  179. emp_h_choice <- join(h_choice, occurences, by=c("GDA_indi", "rating_bins"),type="left")
  180.  
  181. ## percent healthy choices
  182. emp_h_choice$prob.mean<- emp_h_choice$h_choice_sum/emp_h_choice$Number_of_choices
  183.  
  184. ## GDA healthy ## empirical == 1, predicted ==2
  185. healthy_choice_GDA_emppred<-rbind.fill(emp_h_choice,new1) ## rbind fill from plyr package when uneven columns/differnct columns
  186. healthy_choice_GDA_emppred$emppred<-as.factor(rep(c(1:2), each = 14 )) ## 7 rating bins * 2 label types
  187.  
  188. ## PLOT Figure 4
  189. ### to connect lines of only predicted, but not empirical: subset (see below)
  190. fig4 <- ggplot(healthy_choice_GDA_emppred, aes(x=rating_bins, y=prob.mean, group=emppred, color=GDA_indi)) +
  191.   geom_point(aes(shape=emppred),size=4)  + scale_shape_manual(values=c(4,1), labels = c("Empirical", "predicted"), name="") +
  192.   scale_x_discrete(breaks = c(1:7)) + theme_bw(base_size = 14) +
  193.   xlab("Taste unhealthy > healthy                                       Taste healthy > unhealthy") + ylab("Probability of healthy choice")+
  194.   scale_colour_manual(values=c("gray60", "gray20"),name="Label",labels=c("GDA", "TL")) +
  195.   geom_line(data=subset(healthy_choice_GDA_emppred,emppred==2), aes(factor(rating_bins), prob.mean, group=GDA_indi)) +
  196.   ggtitle("Empirical vs. predicted probability of healthy choices") +
  197.   theme_classic(base_size = 14, base_family ="serif")+theme(panel.grid = element_blank())
  198. fig4
  199.  
  200. ggsave("E:/Figure4.tiff",width=8, height=6, dpi=300)
  201. ggsave("E:/Figure4.eps", family ="serif", width=8, height=6, dpi=300)
  202.  
  203.  
  204.  
  205. #######################################################
  206. ############# test for weights effect (APPENDIX) ######
  207. #######################################################
  208.  
  209. ## other factors might contribute to the choice and so in theory the labels could make the subjects
  210. ## put more weight on taste and health, while lowering the weight on other unnamed factors.  
  211. ## One simple way to address this is to just run a logit with only taste predicting choice, once for each label type.  
  212. ## Then comparing the coefficients -- If they go down in the TL condition, we can say that the weight on taste
  213. ## has indeed gone down in the TL condition.
  214. ## Then  the same thing with only health and see if  indeed that coefficient goes up in TL.
  215.  
  216. ## Recode so that we can use choice_left as DV
  217. filelist[[2]]$choice_left[filelist[[2]]$h_left == 1 & filelist[[2]]$h_choice == 1]<-1
  218. filelist[[2]]$choice_left[filelist[[2]]$h_left == 1 & filelist[[2]]$h_choice == 0]<-0
  219. filelist[[2]]$choice_left[filelist[[2]]$h_left == 0 & filelist[[2]]$h_choice == 1]<-0
  220. filelist[[2]]$choice_left[filelist[[2]]$h_left == 0 & filelist[[2]]$h_choice == 0]<-1
  221. filelist[[2]]$choice_left<-as.factor(filelist[[2]]$choice_left)
  222. filelist[[2]]$h_left<-as.factor(filelist[[2]]$h_left)
  223.  
  224. ## now adjust rating differences to left right vs. healthy unhealthy
  225. filelist[[2]]$taste_left <- ifelse(filelist[[2]]$h_left == 1, filelist[[2]]$rating_h , filelist[[2]]$rating_u)
  226. filelist[[2]]$taste_right <- ifelse(filelist[[2]]$h_left == 0, filelist[[2]]$rating_h , filelist[[2]]$rating_u)
  227. filelist[[2]]$taste_left_right <- filelist[[2]]$taste_left - filelist[[2]]$taste_right
  228.  
  229. ## split into TL and GDA, first file is GDA, second is TL
  230. onlyTLGDA <- (split(filelist[[2]], filelist[[2]]$GDA_indi))
  231.  
  232. summary(rev1<- glmer(choice_left ~ taste_left_right  + (1 | ID),data=onlyTLGDA[[1]],family="binomial"))
  233. summary(rev2<- glmer(choice_left ~ taste_left_right  + (1 | ID),data=onlyTLGDA[[2]],family="binomial"))
  234.  
  235. summary(rev3<- glmer(choice_left ~ h_left  + (1 | ID),data=onlyTLGDA[[1]],family="binomial"))
  236. summary(rev4<- glmer(choice_left ~ h_left  + (1 | ID),data=onlyTLGDA[[2]],family="binomial"))
  237.  
  238.  
  239. ###############################################################
  240. ################# fast dm #####################################
  241. ###############################################################
  242.  
  243. ## you need individual dataframes per person for fast-dm
  244. #
  245. # overall<-filelist[[2]][,c("ID", "GDA_indi", "h_choice", "RT_choice" )]
  246. # ID_overall<-split(overall, overall$ID)
  247. # ID_overall<-lapply(ID_overall, function(x)  x[c(2,3,4)])
  248. # ## save as dat file
  249. # sapply(names(ID_overall), function (x) write.table(ID_overall[[x]], file=paste("E:/ID_overall_", x, ".dat", sep=""),
  250. #                                                    sep="\t", col.names = F, row.names = F, quote=F))
  251. #
  252. # ## now you can run the fast-dm code (separate code file used in fast dm, code is provided)
  253.  
  254. ## then read in fast-dm output
  255. setwd("/fast_dm/DDM_v4_overall_drift_and_startingpoint_2_singlesubject/")
  256. results_drift <-read.table("overall_cs_drift.log", header=T)
  257. results_drift_z <- read.table("overall_cs_drift_and_startingpoint.log", header=T)
  258. results_drift_t0 <-read.table("overall_cs_drift_and_t0.log", header=T)
  259. results_drift_z_t0 <-read.table("overall_cs_drift_and_startingpoint_and_t0.log", header=T)
  260.  
  261. ## convert columns to numeric
  262. results_drift<-as.data.frame(apply(results_drift, 2, function(x) as.numeric(as.character(x))))
  263. results_drift_z<-as.data.frame(apply(results_drift_z, 2, function(x) as.numeric(as.character(x))))
  264. results_drift_t0<-as.data.frame(apply(results_drift_t0, 2, function(x) as.numeric(as.character(x))))
  265. results_drift_z_t0<-as.data.frame(apply(results_drift_z_t0, 2, function(x) as.numeric(as.character(x))))
  266.  
  267. ## model drift t test
  268. t.test(results_drift$v_1,results_drift$v_2, paired=TRUE)
  269. ## means and sd
  270. print(mean<-apply(results_drift, 2, function(x) mean(x)), na.rm=TRUE)
  271. sd<- as.data.frame(apply(results_drift, 2, function(x) sd(x)))
  272. print(sem<- apply(sd, 2, function(x) sqrt(x)/43))
  273.  
  274. ## model drift and starting point t test
  275. t.test(results_drift_z$v_1,results_drift_z$v_2, paired=TRUE)
  276. t.test(results_drift_z$zr_1,results_drift_z$zr_2, paired=TRUE)
  277. ## means and sd
  278. print(mean<-apply(results_drift_z, 2, function(x) mean(x)))
  279. sd<- as.data.frame(apply(results_drift_z, 2, function(x) sd(x)))
  280. print(sem<- apply(sd, 2, function(x) sqrt(x)/43))
  281.  
  282. ## model drift and nondecision time t test
  283. t.test(results_drift_t0$v_1,results_drift_t0$v_2, paired=TRUE)
  284. t.test(results_drift_t0$t0_1,results_drift_t0$t0_2, paired=TRUE)
  285. ## means and sd
  286. print(mean<-apply(results_drift_t0, 2, function(x) mean(x)))
  287. sd<- as.data.frame(apply(results_drift_t0, 2, function(x) sd(x)))
  288. print(sem<- apply(sd, 2, function(x) sqrt(x)/43))
  289.  
  290. ## model drift and nondecision and startingpoint t test
  291. t.test(results_drift_z_t0$v_1,results_drift_z_t0$v_2, paired=TRUE)
  292. t.test(results_drift_z_t0$zr_1,results_drift_z_t0$zr_2, paired=TRUE)
  293. t.test(results_drift_z_t0$t0_1,results_drift_z_t0$t0_2, paired=TRUE)
  294. ## means and sd
  295. print(mean<-apply(results_drift_z_t0, 2, function(x) mean(x)))
  296. sd<- as.data.frame(apply(results_drift_z_t0, 2, function(x) sd(x)))
  297. print(sem<- apply(sd, 2, function(x) sqrt(x)/43))
  298.  
  299.  
  300. ## melt for plotting
  301. results_drift_z_t0 <- results_drift_z_t0[, c(8, 5, 3, 2, 9, 6)] ## only columns of interest
  302. results_drift_z_t0_melt<-melt(results_drift_z_t0)
  303. colnames(results_drift_z_t0_melt)<-c("labeling", "Estimate")
  304.  
  305. ## for drift and starting point and t0
  306. results_drift_z_t0_melt$Label[results_drift_z_t0_melt$labeling=="v_1"]<-"GDA"
  307. results_drift_z_t0_melt$Label[results_drift_z_t0_melt$labeling=="zr_1"]<-"GDA"
  308. results_drift_z_t0_melt$Label[results_drift_z_t0_melt$labeling=="t0_1"]<-"GDA"
  309. results_drift_z_t0_melt$Label[results_drift_z_t0_melt$labeling=="v_2"]<-"TL"
  310. results_drift_z_t0_melt$Label[results_drift_z_t0_melt$labeling=="zr_2"]<-"TL"
  311. results_drift_z_t0_melt$Label[results_drift_z_t0_melt$labeling=="t0_2"]<-"TL"
  312.  
  313. results_drift_z_t0_melt$variable<-rep(c("Drift rate", "Starting point", "Non-decision time"), each=nrow(results_drift_z_t0_melt)/3)
  314. results_drift_z_t0_melt$Label <- as.factor(results_drift_z_t0_melt$Label)
  315.  
  316.  
  317. ## zr vs v vs. t0 plot summary
  318. ## FIGURE 5
  319. results_drift_z_t0_melt$Label<-as.factor(results_drift_z_t0_melt$Label)
  320. overall_plot_Fastdm<-ddply(results_drift_z_t0_melt, .(Label, variable), summarize, Estimate.mean=mean(Estimate), Estimate.sd=sd(Estimate),
  321.                     Lower=Estimate.mean - Estimate.sd/sqrt(nrow(results_drift_z_t0_melt)/2),
  322.                     Upper=Estimate.mean + Estimate.sd/sqrt(nrow(results_drift_z_t0_melt)/2))
  323.  
  324. fig5 <- ggplot(overall_plot_Fastdm, aes(y=Estimate.mean, x=Label, fill=Label)) + facet_grid(.~variable)+
  325.   geom_bar(stat="identity",  colour="black", size=.6, width=0.5) +
  326.   scale_fill_manual(values=rep(c("gray80", "gray20"), times = 3)) + labs(fill="")+
  327.   geom_errorbar(aes(ymin=Lower, ymax=Upper), width=0.1)+
  328.   ylab("Mean Parameter Estimates")+ facet_grid(.~variable)+ scale_x_discrete(labels=c("GDA", "TL")) +
  329.   geom_hline(yintercept=0, colour="gray50") + ## adds horizontal line at zero
  330.   theme_classic(base_size = 20, base_family ="serif")+theme(panel.grid = element_blank())+
  331.   theme(axis.title.x=element_blank(), axis.text.x=element_blank(), axis.ticks.x=element_blank())
  332.  
  333. ## add aterisk only to first facet of plot
  334. fig5 + geom_text(data=data.frame(x=1.5, y=0.15, label=c("*", "",""),
  335.                         variable=c("Drift rate", "Non-decision time", "Starting point")),
  336.                         aes(x,y,label=label), inherit.aes=FALSE, size=10)
  337.  
  338. ggsave("E:/Figure5.tiff",width=8, height=6, dpi=300)
  339. ggsave("E:/Figure5.eps", family ="serif", width=8, height=6, dpi=300)
  340.  
  341.  
  342.  
  343. ##################################################################
  344. ################# Drift as function ##############################
  345. ##################################################################
  346. #
  347. # ### function
  348. # rm(list=ls())
  349. # require(RWiener)
  350. # label <- c(1,2) ## 1= GDA; 2= TL
  351. # # pull in the data files
  352. #
  353. # for (y in 1:length(label)){
  354. # labels = label[y]
  355. # data <- read.table('/one_drift_theta/split_df_2_v4.csv', sep=";",  header = TRUE)
  356. #
  357. # ## find min RT (for modeling later)
  358. # minRT <- min(data$RT_choice)
  359. #
  360. # # Personal deviance function
  361. # theta_dev <- function(startPar)
  362. # {
  363. #   alpha <- startPar[1]
  364. #   tau <- startPar[2]
  365. #   beta <- startPar[3]
  366. #   thetaA <- startPar[4]
  367. #   thetaB <- startPar[5]
  368. #  
  369. #   delta <- thetaA+(thetaB*rateDiff)
  370. #  
  371. #   p1 = rep(0, length(dat$RT))
  372. #   p0 = rep(0, length(dat$RT))
  373. #  
  374. #   for (t in 1:length(dat$RT))
  375. #   {
  376. #     RTt <- dat$RT[t]
  377. #     deltat <- delta[t]
  378. #    
  379. #     p1[t] <- dwiener(RTt, alpha, tau, beta, deltat, resp = "upper", give_log = FALSE)
  380. #     p0[t] <- dwiener(RTt, alpha, tau, beta, deltat, resp = "lower", give_log = FALSE)
  381. #   }
  382. #   deviance <- -2*sum((log(p1)*choice)+(log(p0)*(1-choice)))
  383. #   if (is.finite(deviance) == FALSE)
  384. #   {
  385. #     deviance <- 100000
  386. #   }
  387. #   return(deviance)
  388. # }
  389. #
  390. #
  391. # ## trinary rating difference (coarse)
  392. # data$rating_3<-NA
  393. # data$rating_3[data$rating_diff>=-10]<- -1
  394. # data$rating_3[data$rating_diff>=-3]<-0
  395. # data$rating_3[data$rating_diff> 3]<- 1
  396. #
  397. # # Only take the variables we need right now
  398. # data <- data.frame(subjectNum = data$ID, RT = data$RT, Acc = data$h_choice, rateDiff = data$rating_3, GDA_indi=data$GDA_indi)
  399. #
  400. # # create data list to put all the info into
  401. # allProbPar <- data.frame(alpha = numeric(), tau = numeric(), beta = numeric(), a = numeric(), b = numeric(), deviance = numeric(), convergence = numeric())
  402. #
  403. # # First, rename the subject numbers
  404. # currentSub = 1
  405. # for (i in 1:max(data$subjectNum))
  406. # {
  407. #   if (length(data[data$subjectNum == i, ]$subjectNum > 0))
  408. #   {
  409. #     data[data$subjectNum == i, ]$subjectNum <- rep(currentSub, length(data[data$subjectNum == i, "subjectNum"]))
  410. #     currentSub = currentSub+1
  411. #   }
  412. # }
  413. #
  414. # ## split data according to TL and GDA
  415. # data_part<-split(data, data$GDA_indi)
  416. #
  417. # if(labels == 1)
  418. # {data_optim<-data_part[[1]] }
  419. # else if(labels ==2)
  420. #   {data_optim<-data_part[[2]]}  
  421. #
  422. #
  423. # # Now do all calculations on the individual subject level
  424. #
  425. # for (i in 1:max(data_optim$subjectNum))
  426. # {
  427. #   # First pull only the subject you're working on
  428. #   dat <- data_optim[data_optim$subjectNum == i, ]
  429. #  
  430. #   # Now RT cut the data (necessary for estimation purposes, dep. on empirical data)
  431. #  
  432. #   dat <- dat[dat$RT > minRT, ]
  433. #  
  434. #   rateDiff <- dat$rateDiff
  435. #   choice <- dat$Acc    
  436. #  
  437. #   dat$rateDiff <- NULL
  438. #   dat$subjectNum <- NULL
  439. #   dat$GDA_indi <- NULL
  440. #   dat[dat$Acc == 1, 'Acc'] <- 'upper'
  441. #   dat[dat$Acc ==0, 'Acc'] <- 'lower'
  442. #  
  443. #   startPar <- c(sample(seq(.5,4.5,.1), 1), sample(seq(.1, .4, .01), 1), sample(seq(.2, .8, .01), 1), sample(seq(0, 1, .01), 1), sample(seq(0, 1, .01), 1))
  444. #  
  445. #   # Get the paramaters!
  446. #   probPar <- optim(startPar, theta_dev, method = "L-BFGS-B", lower = c(.5, .2, .4, -1, -1), upper = c(4.5, 1.5, 1, 1, 1), control = list(maxit = 5000))
  447. #
  448. #   allProbPar[nrow(allProbPar)+1, ] <- c(probPar$par, probPar$value, probPar$convergence)
  449. # }
  450. # if(labels == 1)
  451. # {write.table(allProbPar, '/one_drift_theta/individualTheta_choiceProb_GDA_L-BFGS.csv', sep = ",")}
  452. # else if(labels ==2)
  453. # {write.table(allProbPar, '/one_drift_theta/individualTheta_choiceProb_TL_L-BFGS.csv', sep = ",")}  
  454. # }
  455.  
  456.  
  457. ## Read OUTPUt form drift-as-function (DAS) after estimation loop
  458. GDA_das_fun <- read.table('/one_drift_theta/individualTheta_choiceProb_GDA_L-BFGS.csv', sep=",", header=TRUE)
  459. GDA_das_fun <- GDA_das_fun[, c("a", "b")]
  460. TL_das_fun <- read.table('/one_drift_theta/individualTheta_choiceProb_TL_L-BFGS.csv', sep=",", header=TRUE)
  461. TL_das_fun <- TL_das_fun[, c("a", "b")]
  462.  
  463. GDA_TL_das_fun <- cbind(GDA_das_fun, TL_das_fun)
  464. colnames(GDA_TL_das_fun) <- c("GDA_a", "GDA_b", "TL_a", "TL_b")
  465.  
  466. ## test
  467. t.test(GDA_TL_das_fun$GDA_a, GDA_TL_das_fun$TL_a, paired=TRUE)
  468. t.test(GDA_TL_das_fun$GDA_b, GDA_TL_das_fun$TL_b, paired=TRUE)
  469.  
  470. ## means and sd
  471. print(mean<-apply(GDA_TL_das_fun, 2, function(x) mean(x)))
  472. sd<- as.data.frame(apply(GDA_TL_das_fun, 2, function(x) sd(x)))
  473. print(sem<- apply(sd, 2, function(x) sqrt(x)/43))
  474.  
  475.  
  476.  
  477. ## Plotting FIGURE 6
  478. data_part1 <- melt(GDA_TL_das_fun)
  479.  
  480. ## add labels and which coefficient
  481. data_part1$GDA_indi<-rep(c("GDA", "TL"), each=88)
  482. data_part1$coeff<-rep(c("a", "b"), each=44)
  483.  
  484. fig6_plot<-ddply(data_part1, c("GDA_indi", "coeff"), summarize, coeff.mean=mean(value, na.rm=T),coeff.sd=sd(value, na.rm=T),
  485.             Lower=coeff.mean - coeff.sd/sqrt(44),
  486.             Upper=coeff.mean + coeff.sd/sqrt(44))
  487.  
  488. fig6 <- ggplot(data=fig6_plot, aes(x=coeff, y=coeff.mean, fill = GDA_indi, group=GDA_indi)) +
  489.   geom_bar(stat="identity", position=position_dodge(), colour="black", size=.6, width=0.5) +
  490.   geom_errorbar(aes(ymin=Lower, ymax=Upper), width=0.1, colour="black", position=position_dodge(width=0.5))+
  491.   scale_fill_manual(values=c("gray80", "gray20", "gray80", "gray20")) + labs(fill="")+
  492.   theme_classic(base_size = 14, base_family ="serif")+theme(panel.grid = element_blank()) +
  493.   coord_cartesian(ylim=c(-0.1, 1.1))+ scale_x_discrete(name="Coefficient", labels=c("health sensitivity", "taste weight"))+
  494.   scale_y_continuous(name="Coefficient mean", expand=c(0,0)) + geom_hline(yintercept=0 )+
  495.   geom_vline(xintercept=0.4) + theme(axis.line=element_blank()) +
  496.   theme(legend.key.width =unit(1.5, "cm"), axis.ticks.x = element_blank())
  497. ## add aterisks
  498. fig6 + annotate("text", x = c(1,2), y = c(0.2, 0.95), label = c("*", "*"), size=10)
  499.  
  500. ggsave("E:/Figure6.tiff",width=8, height=6, dpi=300)
  501. ggsave("E:/Figure6.eps", family ="serif", width=8, height=6, dpi=300)

Paste is for source code and general debugging text.

Login or Register to edit, delete and keep track of your pastes and more.

Raw Paste

Login or Register to edit or fork this paste. It's free.