Type:	Package
Title:	Bayesian Inference in Regression Discontinuity Designs
Version:	0.1.0
Description:	Implementation of the LoTTA (Local Trimmed Taylor Approximation) model described in "Bayesian Regression Discontinuity Design with Unknown Cutoff" by Kowalska, van de Wiel, van der Pas (2024) <doi:10.48550/arXiv.2406.11585>.
License:	GPL-2 | GPL-3 [expanded from: GPL (≥ 2)]
Encoding:	UTF-8
RoxygenNote:	7.3.2
Imports:	stats, bayestestR, utils
Depends:	R (≥ 4.4.0), runjags, ggplot2, ggpubr
NeedsCompilation:	no
Packaged:	2025-07-07 15:29:27 UTC; juliakowalska
Author:	Julia Kowalska [aut, cre, cph]
Maintainer:	Julia Kowalska <j.m.kowalska@vu.nl>
Repository:	CRAN
Date/Publication:	2025-07-11 12:20:02 UTC

Function that evaluates the binary outcome function in a domain x, given the coefficients

Description

Function that evaluates the binary outcome function in a domain x, given the coefficients

Usage

BIN_outcome_function_sample(coef_s, x, d_x, s_x)

Arguments

Value

list with the values of the function

Function that splits the data into bins and computes the average in each bin

Description

Function that splits the data into bins and computes the average in each bin

Usage

Bin_data(y, x, c = NULL, binsize = 0.2)

Arguments

Value

list with elements y_bin: list of average values of y in each bin, x_bin: list of middle points for each bin

Function that evaluates the continuous outcome function in a domain x, given the coefficients

Description

Function that evaluates the continuous outcome function in a domain x, given the coefficients

Usage

CONT_outcome_function_sample(coef_s, x, d_x, s_x, mu_y, sd_y)

Arguments

Value

list with the values of the function

function that samples initial values for fuzzy LoTTA model with a continuous prior and binary outcomes

Description

function that samples initial values for fuzzy LoTTA model with a continuous prior and binary outcomes

Usage

Initial_CONT_BIN(x, t, y, C_start, clb, cub, lb, ubr, ubl, s, jlb = 0.2)

Arguments

Value

list with initial parameters values for LoTTA model with continuous score and binary outcomes, and .RNG.seed value

function that samples initial values for fuzzy LoTTA model with a ontinuous prior and continuous outcomes

Description

function that samples initial values for fuzzy LoTTA model with a ontinuous prior and continuous outcomes

Usage

Initial_CONT_CONT(x, t, y, C_start, clb, cub, lb, ubr, ubl, s, jlb = 0.2)

Arguments

Value

list with initial parameters values for LoTTA model with continuous score and continuous outcomes, and .RNG.seed value

function that samples initial values for LoTTA with a discerete prior and binary ourcomes

Description

function that samples initial values for LoTTA with a discerete prior and binary ourcomes

Usage

Initial_DIS_BIN(x, t, y, Ct_start, cstart, grid, lb, ubr, ubl, s, jlb = 0.2)

Arguments

Value

list with initial parameters values for LoTTA model with discrete score and binary outcomes, and .RNG.seed value

function that samples initial values for fuzzy LoTTA model with a discrete prior and binary outcomes

Description

function that samples initial values for fuzzy LoTTA model with a discrete prior and binary outcomes

Usage

Initial_DIS_CONT(x, t, y, Ct_start, cstart, grid, lb, ubr, ubl, s, jlb = 0.2)

Arguments

Value

list with initial parameters values for LoTTA model with discrete score and continuous outcomes, and .RNG.seed value

function that samples initial values for fuzzy LoTTA model with a known cutoff and binary outcomes

Description

function that samples initial values for fuzzy LoTTA model with a known cutoff and binary outcomes

Usage

Initial_FUZZy_BIN(x, t, y, c, lb, ubr, ubl, s, jlb = 0.2)

Arguments

Value

list with initial parameters values for fuzzy LoTTA model with a known cutoff and binary outcomes, and .RNG.seed value

function that samples initial values for fuzzy LoTTA model with a known cutoff and continuous outcomes

Description

function that samples initial values for fuzzy LoTTA model with a known cutoff and continuous outcomes

Usage

Initial_FUZZy_CONT(x, t, y, c, lb, ubr, ubl, s, jlb = 0.2)

Arguments

Value

list with initial parameters values for LoTTA model with continuous score and continuous outcomes, and .RNG.seed value

function that samples initial values for sharp LoTTA model with binary outcomes

Description

function that samples initial values for sharp LoTTA model with binary outcomes

Usage

Initial_SHARP_BIN(x, y, c, lb, ubr, ubl, s)

Arguments

Value

list with initial parameters values for sharp LoTTA model with binary outcomes, and .RNG.seed value

function that samples initial values for sharp LoTTA model with continuous outcomes

Description

function that samples initial values for sharp LoTTA model with continuous outcomes

Usage

Initial_SHARP_CONT(x, y, c, lb, ubr, ubl, s)

Arguments

Value

list with initial parameters values for sharp LoTTA model with continuous outcomes, and .RNG.seed value

function that samples initial values for the treatment model with a continuous prior

Description

function that samples initial values for the treatment model with a continuous prior

Usage

Initial_treatment_CONT(x, t, C_start, clb, cub, lb, ubr, ubl, s, jlb = 0.2)

Arguments

Value

list with initial parameters values for LoTTA model with continuous score and .RNG.seed value

function that samples initial values for the treatment model with a discrete prior

Description

function that samples initial values for the treatment model with a discrete prior

Usage

Initial_treatment_DIS(x, t, Ct_start, cstart, grid, lb, ubr, ubl, s, jlb = 0.2)

Arguments

Value

list with initial parameters values for treatment model with discrete score and .RNG.seed value

function that samples initial values for the treatment model with a known cutoff

Description

function that samples initial values for the treatment model with a known cutoff

Usage

Initial_treatment_c(x, t, c, lb, ubr, ubl, s, jlb = 0.2)

Arguments

Value

list with initial parameters values for LoTTA model with continuous score and .RNG.seed value

LoTTA_fuzzy_BIN

Description

Function that fits LoTTA model to the fuzzy RD data with binary outcomes with an either known or unknown/suspected cutoff. It supports two types of priors on the cutoff location: a scaled beta distribution of the form beta(alpha,beta)(cub-clb)+clb and a discrete distribution with the support of the form cstart+grid i for i=0,...,(cend-cstart)/grid. The score does NOT have to be normalized beforehand. We recommend NOT to transform the data before imputing it into the function, except for initial trimming of the score which should be done beforehand. The further trimming for the sensitivity analysis can be done through the function, which ensures that the data is normalized before the trimming.

Usage

LoTTA_fuzzy_BIN(
  x,
  t,
  y,
  c_prior,
  jlb = 0.2,
  ci = 0.95,
  trimmed = NULL,
  outcome_prior = list(pr = 1e-04),
  n_min = 25,
  param = c("c", "j", "kl", "kr", "eff", "a0l", "a1l", "a2l", "a3l", "a0r", "a1r", "a2r",
    "a3r", "b1lt", "a1lt", "a2lt", "b2lt", "b1rt", "a1rt", "a2rt", "b2rt", "k1t", "k2t"),
  normalize = TRUE,
  n.chains = 4,
  burnin = 10000,
  sample = 1000,
  adapt = 500,
  inits = NULL,
  method = "parallel",
  seed = NULL,
  ...
)

Arguments

Value

The function returns the list with the elements:

• Effect_estimate: contains a list with MAP estimate and HDI of the treatment effect, the cutoff location (if unknown) and the discontinuity size in the treatment probability function (compliance rate at c) on the original, unnormalized scale;

• JAGS_output: contains output of the run.jags function for the normalized data if normalize=TRUE, based on this output mixing of the chains can be assessed;

• Samples: contains posterior samples of the treatment effect (eff), cutoff location (c) if unknown, and compliance rate (j);

• Normalized_data: contains a list of the normalized data (if normalized=TRUE) and the parameters used to normalize the data (see arg normalize);

• Priors: contains a list of the priors' parameters ;

• Inits contains the list of initial values and .RNG.seed value

Examples

# 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)
}

funB <- function(x) {
  y = x
  x2 = x[x >= 0]
  x1 = x[x < 0]
  y[x < 0] = 1 / (1 + exp(-2 * x1)) - 0.5 + 0.4
  y[x >= 0] = (log(x2 * 2 + 1) - 0.15 * x2^2) * 0.6 - 0.20 + 0.4
  return(y)
}

funB_sample_binary <- function(x) {
  y = x
  P = funB(x)
  for (j in 1:length(x)) {
    y[j] = sample(c(1, 0), 1, prob = c(P[j], 1 - P[j]))
  }
  return(y)
}

## Toy example - for the function check only! ##
N=100
x = sort(runif(N, -1, 1))
t = sample_prob55(x)
y = funB_sample_binary(x)
c_prior=0 # known cutoff c=0

# running LoTTA model on fuzzy RDD with binary outcomes;
out = LoTTA_fuzzy_BIN(x,t,y,c_prior,burnin = 50,sample = 50,adapt=10,n.chains=1
,method = 'simple',inits = NA)

## Use case example ##

  N=500
  x = sort(runif(N, -1, 1))
  t = sample_prob55(x)
  y = funB_sample_binary(x)
  # comment out to try different priors:
   c_prior=list(clb=-0.25,cub=0.25) # uniform prior on the interval [-0.25,0.25]
  # c_prior=list(cstart=-0.25,cend=0.25,grid=0.05) # uniform discrete prior
  # on -0.25, -0.2, ..., 0.25
  # c_prior=0 # known cutoff c=0

  # running LoTTA model on fuzzy RDD with binary outcomes and unknown cutoff;
  # cutoff = 0, compliance rate = 0.55, treatment effect = -0.3636364
  # remember to check convergence and adjust burnin, sample and adapt if needed
  out = LoTTA_fuzzy_BIN(x,t,y,c_prior,burnin = 10000,sample = 5000,adapt=1000)

  # print effect estimate:
  out$Effect_estimate
  # print JAGS output to asses convergence (the output is for normalized data)
  out$JAGS_output
  # plot posterior fit of the outcome function
  LoTTA_plot_outcome(out,nbins = 60)
  # plot posterior fit of the treatment probablity function
  LoTTA_plot_treatment(out,nbins = 60)
  # plot dependence of the treatment effect on the cutoff location
  LoTTA_plot_effect(out)

 

LoTTA_fuzzy_CONT

Description

Function that fits LoTTA model to the fuzzy RD data with continuous outcomes with an either known or unknown/suspected cutoff. It supports two types of priors on the cutoff location: a scaled beta distribution of the form beta(alpha,beta)(cub-clb)+clb and a discrete distribution with the support of the form cstart+grid i for i=0,...,(cend-cstart)/grid. The score does NOT have to be normalized beforehand. We recommend NOT to transform the data before imputing it into the function, except for initial trimming of the score which should be done beforehand. The further trimming for the sensitivity analysis can be done through the function, which ensures that the data is normalized before the trimming.

Usage

LoTTA_fuzzy_CONT(
  x,
  t,
  y,
  c_prior,
  jlb = 0.2,
  ci = 0.95,
  trimmed = NULL,
  outcome_prior = list(pr = 1e-04, shape = 0.01, scale = 0.01),
  n_min = 25,
  param = c("c", "j", "kl", "kr", "eff", "a0l", "a1l", "a2l", "a3l", "a0r", "a1r", "a2r",
    "a3r", "b1lt", "a1lt", "a2lt", "b2lt", "b1rt", "a1rt", "a2rt", "b2rt", "k1t", "k2t",
    "tau1r", "tau2r", "tau1l", "tau2l"),
  normalize = TRUE,
  n.chains = 4,
  burnin = 10000,
  sample = 1000,
  adapt = 500,
  inits = NULL,
  method = "parallel",
  seed = NULL,
  ...
)

Arguments

Value

The function returns the list with the elements:

• Effect_estimate: contains a list with MAP estimate and HDI of the treatment effect, the cutoff location (if unknown) and the discontinuity size in the treatment probability function (compliance rate at c) on the original, unnormalized scale;

• JAGS_output: contains output of the run.jags function for the normalized data if normalize=TRUE, based on this output mixing of the chains can be assessed;

• Samples: contains posterior samples of the treatment effect (eff), cutoff location (c) if unknown, and compliance rate (j);

• Normalized_data: contains a list of the normalized data (if normalized=TRUE) and the parameters used to normalize the data (see arg normalize);

• Priors: contains a list of the priors' parameters ;

• Inits contains the list of initial values and .RNG.seed value

Examples

# 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)
}

funB <- function(x) {
  y = x
  x2 = x[x >= 0]
  x1 = x[x < 0]
  y[x < 0] = 1 / (1 + exp(-2 * x1)) - 0.5 + 0.4
  y[x >= 0] = (log(x2 * 2 + 1) - 0.15 * x2^2) * 0.6 - 0.20 + 0.4
  return(y)
}

funB_sample <- function(x) {
  y = funB(x)+ rnorm(length(x), 0, 0.1)
  return(y)
}

## Toy example - for the function check only! ##
N=100
x = sort(runif(N, -1, 1))
t = sample_prob55(x)
y = funB_sample(x)
c_prior=0 # known cutoff c=0

# running LoTTA model on fuzzy RDD with continuous outcomes;
out = LoTTA_fuzzy_CONT(x,t,y,c_prior,burnin = 50,sample = 50,adapt=10,n.chains=1
,method = 'simple',inits = NA)

## Use case example ##

  N=500
  x = sort(runif(N, -1, 1))
  t = sample_prob55(x)
  y = funB_sample(x)
  # comment out to try different priors:
  c_prior=list(clb=-0.25,cub=0.25) # uniform prior on the interval [-0.25,0.25]
  # c_prior=list(cstart=-0.25,cend=0.25,grid=0.05) # uniform discrete prior
  # on -0.25, -0.2, ..., 0.25
  # c_prior=0 # known cutoff c=0

  # running LoTTA model on fuzzy RDD with continuous outcomes and unknown cutoff;
  # cutoff = 0, compliance rate = 0.55, treatment effect = -0.3636364
  # remember to check convergence and adjust burnin, sample and adapt if needed
  out = LoTTA_fuzzy_CONT(x,t,y,c_prior,burnin = 10000,sample = 5000,adapt=1000)

  # print effect estimate:
  out$Effect_estimate
  # print JAGS output to asses convergence (the output is for normalized data)
  out$JAGS_output
  # plot posterior fit of the outcome function
  LoTTA_plot_outcome(out)
  # plot posterior fit of the treatment probablity function
  LoTTA_plot_treatment(out,nbins = 60)
  # plot dependence of the treatment effect on the cutoff location
  LoTTA_plot_effect(out,nbins = 5)

LoTTA_plot_effect

Description

Function that visualizes the impact of the cutoff location on the treatment effect estimate. It plots too figures. The bottom figure depicts the posterior density of the cutoff location. The top figure depicts the box plot of the treatment effect given the cutoff point. If the prior on the cutoff location was discrete each box corresponds to a distinct cutoff point. If the prior was continuous each box correspond to an interval of cutoff values (the number of intervals can be changed through nbins).

Usage

LoTTA_plot_effect(
  LoTTA_posterior,
  nbins = 10,
  probs = c(0.025, 0.975),
  x_lab = "Cutoff location",
  y_lab1 = "Treatment effect",
  y_lab2 = "Density of cutoff location",
  title = "Cutoff location vs. Treatment effect",
  axis.text = element_text(family = "sans", size = 10),
  text = element_text(family = "serif"),
  plot.theme = theme_classic(base_size = 14),
  plot.title = element_text(hjust = 0.5),
  ...
)

Arguments

Value

ggplot2 object

Examples

# 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)
}

funB <- function(x) {
  y = x
  x2 = x[x >= 0]
  x1 = x[x < 0]
  y[x < 0] = 1 / (1 + exp(-2 * x1)) - 0.5 + 0.4
  y[x >= 0] = (log(x2 * 2 + 1) - 0.15 * x2^2) * 0.6 - 0.20 + 0.4
  return(y)
}

funB_sample <- function(x) {
  y = funB(x)+ rnorm(length(x), 0, 0.1)
  return(y)
}

## Use case example ##

  N=500
  x = sort(runif(N, -1, 1))
  t = sample_prob55(x)
  y = funB_sample(x)
  # comment out to try different priors:
   c_prior=list(clb=-0.25,cub=0.25) # uniform prior on the interval [-0.25,0.25]
  # c_prior=list(cstart=-0.25,cend=0.25,grid=0.05) # uniform discrete prior
  # on -0.25, -0.2, ..., 0.25
  # c_prior=0 # known cutoff c=0

  # running LoTTA model on fuzzy RDD with continuous outcomes and unknown cutoff;
  # cutoff = 0, compliance rate = 0.55, treatment effect = -0.3636364
  # remember to check convergence and adjust burnin, sample and adapt if needed
  out = LoTTA_fuzzy_CONT(x,t,y,c_prior,burnin = 10000,sample = 5000,adapt=1000)

  # print effect estimate:
  out$Effect_estimate
  # print JAGS output to asses convergence (the output is for normalized data)
  out$JAGS_output
  # plot posterior fit of the outcome function
  LoTTA_plot_outcome(out,nbins = 60)
  # plot posterior fit of the treatment probablity function
  LoTTA_plot_treatment(out,nbins = 60)
  # plot dependence of the treatment effect on the cutoff location
  LoTTA_plot_effect(out)

Function that visualizes the impact of the cutoff location on the treatment effect estimate. It plots too figures. The bottom figure depicts the posterior density of the cutoff location. The top figure depicts the box plot of the treatment effect given the cutoff point. If the prior on the cutoff location was discrete each box corresponds to a distinct cutoff point. If the prior was continuous each box correspond to an interval of cutoff values (the number of intervals can be changed through nbins).

Description

Function that visualizes the impact of the cutoff location on the treatment effect estimate. It plots too figures. The bottom figure depicts the posterior density of the cutoff location. The top figure depicts the box plot of the treatment effect given the cutoff point. If the prior on the cutoff location was discrete each box corresponds to a distinct cutoff point. If the prior was continuous each box correspond to an interval of cutoff values (the number of intervals can be changed through nbins).

Usage

LoTTA_plot_effect_CONT(
  LoTTA_posterior,
  nbins = 10,
  probs = c(0.025, 0.975),
  x_lab = "Cutoff location",
  y_lab1 = "Treatment effect",
  y_lab2 = "Density of cutoff location",
  title = "Cutoff location vs. Treatment effect",
  axis.text = element_text(family = "sans", size = 10),
  text = element_text(family = "serif"),
  plot.theme = theme_classic(base_size = 14),
  plot.title = element_text(hjust = 0.5),
  ...
)

Arguments

Value

ggplot2 object

Function that visualizes the impact of the cutoff location on the treatment effect estimate. It plots too figures. The bottom figure depicts the posterior density of the cutoff location. The top figure depicts the box plot of the treatment effect given the cutoff point. If the prior on the cutoff location was discrete each box corresponds to a distinct cutoff point. If the prior was continuous each box correspond to an interval of cutoff values (the number of intervals can be changed through nbins).

Description

Function that visualizes the impact of the cutoff location on the treatment effect estimate. It plots too figures. The bottom figure depicts the posterior density of the cutoff location. The top figure depicts the box plot of the treatment effect given the cutoff point. If the prior on the cutoff location was discrete each box corresponds to a distinct cutoff point. If the prior was continuous each box correspond to an interval of cutoff values (the number of intervals can be changed through nbins).

Usage

LoTTA_plot_effect_DIS(
  LoTTA_posterior,
  probs = c(0.025, 0.975),
  x_lab = "Cutoff location",
  y_lab1 = "Treatment effect",
  y_lab2 = "Density of cutoff location",
  title = "Cutoff location vs. Treatment effect",
  axis.text = element_text(family = "sans", size = 10),
  text = element_text(family = "serif"),
  plot.theme = theme_classic(base_size = 14),
  plot.title = element_text(hjust = 0.5),
  ...
)

Arguments

Value

ggplot2 object

LoTTA_plot_outcome

Description

Function that plots the median (or another quantile) of the LoTTA posterior outcome function along with the quanatile-based credible interval. The function is plotted on top of the binned input data. To bin the data, the score data is divided into bins of fixed length, then the average outcome in each bin is calculated. The average outcomes are plotted against the average values of the score in the corresponding bins. The data is binned separately on each side of the cutoff, the cutoff is marked on the plot with a dotted line. In case of an unknown cutoff, the MAP estimate is used.

Usage

LoTTA_plot_outcome(
  LoTTA_posterior,
  nbins = 100,
  probs = c(0.025, 0.5, 0.975),
  n_eval = 200,
  col_line = "#E69F00",
  size_line = 0.1,
  alpha_interval = 0.35,
  col_dots = "gray",
  size_dots = 3,
  alpha_dots = 0.6,
  col_cutoff = "black",
  title = "Outcome function",
  subtitle = NULL,
  y_lab = "",
  x_lab = expression(paste(italic("x"), " - score")),
  plot.theme = theme_classic(base_size = 14),
  legend.position = "none",
  name_legend = "Legend",
  labels_legend = "median outcome fun.",
  text = element_text(family = "serif"),
  legend.text = element_text(size = 14),
  plot.title = element_text(hjust = 0.5),
  plot.subtitle = element_text(hjust = 0.5),
  ...
)

Arguments

Value

ggplot2 object

Examples

# functions to generate the data

ilogit <- function(x) {
  return(1 / (1 + exp(-x)))
}

funB <- function(x) {
  y = x
  x2 = x[x >= 0]
  x1 = x[x < 0]
  y[x < 0] = 1 / (1 + exp(-2 * x1)) - 0.5 + 0.4
  y[x >= 0] = (log(x2 * 2 + 1) - 0.15 * x2^2) * 0.6 - 0.20 + 0.4
  return(y)
}

funB_sample <- function(x) {
  y = funB(x)+ rnorm(length(x), 0, 0.1)
  return(y)
}
## Toy example - for the function check only! ##
# data generation
N=100
set.seed(1234)
x = sort(runif(N, -1, 1))
y = funB_sample(x)
c = 0

# running LoTTA function on sharp RDD with continuous outcomes;
out = LoTTA_sharp_CONT(x, y, c,normalize=FALSE, burnin = 100, sample = 100, adapt = 100,
n.chains=1, seed = NULL,method = 'simple',inits = NA)
# plot the outcome
LoTTA_plot_outcome(out,n_eval = 100)

## Use case example ##
# data generation
  N=500 # try different dataset size
  x = sort(runif(N, -1, 1))
  y = funB_sample(x)
  c = 0
  # plot the data
  plot(x,y)
  # running LoTTA function on sharp RDD with continuous outcomes;
  # cutoff = 0, treatment effect = -0.2
  # remember to check convergence and adjust burnin, sample and adapt if needed
  out = LoTTA_sharp_CONT(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)

Function that plots the median (or another quantile) of the LoTTA posterior treatment probability function along with the quanatile-based credible interval. The function is plotted on top of the binned input data. To bin the data, the score data is divided into bins of fixed length, then the proportion of treated is calculated in each bin. The proportions are plotted against the average values of the score in the corresponding bins. The data is binned separately on each side of the cutoff, the cutoff is marked on the plot with a dotted line. In case of an unknown cutoff, the MAP estimate is used.

Description

Function that plots the median (or another quantile) of the LoTTA posterior treatment probability function along with the quanatile-based credible interval. The function is plotted on top of the binned input data. To bin the data, the score data is divided into bins of fixed length, then the proportion of treated is calculated in each bin. The proportions are plotted against the average values of the score in the corresponding bins. The data is binned separately on each side of the cutoff, the cutoff is marked on the plot with a dotted line. In case of an unknown cutoff, the MAP estimate is used.

Usage

LoTTA_plot_treatment(
  LoTTA_posterior,
  nbins = 100,
  probs = c(0.025, 0.5, 0.975),
  n_eval = 200,
  col_line = "#E69F00",
  size_line = 0.1,
  alpha_interval = 0.35,
  col_dots = "gray",
  size_dots = 3,
  alpha_dots = 0.6,
  col_cutoff = "black",
  title = "Treatment probability function",
  subtitle = NULL,
  y_lab = "",
  x_lab = expression(paste(italic("x"), " - score")),
  plot.theme = theme_classic(base_size = 14),
  legend.position = "none",
  name_legend = "Legend",
  labels_legend = "median treatment prob. fun.",
  text = element_text(family = "serif"),
  legend.text = element_text(size = 14),
  plot.title = element_text(hjust = 0.5),
  plot.subtitle = element_text(hjust = 0.5),
  ...
)

Arguments

Value

ggplot2 object

Examples

# 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! ##
N=100
x = sort(runif(N, -1, 1))
t = sample_prob55(x)
c_prior=0 # known cutoff

# running LoTTA treatment-only model;
out = LoTTA_treatment(x,t,c_prior,burnin = 50, sample = 50, adapt = 10,n.chains=1
                      ,method = 'simple',inits = NA)
# plot posterior fit of the treatment probablity function
LoTTA_plot_treatment(out,nbins = 60)

## Use case example ##

  N=500
  x = sort(runif(N, -1, 1))
  t = sample_prob55(x)

  # comment out to try different priors:
   c_prior=list(clb=-0.25,cub=0.25) # uniform prior on the interval [-0.25,0.25]
  #  c_prior=list(cstart=-0.25,cend=0.25,grid=0.05) # uniform discrete prior
  # on -0.25, -0.2, ..., 0.25
  # c_prior=0 # known cutoff c=0

  # running LoTTA treatment-only model;
  # cutoff = 0, compliance rate = 0.55
  # remember to check convergence and adjust burnin, sample and adapt if needed
  out = LoTTA_treatment(x,t,c_prior,burnin = 10000,sample = 5000,adapt=1000)

  # print effect estimate:
  out$Effect_estimate
  # print JAGS output to asses convergence (the output is for normalized data)
  out$JAGS_output
  # plot posterior fit of the treatment probablity function
  LoTTA_plot_treatment(out,nbins = 60)

LoTTA_sharp_BIN

Description

Function that fits LoTTA model to the sharp RD data with binary outcomes. The score does NOT have to be normalized beforehand. We recommend NOT to transform the data before imputing it into the function, except for initial trimming of the score which should be done beforehand. The further trimming for the sensitivity analysis can be done through the function, which ensures that the data is normalized before the trimming.

Usage

LoTTA_sharp_BIN(
  x,
  y,
  c,
  ci = 0.95,
  trimmed = NULL,
  outcome_prior = list(pr = 1e-04),
  n_min = 25,
  param = c("eff", "a0l", "a1l", "a2l", "a3l", "a0r", "a1r", "a2r", "a3r", "kl", "kr"),
  normalize = TRUE,
  n.chains = 4,
  burnin = 5000,
  sample = 5000,
  adapt = 1000,
  inits = NULL,
  method = "parallel",
  seed = NULL,
  ...
)

Arguments

Value

The function returns the list with the elements:

• Effect_estimate: contains a list with MAP estimate and HDI of the treatment effect on the original, unnormalized scale;

• JAGS_output: contains output of the run.jags function for the normalized data if normalize=TRUE, based on this output mixing of the chains can be assessed;

• Samples: contains posterior samples of the treatment effect (eff);

• Normalized_data: contains a list of the normalized data (if normalized=TRUE) and the parameters used to normalize the data (see arg normalize);

• Priors: contains a list of the outcome prior parameters ;

• Inits contains the list of initial values and .RNG.seed value

Examples

# 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 ##

  # 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)

LoTTA_sharp_CONT

Description

Function that fits LoTTA model to the sharp RD data with continuous outcomes. The data does NOT have to be normalized beforehand. We recommend NOT to transform the data before imputing it into the function, except for initial trimming which should be done beforehand. The further trimming for the sensitivity analysis can be done through the function, which ensures that the data is normalized before the trimming.

Usage

LoTTA_sharp_CONT(
  x,
  y,
  c,
  ci = 0.95,
  trimmed = NULL,
  outcome_prior = list(pr = 1e-04, shape = 0.01, scale = 0.01),
  n_min = 25,
  param = c("eff", "a0l", "a1l", "a2l", "a3l", "a0r", "a1r", "a2r", "a3r", "tau1r",
    "tau2r", "tau1l", "tau2l", "kl", "kr"),
  normalize = TRUE,
  n.chains = 4,
  burnin = 10000,
  sample = 5000,
  adapt = 1000,
  inits = NULL,
  method = "parallel",
  seed = NULL,
  ...
)

Arguments

Value

The function returns the list with the elements:

• Effect_estimate: contains a list with MAP estimate and HDI of the treatment effect on the original, unnormalized scale;

• JAGS_output: contains output of the run.jags function for the normalized data if normalize=TRUE, based on this output mixing of the chains can be assessed;

• Samples: contains posterior samples of the treatment effect (eff);

• Normalized_data: contains a list of the normalized data (if normalized=TRUE) and the parameters used to normalize the data (see arg normalize);

• Priors: contains a list of the outcome prior parameters ;

• Inits contains the list of initial values and .RNG.seed value

Examples

# functions to generate the data

ilogit <- function(x) {
  return(1 / (1 + exp(-x)))
}

funB <- function(x) {
  y = x
  x2 = x[x >= 0]
  x1 = x[x < 0]
  y[x < 0] = 1 / (1 + exp(-2 * x1)) - 0.5 + 0.4
  y[x >= 0] = (log(x2 * 2 + 1) - 0.15 * x2^2) * 0.6 - 0.20 + 0.4
  return(y)
}

funB_sample <- function(x) {
  y = funB(x)+ rnorm(length(x), 0, 0.1)
  return(y)
}

## Toy example - for the function check only! ##
# data generation
N=100
set.seed(1234)
x = sort(runif(N, -1, 1))
y = funB_sample(x)
c = 0

# running LoTTA function on sharp RDD with continuous outcomes;
out = LoTTA_sharp_CONT(x, y, c,normalize=FALSE, burnin = 50, sample = 50, adapt = 10,
n.chains=1, seed = NULL,method = 'simple',inits = NA)

## Use case example ##
# data generation
  N=500 # try different dataset size
  x = sort(runif(N, -1, 1))
  y = funB_sample(x)
  c = 0
  # plot the data
  plot(x,y)
  # running LoTTA function on sharp RDD with continuous outcomes;
  # cutoff = 0, treatment effect = -0.2
  # remember to check convergence and adjust burnin, sample and adapt if needed
  out = LoTTA_sharp_CONT(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)

LoTTA_treatment

Description

Function that fits LoTTA treatment model to the fuzzy RD treatment data with an either known or unknown/suspected cutoff. It supports two types of priors on the cutoff location: a scaled beta distribution of the form beta(alpha,beta)(cub-clb)+clb and a discrete distribution with the support of the form cstart+grid i for i=0,...,(cend-cstart)/grid. The score does NOT have to be normalized beforehand. We recommend NOT to transform the data before imputing it into the function, except for initial trimming of the score which should be done beforehand. The further trimming for the sensitivity analysis can be done through the function, which ensures that the data is normalized before the trimming.

Usage

LoTTA_treatment(
  x,
  t,
  c_prior,
  jlb = 0.2,
  ci = 0.95,
  trimmed = NULL,
  n_min = 25,
  param = c("c", "j", "b1lt", "a1lt", "a2lt", "b2lt", "b1rt", "a1rt", "a2rt", "b2rt",
    "k1t", "k2t"),
  normalize = TRUE,
  n.chains = 4,
  burnin = 5000,
  sample = 1000,
  adapt = 500,
  inits = NULL,
  method = "parallel",
  seed = NULL,
  ...
)

Arguments

Value

The function returns the list with the elements:

• Effect_estimate: contains a list with MAP estimate and HDI of the cutoff location (if unknown) and the discontinuity size in the treatment probability function (compliance rate at c) on the original, unnormalized scale;

• JAGS_output: contains output of the run.jags function for the normalized data if normalize=TRUE, based on this output mixing of the chains can be assessed;

• Samples: contains posterior samples of the cutoff location (c) if unknown, and compliance rate (j);

• Normalized_data: contains a list of the normalized data (if normalized=TRUE) and the parameters used to normalize the data (see arg normalize);

• Priors: contains a list of the priors' parameters ;

• Inits contains the list of initial values and .RNG.seed value

Examples

# 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! ##
N=100
x = sort(runif(N, -1, 1))
t = sample_prob55(x)
c_prior=0 # known cutoff

# running LoTTA treatment-only model;
out = LoTTA_treatment(x,t,c_prior,burnin = 50, sample = 50, adapt = 10,n.chains=1
                      ,method = 'simple',inits = NA)

## Use case example ##

  N=500
  x = sort(runif(N, -1, 1))
  t = sample_prob55(x)

  # comment out to try different priors:
   c_prior=list(clb=-0.25,cub=0.25) # uniform prior on the interval [-0.25,0.25]
  #  c_prior=list(cstart=-0.25,cend=0.25,grid=0.05) # uniform discrete prior
  # on -0.25, -0.2, ..., 0.25
  # c_prior=0 # known cutoff c=0

  # running LoTTA treatment-only model;
  # cutoff = 0, compliance rate = 0.55
  # remember to check convergence and adjust burnin, sample and adapt if needed
  out = LoTTA_treatment(x,t,c_prior,burnin = 10000,sample = 5000,adapt=1000)

  # print effect estimate:
  out$Effect_estimate
  # print JAGS output to asses convergence (the output is for normalized data)
  out$JAGS_output
  # plot posterior fit of the treatment probablity function
  LoTTA_plot_treatment(out,nbins = 60)

function that finds maximum widow size to searxh for a cutoff

Description

function that finds maximum widow size to searxh for a cutoff

Usage

bounds(x, ns = 25)

Arguments

Value

list with ubl - minimum value of the window's left boundary point, ubr - maximum value of the window's right boundary point

inverse logit function

Description

inverse logit function

Usage

invlogit(x)

Arguments

Value

value of inverse logit function at x

logit function

Description

logit function

Usage

logit(x)

Arguments

Value

value of logit function at x

normalize continuous score function

Description

normalize continuous score function

Usage

normalize_cont_x(x, trimmed = NULL, normalize = TRUE)

Arguments

Value

list with x - normalized score, d - parameter used to shift the score

normalize continuous outcome function

Description

normalize continuous outcome function

Usage

normalize_cont_y(y, x, trimmed = NULL, normalize = TRUE)

Arguments

Value

list with y - normalized outcome, mu - shift parameter, sd - scaling parameter

normalize discrete score function

Description

normalize discrete score function

Usage

normalize_dis_x(x, trimmed = NULL, normalize = TRUE)

Arguments

Value

list with x - normalized score, s - parameter used to scale the score

function that searches for initial parameters of outcome function to initiate the sampler

Description

function that searches for initial parameters of outcome function to initiate the sampler

Usage

optimal_k(x, y, c, lb, ubr, ubl)

Arguments

Value

• a list with kl and kr and cubic functions parameters

function that searches for initial parameters of binary outcome function to initiate the sampler

Description

function that searches for initial parameters of binary outcome function to initiate the sampler

Usage

optimal_k_bin(x, y, c, lb, ubr, ubl)

Arguments

Value

• a list with kl and kr and cubic functions parameters

Function that plots the median (or another quantile) of the posterior function with binary outcome along with the quanatile-based credible interval. The function is plotted on top of the binned input data. To bin the data, the score data is divided into bins of fixed length, then the average outcome in each bin is calculated. The average outcomes are plotted against the average values of the score in the corresponding bins.

Description

Function that plots the median (or another quantile) of the posterior function with binary outcome along with the quanatile-based credible interval. The function is plotted on top of the binned input data. To bin the data, the score data is divided into bins of fixed length, then the average outcome in each bin is calculated. The average outcomes are plotted against the average values of the score in the corresponding bins.

Usage

plot_outcome_BIN(
  LoTTA_posterior,
  nbins = 100,
  probs = c(0.025, 0.5, 0.975),
  n_eval = 200,
  col_line = "#E69F00",
  size_line = 0.1,
  alpha_interval = 0.35,
  col_dots = "gray",
  size_dots = 3,
  alpha_dots = 0.6,
  col_cutoff = "black",
  title = "Outcome function",
  subtitle = NULL,
  y_lab = "",
  x_lab = expression(paste(italic("x"), " - score")),
  plot.theme = theme_classic(base_size = 14),
  legend.position = "none",
  name_legend = "Legend",
  labels_legend = "median outcome fun.",
  text = element_text(family = "serif"),
  legend.text = element_text(size = 14),
  plot.title = element_text(hjust = 0.5),
  plot.subtitle = element_text(hjust = 0.5),
  ...
)

Arguments

Value

ggplot2 object

Function that plots the median (or another quantile) of the posterior function of a continous outcome along with the quanatile-based credible interval. The function is plotted on top of the binned input data. To bin the data, the score data is divided into bins of fixed length, then the average outcome in each bin is calculated. The average outcomes are plotted against the average values of the score in the corresponding bins.

Description

Function that plots the median (or another quantile) of the posterior function of a continous outcome along with the quanatile-based credible interval. The function is plotted on top of the binned input data. To bin the data, the score data is divided into bins of fixed length, then the average outcome in each bin is calculated. The average outcomes are plotted against the average values of the score in the corresponding bins.

Usage

plot_outcome_CONT(
  LoTTA_posterior,
  nbins = 100,
  probs = c(0.025, 0.5, 0.975),
  n_eval = 200,
  col_line = "#E69F00",
  size_line = 0.1,
  alpha_interval = 0.35,
  col_dots = "gray",
  size_dots = 3,
  alpha_dots = 0.6,
  col_cutoff = "black",
  title = "Outcome function",
  subtitle = NULL,
  y_lab = "",
  x_lab = expression(paste(italic("x"), " - score")),
  plot.theme = theme_classic(base_size = 14),
  legend.position = "none",
  name_legend = "Legend",
  labels_legend = "median outcome fun.",
  text = element_text(family = "serif"),
  legend.text = element_text(size = 14),
  plot.title = element_text(hjust = 0.5),
  plot.subtitle = element_text(hjust = 0.5),
  ...
)

Arguments

Value

ggplot2 object

function that checks the type of a prior and whether it is correct

Description

function that checks the type of a prior and whether it is correct

Usage

read_prior(prior, minx, maxx, ubl, ubr, lb)

Arguments

Value

list with type of prior and prior parameters

Function that evaluates the treatment probability function in a domain x, given the coefficients

Description

Function that evaluates the treatment probability function in a domain x, given the coefficients

Usage

treatment_function_sample(coef_s, x, d_x, s_x)

Arguments

Value

list with the values of the function

Binary outcomes for trimmed score

Description

Binary outcomes for trimmed score

Usage

trim_dis_y(y, x, trimmed = NULL)

Arguments

Value

list with y - outcomes for trimmed score