DAG reconstruction

pacman::p_load(tidyverse, brms)

seed <- 1981 # Defining a seed so the results are always the same
samplesize <- 300 # Defining the amount of datapoints
n <- 1000 # Defining number of simulations
MaxScore <- 8 # Defining max rating
MinScore <- 1 # Defining min rating

## Regression to the mean?
Regression2Mean <- 0.7 # 1st to 2nd in the empirical data : 0.668444
## Defining the correlation coefficients
Conformity <- 0.5 # Defining the true correlation of feedback and change

SimulateData <- function(samplesize,
                         Regression2Mean,
                         Conformity,
                         MinScore = 1,
                         MaxScore = 8){
  
  FirstRating <- round(runif(samplesize, MinScore, MaxScore), 0)
  Feedback <- round(runif(samplesize, -3, 3), 0)
  SecondRating <- round(Regression2Mean * FirstRating + Conformity * Feedback)
  SecondRating <- ifelse(SecondRating > MaxScore, MaxScore,ifelse(SecondRating < MinScore, MinScore, SecondRating))
  Change <- SecondRating - FirstRating
    
  d1 <- data.frame(FirstRating, Feedback, SecondRating, Change) %>% 
    subset(FirstRating + Feedback < MaxScore & FirstRating + Feedback > MinScore) %>%
    mutate(
      FirstRatingC <- FirstRating - 4.5,
      SecondRatingC <- SecondRating - 4.5,
    )
  return(d1)
}

## Simulating data
d0 <- SimulateData(samplesize, Regression2Mean = .7, Conformity = 0, MinScore = 1, MaxScore = 7)
d1 <- SimulateData(samplesize, Regression2Mean = .7, Conformity = .1, MinScore = 1, MaxScore = 7)
d2 <- SimulateData(samplesize, Regression2Mean = .7, Conformity = .2, MinScore = 1, MaxScore = 7)
d3 <- SimulateData(samplesize, Regression2Mean = .7, Conformity = .3, MinScore = 1, MaxScore = 7)
d4 <- SimulateData(samplesize, Regression2Mean = .7, Conformity = .4, MinScore = 1, MaxScore = 7)
d5 <- SimulateData(samplesize, Regression2Mean = .7, Conformity = .5, MinScore = 1, MaxScore = 7)

# Creating possible models

Conformity_f1 <- bf(Change ~ 1 + Feedback)
Conformity_f2 <- bf(Change ~ 1 + Feedback + FirstRating)
Conformity_f3 <- bf(Feedback ~ 1 + FirstRating) + bf(Change ~ 1 + Feedback + FirstRating) + set_rescor(F)

prior <- c(
  prior(normal(0,1), class = Intercept),
  prior(normal(0,1), class = b),
  prior(normal(0,1), class = sigma)
)

Conformity_m1_prior <- brm(
  Conformity_f1,
  d1,
  family=gaussian(),
  prior=prior,
  sample_prior="only",
  cores=1,
  chains=2,
  control = list(
    adapt_delta=0.99,
    max_treedepth = 20
  )
)

pp_check(Conformity_m1_prior, nsamples=100)

Conformity_m1_0 <- brm(
  Conformity_f1,
  d0,
  family=gaussian(),
  prior=prior,
  sample_prior=T,
  cores=1,
  chains=2,
  control = list(
    adapt_delta=0.99,
    max_treedepth = 20
  )
)

Conformity_m1_1 <- update(Conformity_m1_0, newdata = d1)
Conformity_m1_2 <- update(Conformity_m1_0, newdata = d2)
Conformity_m1_3 <- update(Conformity_m1_0, newdata = d3)
Conformity_m1_4 <- update(Conformity_m1_0, newdata = d4)
Conformity_m1_5 <- update(Conformity_m1_0, newdata = d5)

Conformity_m2_0 <- update(Conformity_m1_0, formula = Conformity_f2, newdata = d0)
Conformity_m2_1 <- update(Conformity_m2_0, newdata = d1)
Conformity_m2_2 <- update(Conformity_m2_0, newdata = d2)
Conformity_m2_3 <- update(Conformity_m2_0, newdata = d3)
Conformity_m2_4 <- update(Conformity_m2_0, newdata = d4)
Conformity_m2_5 <- update(Conformity_m2_0, newdata = d5)

prior_SEM <- c(
  prior(normal(0,1), class = Intercept, resp=Change),
  prior(normal(0,1), class = Intercept, resp=Feedback),
  prior(normal(0,1), class = b, resp=Change),
  prior(normal(0,1), class = b, resp=Feedback),
  prior(normal(0,1), class = sigma, resp=Change),
  prior(normal(0,1), class = sigma, resp=Feedback)
)

Conformity_m3_0 <- brm(
  Conformity_f3,
  d0,
  family=gaussian(),
  prior=prior_SEM,
  sample_prior=T,
  cores=1,
  chains=2,
  control = list(
    adapt_delta=0.99,
    max_treedepth = 20
  )
)

Conformity_m3_1 <- update(Conformity_m3_0, newdata = d1)
Conformity_m3_2 <- update(Conformity_m3_0, newdata = d2)
Conformity_m3_3 <- update(Conformity_m3_0, newdata = d3)
Conformity_m3_4 <- update(Conformity_m3_0, newdata = d4)
Conformity_m3_5 <- update(Conformity_m3_0, newdata = d5)



paste0(
  round(fixef(Conformity_m1_0)[,"Estimate"][[2]],2), 
  " SE: ", 
  round(fixef(Conformity_m1_0)[,"Est.Error"][[2]],2))
paste0(
  round(fixef(Conformity_m2_0)[,"Estimate"][[2]],2), 
  " SE: ", 
  round(fixef(Conformity_m2_0)[,"Est.Error"][[2]],2))
paste0(
  round(fixef(Conformity_m3_0)[,"Estimate"][[4]],2), 
  " SE: ", 
  round(fixef(Conformity_m3_0)[,"Est.Error"][[4]],2))

paste0(
  round(fixef(Conformity_m1_1)[,"Estimate"][[2]],2), 
  " SE: ", 
  round(fixef(Conformity_m1_1)[,"Est.Error"][[2]],2))
paste0(
  round(fixef(Conformity_m2_1)[,"Estimate"][[2]],2), 
  " SE: ", 
  round(fixef(Conformity_m2_1)[,"Est.Error"][[2]],2))
paste0(
  round(fixef(Conformity_m3_1)[,"Estimate"][[4]],2), 
  " SE: ", 
  round(fixef(Conformity_m3_1)[,"Est.Error"][[4]],2))

paste0(
  round(fixef(Conformity_m1_2)[,"Estimate"][[2]],2), 
  " SE: ", 
  round(fixef(Conformity_m1_2)[,"Est.Error"][[2]],2))
paste0(
  round(fixef(Conformity_m2_2)[,"Estimate"][[2]],2), 
  " SE: ", 
  round(fixef(Conformity_m2_2)[,"Est.Error"][[2]],2))
paste0(
  round(fixef(Conformity_m3_2)[,"Estimate"][[4]],2), 
  " SE: ", 
  round(fixef(Conformity_m3_2)[,"Est.Error"][[4]],2))