Commit 990f82b2 authored by Simon Ciranka's avatar Simon Ciranka

New Report

parent 47284d71
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
......@@ -123,7 +123,7 @@ for(i in 1:length(TidyMarble$rt)){
if(as.numeric(as.character(TidyMarble$age[i]))<18){
TidyMarble$Agegroup[i]=1#are they adolescents?
if(as.numeric(as.character(TidyMarble$age[i]))<12){
if(as.numeric(as.character(TidyMarble$age[i]))<=12){
TidyMarble$Agegroup[i]=0#are they kids?
}#enddoubleif.
}else TidyMarble$Agegroup[i]=2#ore adults?
......@@ -287,7 +287,7 @@ for(i in 1:length(TidyMarble$rt)){
if(as.numeric(as.character(TidyMarble$age[i]))<18){
TidyMarble$Agegroup[i]=1#are they adolescents?
if(as.numeric(as.character(TidyMarble$age[i]))<12){
if(as.numeric(as.character(TidyMarble$age[i]))<=12){
TidyMarble$Agegroup[i]=0#are they kids?
}#enddoubleif.
}else TidyMarble$Agegroup[i]=2#ore adults?
......
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
#here you need to specify how many instances you want to run simulations on.
#currently it si set up as such that for any "real" collected data from our pilot there will be Simulations.
#this is not to fit a model yet. this makes extensive Simulations. I fit the model back at another place.
HowManyGroups=3
HowManyModels=1
#j=2;
#i=1;
#k=5;
for i in `seq 1 $HowManyGroups`;
do
for j in `seq 1 $HowManyModels`;
do
echo "$i $Subjects";
#here i specify the job for the cluster.
#for input_file in INPUT/* ; do
#echo "#PBS -m n" > job.pbs
echo "#PBS -N BayesianLearning$i$j" > job.pbs
echo "#PBS -o /logs/" >> job.pbs
echo "#PBS -l mem=15gb" >> job.pbs
echo "#PBS -j oe" >> job.pbs
echo "#PBS -l walltime=100:00:0" >> job.pbs
echo "#PBS -d ." >> job.pbs
echo "#PBS -l nodes=1:ppn=4" >>job.pbs
echo "module load R/3.6.0; Rscript stanScriptMarbles.R $i $j" >> job.pbs
qsub job.pbs
rm -f job.pbs
done
done
# this script subsets the data and fits different learning models on the PBS cluster of the MPI
library(tidyverse)
library(rstan)
bashInput <- commandArgs(trailingOnly = TRUE)
#bashInput=c(1,1)
####### Subsetting and Model Determination Block
LearningModelPath<-"../B_Model_Code/LearningModels/ocu_hier_RiskUnc.stan"
# "")
Agegroups=c("0","1","2")
MarbleData<-read.csv("../TidyMarbleNew.csv")
hist(MarbleData$subject,breaks=length(unique(MarbleData$subject))*4)
MarbleData$Agegroup=as.factor(MarbleData$Agegroup)
MarbleData$DFE1DFD0=as.factor(MarbleData$DFE1DFD0)
MarbleData$Social1Ind0=as.factor(MarbleData$Social1Ind0)
MarbleData$PercentBlueEstimate=round(as.numeric(MarbleData$PercentBlueEstimate))
MarbleData$PercentBlueShownRel=round(as.numeric(MarbleData$PercentBlueShownRel),2)
MarbleData$PercentBlueShownRel=as.factor(as.character(MarbleData$PercentBlueShownRel))
MarbleData$HowSure=as.numeric(MarbleData$HowSure)
# Set Priors and get everything ready for stan
MarbleData<-MarbleData%>%filter(Agegroup==Agegroups[as.numeric(bashInput[1])])#Bin_Update_Data<-Bin_Update_Data[Bin_Update_Data$typeRA=="1",]# keep only The Risk trails.
Bin_Update_Data<-MarbleData%>%arrange(Agegroup)# i order it first so i can make sure that
Agegroups<-unique(Bin_Update_Data$Age.bins)# for indexing my agegroups.
#make it fit with my stan file.
Bin_Update_Data<-Bin_Update_Data%>%dplyr::mutate(
OtherChoseRisk = case_when(
OtherChoseRisk=="NULL" ~ 2,# i dont need this but i restricted the numbers in stan between 1 and 3
OtherChoseRisk=="1" ~ 3,# risky choices are coded as 3 in my stan code
OtherChoseRisk=="0" ~ 1,# safe choices are coded as 1 in my stan code
TRUE~0 # keep the rest.
)
)#end PeerChoice.
#subset it to fit on only one agegroup.
#Bin_Update_Data%>%filter(Agegroup==2)->Bin_Update_Data
#Bin_Update_Data%>%filter(DFE1DFD0==0)->Bin_Update_Data
# now check how many participants we have.
Bin_Update_Data$subject<-as.numeric(Bin_Update_Data$subject)
#change colname of subject into subjID
numSubjs<-length(unique(Bin_Update_Data$subject))#Total Number of Subs
subjList <- unique(Bin_Update_Data$subject)######
Sequence_Length <- unique(Bin_Update_Data[Bin_Update_Data$TotalNShown<99,]$TotalNShown)
Sequence_Length<-9
####### number of trials and see which group he is in for each subject######
Tsubj <- as.vector(rep(0, numSubjs))
for (sIdx in 1:numSubjs) {
curSubj <- subjList[sIdx]
Tsubj[sIdx] <- length(Bin_Update_Data[Bin_Update_Data$subject==curSubj,]$subject) # How many entries per Subject?
#GSubj[sIdx] <- Simulations[[i]]$group
}
maxTrials <- max(Tsubj)
# Information for user continued
cat(" # of (max) trials per subject = ", maxTrials, "\n\n")
# for multiple subjects
condition <- array(0, c(numSubjs, maxTrials))#Other Chose Risk?
p_gamble <- array(0, c(numSubjs, maxTrials))
choice <- array(0, c(numSubjs, maxTrials))
safe_payoff<- array(0, c(numSubjs, maxTrials))
risky_payoff<- array(0, c(numSubjs, maxTrials))
risk1Unc0 <- array(0, c(numSubjs, maxTrials))
#specs for the uncertain trails
Sucess <- array(0, c(numSubjs, maxTrials,Sequence_Length))
Fail <- array(0, c(numSubjs, maxTrials,Sequence_Length))
p_gamble_est<- array(0, c(numSubjs, maxTrials))
Blue_Sum_All<- array(0, c(numSubjs, maxTrials))
Red_Sum_All<- array(0, c(numSubjs, maxTrials))
#Agegroup<- array(0, c(numSubjs, maxTrials))
# generate the data Lists to be passed to stan
# concatenate different groups in the third dimension.
for (i in 1:numSubjs) {
curSubj <- subjList[i]
useTrials <- Tsubj[i]
condition[i, 1:useTrials] <- Bin_Update_Data[Bin_Update_Data$subject==curSubj,]$OtherChoseRisk
p_gamble[i, 1:useTrials] <- Bin_Update_Data[Bin_Update_Data$subject==curSubj,]$probGamb
choice[i, 1:useTrials] <- Bin_Update_Data[Bin_Update_Data$subject==curSubj,]$ChooseRisk
risky_payoff[i, 1:useTrials] <- Bin_Update_Data[Bin_Update_Data$subject==curSubj,]$valueGamble
safe_payoff[i, 1:useTrials] <- 5
risk1Unc0[i, 1:useTrials] <- as.numeric(as.character(Bin_Update_Data[Bin_Update_Data$subject==curSubj,]$DFE1DFD0))
p_gamble_est[i, 1:useTrials]<-as.double(as.character(Bin_Update_Data[Bin_Update_Data$subject==curSubj,]$PercentBlueEstimate))/100
#ambigLevel[i, 1:useTrials] <- Bin_Update_Data[Bin_Update_Data$subject==curSubj,]$Ambiguity
# here i tear up the sequence. its saved as factor so i need to make it into a character string and split it up according to the seperator,
for (t in 1:useTrials){
Blue<-strsplit(toString((Bin_Update_Data[Bin_Update_Data$subject==curSubj,]$blue_marbles[t])),split=",")
BlueSum<-0;
# now i add up the number of sucesses
for (k in 1:length(Blue[[1]])){
BlueSum=BlueSum+as.numeric(Blue[[1]][k])
}
Blue_Sum_All[i, t]<-BlueSum
#here i add up the number of failures
Red<-strsplit(toString((Bin_Update_Data[Bin_Update_Data$subject==curSubj,]$red_marbles[t])),split=",")
RedSum<-0;
for (k in 1:length(Blue[[1]])){
RedSum=RedSum+as.numeric(Red[[1]][k])
}
Red_Sum_All[i, t]<-RedSum
# this isfor sequential updating where binary sucesses and fails are stored in the thrid dimension of the matrix
if(RedSum!=99){
Sucess[i, t,1:Sequence_Length]<-c(array(1, c(BlueSum)), array(0, c(Sequence_Length-BlueSum)) )
Fail[i, t,1:Sequence_Length]<-c(array(1, c(RedSum)), array(0, c(Sequence_Length-RedSum)) )
}else{
# if it was a risk trial, then either put in 99 or make sure that you sample from uniform.
Sucess[i, t,1:Sequence_Length]<-99
Red_Sum_All[i, t]<-1
Fail[i, t,1:Sequence_Length]<-99
Blue_Sum_All[i, t]<-1
}
}
}
#cant be bigger than 1.
p_gamble_est[p_gamble_est>=1]<-0.99
# Specify the number of parameters and parameters of interest
numPars <- 5
POI <- c("mu_rho", "mu_tau","mu_ocu_Uncertainty","mu_ocu_Risk",
"mu_alpha_add","mu_beta_add",
"ocu_Uncertainty","ocu_Risk","alpha_add","beta_add",
"rho","tau") #"log_lik","y_pred")
#this is all i need for stan.
dataList <- list(
N = numSubjs,## number of subjects in each group.
T = maxTrials,
#Seq = Sequence_Length,
Tsubj = Tsubj,
# numPars = numPars,
safe_payoff = safe_payoff,
risky_payoff = risky_payoff,
p_gamble_est= p_gamble_est,
Sucess = Blue_Sum_All,
Fail = Red_Sum_All,
condition = condition, # condition is 0= solo ,1 = safe ,3 = risky
p_gamble = p_gamble,
choice = choice,
risk1Unc0=risk1Unc0
# agegroup=agegroup
)
inits<-"fixed"
# priors
if (inits[1] != "random") {
if (inits[1] == "fixed") {
inits_fixed <- c(0.5, 0.9, 0.0 ,0.5, 0.5)
} else {
if (length(inits) == numPars) {
inits_fixed <- inits
# mu_ocu =rep(inits_fixed[3],2)
} else {
stop("Check your inital values!")
}
}
genInitList <- function() {
list(
#initial values.
mu_rho =rep(inits_fixed[1]),
mu_tau =rep(inits_fixed[2]),
mu_ocu_Risk =c(0.1),
mu_ocu_Uncertainty =c(0.1),
mu_alpha_add=inits_fixed[4],
mu_beta_add=inits_fixed[5],
sigma_rho= rep(c(1.0)),
sigma_tau= rep(c(1.0)),
sigma_ocu_Risk= rep(c(1.0)),
sigma_ocu_Uncertainty= rep(c(1.0)),
sigma_alpha_add=1.0,
sigma_beta_add=1.0,
rho_p = rep(qnorm(inits_fixed[1]/2),numSubjs),
tau_p = rep(log(inits_fixed[2]),numSubjs),
ocu_Risk_p = rep(0.1,numSubjs),
ocu_Uncertainty_p = rep(0.1,numSubjs),
alpha_add_p = rep(qnorm(inits_fixed[4]/2),numSubjs),
beta_add_p = rep(qnorm(inits_fixed[5]/2),numSubjs)
#ocu_p = matrix(inits_fixed[3],numSubjsG, nGroups)
)
}
} else {
genInitList <- "random"
}
ModelFit = stan(LearningModelPath[as.numeric(bashInput[2])],
data = dataList,
pars = POI,
init = genInitList,
iter = 10000,
cores = 4,
chains =4,
control = list(adapt_delta = 0.999)
)
saveRDS(ModelFit,file=paste0("../C_ModelFits/Marbles_Model_",bashInput[2],"Age_",Agegroups[as.numeric(bashInput[1])],".rds"))
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment