This notebook contains simulations to examine different options for modeling the design of the NARPS study.

The goal of the study is to present subjects with one of two different gamble distributions:

The statistical model for the first level responses (after Tom et al.) will be:

\[y = \beta_0*intercept + \beta_1*gain + \beta_2*loss\] where the intercept can be interpreted as the overall tendency towards gambling.

The goal of these simulations is to determine how to arrange the gambles and the model such that a subject with equal utility and risk preference would exhibit equal model parameters across the two gamble distributions.

First we set up the gamble structures. For this example we assume expected value computed via prospect theory with a loss aversion lambda of 2 and linear value functions for gains and losses. We start using identical gain distributions across conditions.

msize=16
pt_lambda=2
gains_indiff=seq(from=10,length.out=16,by=2)
loss_indiff=seq(from=5,length.out=16,by=1)

gains_equal=gains_indiff
loss_equal=gains_equal

df=data.frame(gain_indiff=kronecker(array(1,dim=16),gains_indiff),
                  loss_indiff=kronecker(loss_indiff,array(1,dim=16)),
              gain_equal=kronecker(array(1,dim=16),gains_equal),
                  loss_equal=kronecker(loss_equal,array(1,dim=16))
              )

Now let’s simulate the subject’s brain responses to these gambles. For now we just assume that the brain response is a function of prospect theory EV, and we fit a linear regression model to each dataset to assess parameter recovery. First generate a function to do this.

library(QuantPsyc) # for lm.beta
## Loading required package: boot
## Loading required package: MASS
## 
## Attaching package: 'QuantPsyc'
## The following object is masked from 'package:base':
## 
##     norm
domodel <- function(df,noise_sd=1,pt_lambda=2,stdize=TRUE) {
  df$ev_indiff=df$gain_indiff - df$loss_indiff*pt_lambda
  df$ev_equal=df$gain_equal - df$loss_equal*pt_lambda
  df$resp_indiff=df$ev_indiff+rnorm(dim(df)[1])*noise_sd
  df$resp_equal=df$ev_equal+rnorm(dim(df)[1])*noise_sd
  lm.indiff=lm(resp_indiff~gain_indiff + loss_indiff, data=df)
  lm.equal=lm(resp_equal~gain_equal + loss_equal, data=df)
  la_equal=(-1*lm.equal$coefficients[3])/lm.equal$coefficients[2]
  la_indiff=(-1*lm.indiff$coefficients[3])/lm.indiff$coefficients[2]
  if (stdize){
     return(c(lm.beta(lm.equal),
           lm.beta(lm.indiff),
           la_equal,la_indiff))
  } else {
      return(c(lm.equal$coefficients[2],lm.equal$coefficients[3],
           lm.indiff$coefficients[2],lm.indiff$coefficients[3],
           la_equal,la_indiff))

  }
}

Now let’s loop through a range of plausible lambda values and see how well each model can recover the parameters. First do it using unstandardized regression parameters.

lamvals=seq(1,5,by=0.25)
nruns=100
temp=1.5
simresults=c()
for (lam in lamvals){
  for (i in 1:nruns){
    tmp=domodel(df,temp,lam,stdize=FALSE)
    simresults=cbind(simresults,c(lam,tmp))
  }
}

sim.df=data.frame(t(simresults))
names(sim.df)=c('lam',"gain_equal","loss_equal","gain_indiff","loss_indiff",
  "la_equal","la_indiff")

lam_ests=aggregate(la_equal~lam, data=sim.df, FUN=function(x) c(mean=mean(x)))
lam_ests$la_indiff=aggregate(la_indiff~lam, data=sim.df, FUN=function(x) c(mean=mean(x)))$la_indiff
par(mfrow=c(2,2))
plot(sim.df$lam,sim.df$la_indiff,pch=20,main='lambda: indifference',xlab='true lambda',ylab='estimated lambda')
lines(c(1,5),c(1,5))
plot(sim.df$lam,sim.df$la_equal,pch=20,main='lambda: equal',xlab='true lambda',ylab='estimated lambda')
lines(c(1,5),c(1,5))

plot(sim.df$gain_equal,sim.df$gain_indiff,pch=20,main='gain responses',xlab='equal',ylab='indiff',xlim=c(0.9,1.1),ylim=c(0.9,1.1))
#lines(c(0,1),c(0,1))

plot(sim.df$loss_equal,sim.df$loss_indiff,pch=20,main='loss responses',xlab='equal',ylab='indiff',xlim=c(-5,-1),ylim=c(-5,-1))
lines(c(-5,-1),c(-5,-1))

Compare slopes of loss responses between models

loss_lm=lm(loss_equal~loss_indiff,data=sim.df)
print(loss_lm)
## 
## Call:
## lm(formula = loss_equal ~ loss_indiff, data = sim.df)
## 
## Coefficients:
## (Intercept)  loss_indiff  
##  -7.024e-06    9.996e-01

We see that the parameters match well when using unstandardized regression coefficients:

Now do the same using standardized regression parameters.

simresults=c()
for (lam in lamvals){
  for (i in 1:nruns){
    tmp=domodel(df,temp,lam,stdize=TRUE)
    simresults=cbind(simresults,c(lam,tmp))
  }
}

sim.df.std=data.frame(t(simresults))
names(sim.df.std)=c('lam',"gain_equal","loss_equal","gain_indiff","loss_indiff",
  "la_equal","la_indiff")

lam_ests=aggregate(la_equal~lam, data=sim.df.std, FUN=function(x) c(mean=mean(x)))
lam_ests$la_indiff=aggregate(la_indiff~lam, data=sim.df.std, FUN=function(x) c(mean=mean(x)))$la_indiff
par(mfrow=c(2,2))
plot(sim.df.std$lam,sim.df.std$la_indiff,pch=20,main='lambda: indifference',xlab='true lambda',ylab='estimated lambda')
lines(c(1,5),c(1,5))
plot(sim.df.std$lam,sim.df.std$la_equal,pch=20,main='lambda: equal',xlab='true lambda',ylab='estimated lambda')
lines(c(1,5),c(1,5))

plot(sim.df.std$gain_equal,sim.df.std$gain_indiff,pch=20,main='gain responses',xlab='equal',ylab='indiff',xlim=c(0,1),ylim=c(0,1))
lines(c(0,1),c(0,1))

plot(sim.df.std$loss_equal,sim.df.std$loss_indiff,pch=20,main='loss responses',xlab='equal',ylab='indiff',xlim=c(-1,0),ylim=c(-1,0))
lines(c(-1,0),c(-1,0))

loss_lm=lm(loss_equal~loss_indiff,data=sim.df.std)
print(loss_lm)
## 
## Call:
## lm(formula = loss_equal ~ loss_indiff, data = sim.df.std)
## 
## Coefficients:
## (Intercept)  loss_indiff  
##     -0.5010       0.5317

This shows us that the analysis performs incorrectly when using standardized regression parameters:

However, loss aversion coefficients are still properly estimated.