Systematic Error

Load file

load(file = "C:/Users/DELL/Documents/integration-simulations/index.RData")
library(pracma)
library(ggplot2)

Estimate Parameters

plottingInfoAll = data.frame(Outcome = rep(c("Incorrect & N.S.", 
                                             "Correct & N.S.", 
                                             "Incorrect & p < 0.05",
                                             "Correct & p < 0.05"), times = 5),
                             Count = 0,
                             factor = rep(c(0, 0.025, 0.05, 0.1, 0.15), each = 4))

sampleSize = 60
Factors = c(0, 0.025, 0.05, 0.1, 0.15)
nTrials = 60
for (factor in Factors){
  df = n1000[sample(200), ]
  df = as.data.frame(df)
  df$NLL = 0
  df$NLL_alt = 0
  for (i in 1:nrow(df)){
    row = which(round(freeParameters$theta, 2) == round(df$theta[i], 2) & 
                  round(freeParameters$phi, 2) == round(df$phi[i], 2))
    preds = as.numeric(predictions[row, ])
    preds = rep(preds, times = nTrials/length(preds))
    data = rbind(trialSet, trialSet)
    Xmax = max(data$Endowment/2)
    data$Transferred = noise(preds, data$Endowment)
    data$Transferred = round(data$Transferred + 
      (((Xmax - (data$Endowment/2))/Xmax) * factor))
    
    result = tryCatch({fmincon(obj_function,x0 = initial_params,
                       lb = lower_bounds, ub = upper_bounds,
                       df = data)}, error = function(e){data.frame(par = NA)})
    resultAlt = tryCatch({optimizePerCondition(df = data)}, 
                         error = function(e){data.frame(par = NA)})
    if (is.na(resultAlt$par[1]) | is.na(result$par[1])){
      df$NLL[i] = NA; df$NLL_alt[i] = NA; next
    }
    data$Preds = generatePredictions(result$par, data)
    data$PredsAlt = generateAltPerCondition(resultAlt$par, data)
    
    if (sum(data$Preds == data$Transferred)==nrow(data)){
      data$Preds[sample(nrow(data), 1)] =+ 1
    }
    if (sum(data$PredsAlt == data$Transferred)==nrow(data)){
      data$PredsAlt[sample(nrow(data), 1)] =+ 1
    }
    
    df$NLL[i] = -sum(dnorm(data$Transferred/data$Endowment, 
                           mean = data$Preds/data$Endowment, 
                           sd = sd((data$Transferred - data$Preds)/data$Endowment),
                           log = T))
    df$NLL_alt[i] = -sum(dnorm(data$Transferred/data$Endowment, 
                               mean = data$PredsAlt/data$Endowment, 
                           sd = sd((data$Transferred - data$PredsAlt)/data$Endowment), log = T))
  }
  df$BIC = 2 * df$NLL + 2 * log(nrow(data))
  df$BIC_alt = 2 * df$NLL_alt + 3 * log(nrow(data))
  if ((sum(is.na(df$NLL)) + sum(is.na(df$NLL_alt))) > 0){
    df = df[-which(is.na(df$NLL) | is.na(df$NLL_alt)), ]
  }
  test = vector('logical', 100000)
  better = vector('logical', 100000)
  for (i in 1:length(test)){
    subsample = df[sample(1:nrow(df), sampleSize, replace = T), ]
    better[i] = sum(subsample$BIC) < sum(subsample$BIC_alt)
    test[i] = t.test(subsample$BIC, subsample$BIC_alt, paired = T)$p.value < 0.05
  }
  plottingInfoAll$Count[which(plottingInfoAll$factor == factor & 
                              plottingInfoAll$Outcome == "Incorrect & N.S.")] = 
    sum(better==F & test==F)
  plottingInfoAll$Count[which(plottingInfoAll$factor == factor & 
                              plottingInfoAll$Outcome == "Correct & N.S.")] = 
    sum(better==T & test==F)
  plottingInfoAll$Count[which(plottingInfoAll$factor == factor & 
                              plottingInfoAll$Outcome == "Incorrect & p < 0.05")] = 
    sum(better==F & test==T)
  plottingInfoAll$Count[which(plottingInfoAll$factor == factor & 
                              plottingInfoAll$Outcome == "Correct & p < 0.05")] = 
    sum(better==T & test==T)
}

And plot

plottingInfoAll$factor = factor(plottingInfoAll$factor, 
                                levels = c(0.15, 0.1, 0.05, 0.025, 0))
plottingInfoAll$Outcome = factor(plottingInfoAll$Outcome,
                                 levels = c("Correct & p < 0.05",
                                            "Correct & N.S.", 
                                            "Incorrect & N.S.",
                                            "Incorrect & p < 0.05"))
plottingInfoAll$Proportion = plottingInfoAll$Count/length(test)
plottingInfoAll$Position = plottingInfoAll$Proportion/2

these = which(plottingInfoAll$Outcome == "Correct & p < 0.05")
plottingInfoAll$Position[these] = 1 - plottingInfoAll$Proportion[these]/2

those = which(plottingInfoAll$Outcome == "Correct & N.S.")
plottingInfoAll$Position[those] = 1 - plottingInfoAll$Proportion[these] -
  plottingInfoAll$Proportion[those]/2

these = which(plottingInfoAll$Outcome == "Incorrect & p < 0.05")
those = which(plottingInfoAll$Outcome == "Incorrect & N.S.")
plottingInfoAll$Position[those] = max(plottingInfoAll$Proportion[these] +
  plottingInfoAll$Proportion[those]/2, 0.2)

final = plottingInfoAll[-which(plottingInfoAll$Proportion == 0),]

ggplot(data = final, aes(x = Proportion, y = factor, group = Outcome,
                                   fill = Outcome)) + 
  geom_col() + theme_minimal() + 
  scale_fill_manual(values = c("green", "forestgreen", "firebrick4", "red")) +
  labs(x = "Proportion", y = "Bias Factor", 
       fill = "Result",
       label = NULL) + 
  geom_label(aes(label = Proportion, x = Position)) + lims(x = c(-0.025, 1.25))

Now save workspace for later reference

save.image(file = "C:/Users/DELL/Documents/integration-simulations/bias.RData")