# Numerical model for predicting effects of timing imprecision on impulsivity (delay discounting) library(ggplot2) library(cowplot) # Valuer function gives you v(t) given k and t valuer=function(t, k=0.54){return(1/(1 + k*t))} # Function to run the model # Considers two individuals, a precise one and an imprecise one # ap and ai are the alpha values of the precise and imprecise individual # bp and bi are the beta values of the precise and imprecise # ts and tl are the delays to the long and short reward model=function(ap, ai, bp, bi, ts=3, tl=8, k=0.54, n=100000) { precise.estimates.short=rnorm(n, mean=ts, sd=sqrt(ap+bp*ts)) precise.estimates.short[precise.estimates.short<0]=0 imprecise.estimates.short=rnorm(n, mean=ts, sd=sqrt(ai+bi*ts)) imprecise.estimates.short[imprecise.estimates.short<0]=0 precise.estimates.long=rnorm(n, mean=tl, sd=sqrt(ap+bp*tl)) precise.estimates.long[precise.estimates.long<0]=0 imprecise.estimates.long=rnorm(n, mean=tl, sd =sqrt(ai+bi*tl)) imprecise.estimates.long[imprecise.estimates.long<0]=0 type=c("precise", "imprecise") values.short=c(mean(valuer(precise.estimates.short, k=k)), mean(valuer(imprecise.estimates.short, k=k))) values.long=c(mean(valuer(precise.estimates.long, k=k)), mean(valuer(imprecise.estimates.long, k=k))) ratio=values.short/values.long output=data.frame(type, values.short, values.long, ratio) return(ratio[2]/ratio[1]) } ###Figure 2##### # The outcome variable is the relative valuation of the smaller, sooner reward # i.e. more than 1 indicates that the imprecise individual relatively values the smaller, sooner reward more # Less than 1, the opposite ap=0 bp=1 ais=c(0, 1, 2) bis=seq(0.5, 2, by=0.05) answer=NULL current.ai=NULL current.bi=NULL for (i in 1:length(ais)){ for (j in 1: length(bis)){ answer[(i-1)*length(bis) +j] = model(ap=ap, bp=bp, ai=ais[i], bi=bis[j], k=.54) current.ai[(i-1)*length(bis) +j]=ais[i] current.bi[(i-1)*length(bis) +j]=bis[j] } } d1=data.frame(current.ai, current.bi, answer) # Plot figure # Figure 2A f1=ggplot(d1, aes(x=current.bi, y=answer, group=factor(current.ai))) + geom_line(aes(linetype=factor(current.ai))) + theme_bw() + geom_hline(yintercept = 1, linetype="twodash") + ylab("Relative valuation of SS") + xlab(expression(beta)) + labs(linetype=expression(alpha)) + coord_cartesian(ylim=c(0.9, 1.1)) + theme(axis.title = element_text(size = 14), legend.text = element_text(size=14), legend.title = element_text(size=14), axis.text=element_text(size=12)) # Lower panels: Make a plot of a hyperbolic discount function delay=seq(1, 10, by=.1) value=1/(1+0.54*delay) p=data.frame(delay, value) fb=ggplot(p, aes(x=delay, y=value)) + geom_line() + xlab("Delay") + ylab("Value") + theme(axis.text = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.title=element_text(size =14)) fb bottom.row=plot_grid(fb, fb, labels=c("B", "C")) fig2=plot_grid(f1, bottom.row, ncol=1, labels=c("A", "")) fig2 # Make figure 2 (normal curves are added on in Powerpoint) win.metafile("figure2.wmf", width=8, height=8) fig2 dev.off() # Try higher and lower values of k # Higher value answer=NULL current.ai=NULL current.bi=NULL for (i in 1:length(ais)){ for (j in 1: length(bis)){ answer[(i-1)*length(bis) +j] = model(ap=ap, bp=bp, ai=ais[i], bi=bis[j], k=1) current.ai[(i-1)*length(bis) +j]=ais[i] current.bi[(i-1)*length(bis) +j]=bis[j] } } d2=data.frame(current.ai, current.bi, answer) # Figure 2A with higher value of k f2Abis=ggplot(d2, aes(x=current.bi, y=answer, group=factor(current.ai))) + geom_line(aes(linetype=factor(current.ai))) + theme_bw() + geom_hline(yintercept = 1, linetype="twodash") + ylab("Relative valuation of SS") + xlab(expression(beta)) + coord_cartesian(ylim=c(0.9, 1.1)) + guides(linetype=FALSE) f2Abis # Lower value of k answer=NULL current.ai=NULL current.bi=NULL for (i in 1:length(ais)){ for (j in 1: length(bis)){ answer[(i-1)*length(bis) +j] = model(ap=ap, bp=bp, ai=ais[i], bi=bis[j], k=0.1) current.ai[(i-1)*length(bis) +j]=ais[i] current.bi[(i-1)*length(bis) +j]=bis[j] } } d3=data.frame(current.ai, current.bi, answer) f2Abis2=ggplot(d3, aes(x=current.bi, y=answer, group=factor(current.ai))) + geom_line(aes(linetype=factor(current.ai))) + theme_bw() + geom_hline(yintercept = 1, linetype="twodash") + ylab("Relative valuation of SS") + xlab(expression(beta)) + coord_cartesian(ylim=c(0.9, 1.1)) + guides(linetype=FALSE) f2Abis2