LoTTA_plot_treatment: Function that plots the median (or another quantile) of the...

In LoTTA: Bayesian Inference in Regression Discontinuity Designs

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

LoTTA_posterior	output of one of the LoTTA functions (LoTTA_fuzzy_CONT, LoTTA_fuzzy_BIN, LoTTA_treatment) with all the parameters sampled (the default option in those functions)
nbins	number of bins to aggregate the input data
probs	list of three quantiles, the first and the last one define the quanatile-based credible interval, the middle value defines the quantile of the posterior function to plot; by default the quantiles correspond to the median posterior function and 95% credible interval probs=c(0.025,0.5,0.975)
n_eval	n_eval*range(x) is the number of points at which each posterior function is evaluated, the higher number means slower computing time and a smoother plot; default n_eval=200
col_line	the color of the line and the band
size_line	thickness of the line
alpha_interval	alpha value of the band, lower values correspond to a more transparent color
col_dots	color of the dots that correspond to the binned data
size_dots	size of the dots that correspond to the binned data
alpha_dots	transparency of the dots that correspond to the binned data, lower values correspond to a more transparent color
col_cutoff	color of the dotted line at the cutoff
title	title of the plot
subtitle	subtitle of the plot
y_lab	label of the y-axis
x_lab	label of the x-axis
plot.theme	ggplot2 plot theme (see https://ggplot2.tidyverse.org/reference/ggtheme.html) possibly with additional arguments, it takes the default value plot.theme=theme_classic(base_size = 14),
legend.position	position of the legend, refer to ggplot2 manual for the possible values; by default legend is not printed legend.position='none'
name_legend	title of the legend
labels_legend	the label of the plotted function in the legend
text	can be any value that is accepted in the argument text in the theme function of ggplot2 package,refer to ggplot2 manual for the possible values; by default is changes font to a serif one text=element_text(family='serif')
legend.text	can be any value that is accepted in the argument legend.text in the theme function of ggplot2 package,refer to ggplot2 manual for the possible values; by default is changes font size to the legend to 14 legend.text=element_text(size = 14)
plot.title	can be any value that is accepted in the argument plot.title in the theme function of ggplot2 package,refer to ggplot2 manual for the possible values; by default it centers the plot title plot.title = element_text(hjust = 0.5)
plot.subtitle	can be any value that is accepted in the argument plot.subtitle in the theme function of ggplot2 package,refer to ggplot2 manual for the possible values; by default it centers the plot subtitle plot.title = element_text(hjust = 0.5)
...	other arguments of the theme function, refer to ggplot2 manual

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 documentation built on Aug. 8, 2025, 6:27 p.m.