### example from the last lecture using rethinking package
library(devtools) # install.packages("devtools")
devtools::install_github("rmcelreath/rethinking") # https://github.com/rmcelreath/rethinking
library(rethinking)

# some Stan settings
options(mc.cores = parallel::detectCores())
rstan_options(auto_write = TRUE)
Sys.setenv(LOCAL_CPPFLAGS = '-march=native')
# 69/81 reported white in CT
# 54/79 reported white in RT
data_cheating <- list(
  whites = c(69,54),
  tosses = c(81,79),
  condition = c(1,2)
)

fit_cheating <- map( 
  alist(
    whites ~ dbinom( tosses , p ),
    p <- cond[condition],
    cond[condition] ~ dunif( 0 , 1 )
  ), 
  data=data_cheating
)
precis(fit_cheating, depth = 2)

# extract posterior distribution
posterior_cheating <- extract.samples(fit_cheating)
str(posterior_cheating)

# plot the results
par(mfrow = c(1,2))
hist(posterior_cheating$cond[,1], xlim = c(0,1), main = "probability of white in CT")
hist(posterior_cheating$cond[,2], xlim = c(0,1), main = "probability of white in RT")
dev.off()

# with stan
data_cheating2 <- list(
  whites = as.integer(c(69,54)),
  tosses = as.integer(c(81,79)),
  condition = as.integer(c(1,2))
)

fit_cheating2 <- map2stan( 
  alist(
    whites ~ dbinom( tosses , p ),
    p <- cond[condition],
    cond[condition] ~ dunif( 0, 1 )
  ), 
  data=data_cheating2,
  start = list(cond = c(.5,.5)),
  iter = 2000, warmup = 1000, chains = 3
)
plotchains(fit_cheating2)
precis(fit_cheating2, depth = 2)
stancode(fit_cheating2)

precis(fit_cheating, depth = 2)

# extract posterior distribution
posterior_cheating2 <- extract.samples(fit_cheating2)
str(posterior_cheating2)

# plot the results
par(mfrow = c(1,2))
hist(posterior_cheating2$cond[,1], xlim = c(0,1), main = "probability of white in CT")
hist(posterior_cheating2$cond[,2], xlim = c(0,1), main = "probability of white in RT")
dev.off()

dif_cond <- posterior_cheating2$cond[,1] - posterior_cheating2$cond[,2]
hist(dif_cond)

HPDI(dif_cond, prob = .95)
quantile(dif_cond, c(.025,.5,.975))
  

### another example
data_overconfidence <- read.table("Atir Rosenzweig Dunning 2015 Study 1b.csv", sep = ",", header = T)
colnames(data_overconfidence)[1] <- "id"
head(data_overconfidence)

data_over <- list(
  "overclaiming"        = data_overconfidence$overclaiming_proportion,
  "perceived_knowledge" = data_overconfidence$self_perceived_knowledge,
  "FINRA"               = data_overconfidence$FINRA_score,
  "condition"           = data_overconfidence$order_of_tasks
)

data_over$overclaiming <- as.vector(scale(data_over$overclaiming))
data_over$perceived_knowledge <- as.vector(scale(data_over$perceived_knowledge))
data_over$FINRA <- as.vector(scale(data_over$FINRA))
data_over$condition <- data_over$condition - 1

summary(lm(overclaiming~perceived_knowledge+FINRA+condition, data = data_over))
summary(lm(overclaiming~perceived_knowledge, data = data_over))

fit_overclaiming <- map2stan( 
  alist(
    overclaiming ~ dnorm( mu , sigma ),
    mu <- a + b_perceived_knowledge * perceived_knowledge + b_FINRA * FINRA + b_condition * condition,
    
    a ~ dnorm(0, 1),
    b_perceived_knowledge ~ dnorm(0, 1),
    b_FINRA  ~ dnorm(0, 1),
    b_condition  ~ dnorm(0, 1),
    
    sigma ~ dnorm( 0, 1 )
  ), 
  data=data_over,
  iter = 2000, warmup = 1000, chains = 3
)

plotchains(fit_overclaiming)
precis(fit_overclaiming, prob = .95)
stancode(fit_overclaiming)


### Z-curve - p-values and etc...
library("nleqslv")
# computing z-scores acording to power
solve_Z <- function(x,p){
  y = numeric(1)
  y = 1 - pnorm(qnorm(0.975), x) + pnorm(qnorm(0.025), x) - p
  y
}

pwr = .8
nleqslv(.05,solve_Z, p = pwr, control = list(stepmax = .1))


mu = 2.81
pnorm(1.96, mu, 1, lower.tail = F) + pnorm(-1.96, mu, 1)


zval <- rnorm(10000,2.81)
zval <- abs(zval)
hist(zval)

p.val <- pnorm(zval, lower.tail = F)*2
hist(p.val)

mean(p.val < .05)

z_obs <- zval[zval > 1.96]
hist(z_obs)

z_data <- list(
  "z" = z_obs,
  "N" = length(z_obs)
)

fit_zdist <- stan("z_dist.stan", data = z_data, chains = 3, warmup = 1000, iter = 2000)
traceplot(fit_zdist)
summary(fit_zdist)$summary
round(summary(fit_zdist)$summary,3)

# adding zero effects
z_null <- rnorm(1000)
z_null <- abs(z_null)
hist(z_null)

z1 <- rnorm(1000,2.81)
z1 <- abs(z1)
hist(z1)

zval <- c(z_null,z1)
hist(zval)

z_obs <- zval[zval > 1.96]
hist(z_obs)

z_data <- list(
  "z" = z_obs,
  "N" = length(z_obs)
)
fit_zdist2 <- stan("z_dist2.stan", data = z_data, chains = 1, warmup = 1000, iter = 2000)
traceplot(fit_zdist2)
summary(fit_zdist2)$summary
round(summary(fit_zdist2)$summary,3)