-
Notifications
You must be signed in to change notification settings - Fork 9
/
Copy pathrdmcplot_illustration.R
63 lines (48 loc) · 2.32 KB
/
rdmcplot_illustration.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
###################################################################
# rdmulti: analysis of RD designs with multiple cutoffs or scores
# RDMCPLOT illustration file
# Authors: Matias Cattaneo, Rocio Titiunik, Gonzalo Vazquez-Bare
###################################################################
rm(list = ls())
library(rdmulti)
library(ggplot2)
###################################################################
# Load data and generate plot variables (omitting plot)
###################################################################
data <- read.csv('simdata_multic.csv')
Y <- data$y
X <- data$x
C <- data$c
aux <- rdmcplot(Y,X,C,ci=95,nodraw=TRUE)
Xmean <- aux$Xmean
Ymean <- aux$Ymean
X0 <- aux$X0
X1 <- aux$X1
Yhat0 <- aux$Yhat0
Yhat1 <- aux$Yhat1
CI_l <- aux$CI_l
CI_r <- aux$CI_r
###################################################################
# Replicate default plot from rdmcplot output directly
###################################################################
aux$rdmc_plot
###################################################################
# Replicate default plot from rdmcplot by hand
###################################################################
# First cutoff: C=33
rdmc_plot <- ggplot() + theme_bw() + labs(x="Running variable", y="Outcome")
rdmc_plot <- rdmc_plot + geom_point(aes(x=Xmean[,1],y=Ymean[,1]),col="darkblue",shape=1,na.rm=TRUE) +
geom_line(aes(x=X0[,1],y=Yhat0[,1]),col="darkblue",linetype=1,na.rm=TRUE) +
geom_line(aes(x=X1[,1],y=Yhat1[,1]),col="darkblue",linetype=1,na.rm=TRUE) +
geom_vline(xintercept=33,col="darkblue",linetype=1)
rdmc_plot
# Add second cutoff: C=66
rdmc_plot <- rdmc_plot + geom_point(aes(x=Xmean[,2],y=Ymean[,2]),col="darkred",shape=1,na.rm=TRUE) +
geom_line(aes(x=X0[,2],y=Yhat0[,2]),col="darkred",linetype=1,na.rm=TRUE) +
geom_line(aes(x=X1[,2],y=Yhat1[,2]),col="darkred",linetype=1,na.rm=TRUE) +
geom_vline(xintercept=66,col="darkred",linetype=1)
rdmc_plot
# Add confidence intervals for both cutoffs
rdmc_plot <- rdmc_plot + geom_errorbar(aes(x=Xmean[,1], ymin=CI_l[,1], ymax=CI_r[,1]),col="darkblue",linetype=1) +
geom_errorbar(aes(x=Xmean[,2], ymin=CI_l[,2], ymax=CI_r[,2]),col="darkred",linetype=1)
rdmc_plot