I am new to R and I am having some struggle trying to manipulate a model published in a scientific article to see what predictions it makes for certain situations.
This is a model of 'causal judgment' published in a psych paper. It tries to predict in what situations we are inclined to say that 'X caused Y' depending on the priors of relevant variables and the causal function linking them to the outcome.
In a nutshell, the model says that when producing a causal judgments regarding a real situation: 1- Sample a series of counterfactual worlds, where every variable connected to the outcome of interest is given a value that depends on the priors for that variable, and the value that it takes in the real world (with a certain probability s, the variable is given the value that it takes in the real world, and with probability 1-s, its value is sampled from the prior). 2- Compute a measure of correlation between each causal variable and the outcome across counterfactual worlds. This is going to give you a measure of the 'causal strength' that links the cause and the outcome, which predicts the extent to which you are inclined to say that 'X caused Y' in a given situation.
The model is used to predict the results of experiments involving urns from which a player draws colored/uncolored balls (according to certain priors), where the player wins or lose depending on the rule of that game (the outcome).
Below is the code that implements the simulation model:
draw <- function(n, prob_x){
x_prior <- rbinom(n, 1, prob_x)
return(x_prior)
}
draw(3, c(.1,.5,.5) )
#the function to run a simulation
simulation_altS <- function(s, actual_world, threshold){
#draw a ball from each urn and record what you have
urn_draws_prior <- draw(length(probs), probs)
where_to_draw <- rbinom(1,1,s)
if(where_to_draw){urn_draws <- actual_world}
else{urn_draws <- urn_draws_prior}
#compute how many points the player gets
contributions <- urn_draws
#if points are above a threshold, the player wins
#outcome <- sum(contributions) >= threshold
outcome <- (a&!b)|(!a&b)
#return the list of all balls that were drawn, and
#whether the player won,
return(list("draws" = urn_draws, "outcome" = outcome))
}
simulation_altS(.5, actual_world, threshold)
###a function to run simulations and compute causal scores:
run_simulations_altS <- function(actual_world, probs, threshold, s, m=4){
##lists where we will store the simulation data
#this list contains the relevant data for the CESM
simulations <- vector("list", m)
##run the simulations
for (i in 1:n){
sim_global <- simulation_altS(s, actual_world, threshold)
#cycle through all urns
for (urn in 1:(m-1)){
#record whether a colored ball was drawn from the urn
simulations[[urn]][i] <- sim_global$draws[urn]
}
#record whether the player won
simulations[[m]][i] <- sim_global$outcome
}
##compute the causal score for each urn, CESM
cor_list <- integer(length(actual_world))
#cycle through all urns
for (a in 1:(m-1)){
#if a colored ball was drawn from the urn, compute a causal score
#otherwise return NA
cor_list[a] <- ifelse(actual_world[a] == 1,
cor(simulations[[a]], simulations[[m]]), NA)
}
#store the list of causal scores in different vectors
mean_outcome = mean(simulations[[m]])
#gather our causal scores in a dataframe
df <- data.frame(1:(m-1), cor_list)
colnames(df) <- c("urn", "cesm")
#return the dataframe
return(list("df"=df,
"mean_outcomes"=mean_outcome))
}
And below is an example of a study where certain prior probabilities and real-world values are being fed to the model:
draw <- function(n, prob_x){
x_prior <- rbinom(n, 1, prob_x)
return(x_prior)
}
draw(3, c(.1,.5,.5) )
#the function to run a simulation
simulation_altS <- function(s, actual_world, threshold){
#draw a ball from each urn and record what you have
urn_draws_prior <- draw(length(probs), probs)
where_to_draw <- rbinom(1,1,s)
if(where_to_draw){urn_draws <- actual_world}
else{urn_draws <- urn_draws_prior}
#compute how many points the player gets
contributions <- urn_draws
#if points are above a threshold, the player wins
#outcome <- sum(contributions) >= threshold
outcome <- (a&!b)|(!a&b)
#return the list of all balls that were drawn, and
#whether the player won,
return(list("draws" = urn_draws, "outcome" = outcome))
}
simulation_altS(.5, actual_world, threshold)
###a function to run simulations and compute causal scores:
run_simulations_altS <- function(actual_world, probs, threshold, s, m=4){
##lists where we will store the simulation data
#this list contains the relevant data for the CESM
simulations <- vector("list", m)
##run the simulations
for (i in 1:n){
sim_global <- simulation_altS(s, actual_world, threshold)
#cycle through all urns
for (urn in 1:(m-1)){
#record whether a colored ball was drawn from the urn
simulations[[urn]][i] <- sim_global$draws[urn]
}
#record whether the player won
simulations[[m]][i] <- sim_global$outcome
}
##compute the causal score for each urn, CESM
cor_list <- integer(length(actual_world))
#cycle through all urns
for (a in 1:(m-1)){
#if a colored ball was drawn from the urn, compute a causal score
#otherwise return NA
cor_list[a] <- ifelse(actual_world[a] == 1,
cor(simulations[[a]], simulations[[m]]), NA)
}
#store the list of causal scores in different vectors
mean_outcome = mean(simulations[[m]])
#gather our causal scores in a dataframe
df <- data.frame(1:(m-1), cor_list)
colnames(df) <- c("urn", "cesm")
#return the dataframe
return(list("df"=df,
"mean_outcomes"=mean_outcome))
} # study 2a
probs <- c(.05,.5,.95)
actual_world <- c(1, 1, 1)
threshold <- 2
Here the study would be one involving three urns, one with very low probability (0.05), one intermediate (.5), and one with very high probability (0.95) of drawing a colored ball from that urn, where the rule of the game is that you need to draw at least 2 balls (threshold <- 2) to win, and where in the real world situation presented to the experimental subject, the players happens to have drawn a colored ball from each of the three urns ( actual_world <- c(1, 1, 1) ).
In general, the model makes very accurate predictions wrt. subjects behavior for monotonic causal functions. For example it predicts that subject will tend to think of low-probability events as more causal when the causal function linking some variables to an outcome is conjunctive (you need A and B to get E), whereas they would favor high-probability events in disjunctive settings.
Now what I want to test what are the predictions of this model for non-monotonic causal functions, such as for example a game where a player needs to draw 'exactly 1' colored ball to win the game (meaning they lose if they draw, 0, 2 or 3 balls instead).
My problem is that because all outcome functions defined so far are monotonic, the relationship between draws and outcomes in the model is defined by a 'threshold' parameter, which basically counts the number of colored balls drawn from urns in a given experimental situation, and makes the draw a winning one if that number is above a certain threshold:
(snippet from code above:)
In general, the model makes very nice predictions of subjects behavior #compute how many points the player gets
contributions <- urn_draws
#if points are above a threshold, the player wins
outcome <- sum(contributions) >= threshold
#return the list of all balls that were drawn, and
#whether the player won,
return(list("draws" = urn_draws, "outcome" = outcome))
}
simulation_altS(.5, actual_world, threshold)
I would like to change that to make the outcome function a non monotonic one, like the 'exactly one' function for three urns.
One way I imagine this could be done is by giving a label (a, b, c) to the urns in the experiment, and then writing the outcome as logical formula, like
outcome <- (a&!b&!c)|(!a&b&!c)|(!a&!b&c)
in order to indirectly implement the 'exactly 1' function, by this would require making substantial changes upstream of the code, which is now tailored to count the number of draws and define the outcome as a function of that.
So I am at loss to find a solution to fix this without messing up the entire function, all the more so that I have very little experiment with R, so I wouldn't want to inadvertently twist the code in a way that would make whatever results I get non-representative of the model.
If you ever need more detailed info about the code or are interested in the article, the preprint is available here: https://psyarxiv.com/ts76y/ together with supplemental materials (which include all of the code used for their model).
Any help would be much appreciated, thanks a lot in advance !!
CK