# functions to generate the data
ilogit <- function(x) {
return(1 / (1 + exp(-x)))
}
fun_prob55 <- function(x) {
P = rep(0, length(x))
P[x >= 0.] = ilogit((8.5 * x[x >= 0.] - 1.5)) / 10.5 + 0.65 - 0.0007072
P[x < 0.] = (x[x < 0.] + 1)^4 / 15 + 0.05
return(P)
}
sample_prob55 <- function(x) {
P = rep(0, length(x))
P[x >= 0.] = ilogit((8.5 * x[x >= 0.] - 1.5)) / 10.5 + 0.65 - 0.0007072
P[x < 0.] = (x[x < 0.] + 1)^4 / 15 + 0.05
t = rep(0, length(x))
for (j in 1:length(x)) {
t[j] = sample(c(1, 0), 1, prob = c(P[j], 1 - P[j]))
}
return(t)
}
## Toy example - for the function check only! ##
# data generation
N=100
x = sort(runif(N, -1, 1))
y = sample_prob55(x)
c = 0
# running LoTTA model on a sharp RDD with a binary outcome
out = LoTTA_sharp_BIN(x, y, c, burnin = 50, sample = 50, adapt = 10,n.chains=1
,method = 'simple',inits = NA)
## Use case example ##
# \donttest{
# data generation
N=1000 # try different dataset size
x = sort(runif(N, -1, 1))
y = sample_prob55(x)
c = 0
# running LoTTA function on sharp RDD with binary outcomes;
# cutoff = 0, treatment effect = 0.55
# remember to check convergence and adjust burnin, sample and adapt if needed
out = LoTTA_sharp_BIN(x, y, c, burnin = 10000, sample = 5000, adapt = 1000,n.chains=4)
# print effect estimate:
out$Effect_estimate
# print JAGS output to asses convergence (the output is for normalized data)
out$JAGS_output
# plot posterior fit
LoTTA_plot_outcome(out,nbins = 60)
# }
Run the code above in your browser using DataLab