-
Notifications
You must be signed in to change notification settings - Fork 2
/
Copy pathrdpower_illustration.R
218 lines (141 loc) · 5.28 KB
/
rdpower_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
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
###################################################################
# rdpower: power calculations for RD designs
# Illustration file
# !version 2.2 20-Jun-2022
# Authors: Matias Cattaneo, Rocio Titiunik, Gonzalo Vazquez-Bare
###################################################################
## NOTE: if you are using rdrobust version 2020 or newer, the option
## masspoints="off" and stdvars="on" may be needed in order to replicate the
## results in the Stata journal article.
## For example, line 39:
## aux = rdpower(data=Z,tau=5)
## should be replaced by:
## aux = rdpower(data=Z,tau=5,masspoints="off",stdvars="on")
###################################################################
#################################################################
# Setup
#################################################################
rm(list = ls())
library(rdpower)
library(rdrobust)
data <- read.csv('rdpower_senate.csv')
Y <- data$demvoteshfor2
R <- data$demmv
names(data)
# Outcome and running variables
Z <- data[c('demvoteshfor2','demmv')]
# Covariates
covs <- data[c('population','dopen','dmidterm')]
# Cluster var
clustvar <- data$state
#################################################################
# RDPOWER
#################################################################
# rdpower against tau = 5
aux <- rdpower(data=Z,tau=5)
# rdpower with covariates
aux <- rdpower(data=Z,tau=5,covs=covs)
# rdpower with plot
aux <- rdpower(data=Z,tau=5,plot=TRUE)
# rdpower with rdrobust options
aux <- rdpower(data=Z,tau=5,h=c(16,18),b=c(18,20))
aux <- rdpower(data=Z,tau=5,kernel='uniform',cluster=data$state,vce='hc0')
aux <- rdpower(data=Z,tau=5,bwselect='certwo',vce='hc3',scaleregul=0,rho=1)
# rdpower with conventional inference
aux <- rdpower(data=Z,tau=5,all=TRUE)
# rdpower with user-specified bias and variance
rdr <- rdrobust(Y,R)
samph <- rdr$bws[1]
sampsi.l <- rdr$N_h_l
sampsi.r <- rdr$N_h_r
bias.l <- rdr$bias[1]/(rdr$bws[1]^2)
bias.r <- rdr$bias[2]/(rdr$bws[2]^2)
VL <- rdr$V_rb_l
VR <- rdr$V_rb_r
N <- sum(rdr$N)
Vl.rb <- N*rdr$bws[1]*VL[1,1]
Vr.rb <- N*rdr$bws[1]*VR[1,1]
aux <- rdpower(data=Z,tau=5,bias=c(bias.l,bias.r),variance=c(Vl.rb,Vr.rb),samph=samph,sampsi=c(sampsi.l,sampsi.r))
# rdpower manually increasing variance by 20%
rdr <- rdrobust(Y,R)
VL <- rdr$V_rb_l
VR <- rdr$V_rb_r
N = sum(rdr$N)
Vl.rb = N*rdr$bws[1]*VL[1,1]*1.2
Vr.rb = N*rdr$bws[1]*VR[1,1]*1.2
aux <- rdpower(data=Z,tau=5,variance=c(Vl.rb,Vr.rb))
# rdpower without data
aux1 <- rdpower(data=Z,tau=5)
aux <- rdpower(tau=5,
nsamples=c(aux1$N.l,aux1$Nh.l,aux1$N.r,aux1$Nh.r),
bias=c(aux1$bias.l,aux1$bias.r),
variance=c(aux1$Vl.rb,aux1$Vr.rb),
sampsi=c(aux1$sampsi.l,aux1$sampsi.r),
samph=c(aux1$samph.l,aux1$samph.r))
# comparing exp-post power across specifications
aux1 <- rdpower(data=Z,tau=5,p=1,h=20,plot=TRUE)
aux2 <- rdpower(data=Z,tau=5,p=2,h=20,plot=TRUE)
aux3 <- rdpower(data=Z,tau=5,p=1,plot=TRUE)
aux4 <- rdpower(data=Z,tau=5,p=2,plot=TRUE)
# rdpower with clustering
aux <- rdpower(data=Z,tau=5,cluster=clustvar)
#################################################################
## RDSAMPSI
#################################################################
# rdsampsi with tau = 5
aux <- rdsampsi(data=Z,tau=5)
# rdsampsi setting bandwidth and nratio with plot
aux <- rdsampsi(data=Z,tau=5,beta=.9,samph=c(18,19),nratio=.5,plot=TRUE)
# rdsampsi with conventional inference
aux <- rdsampsi(data=Z,tau=5,all=TRUE)
# rdsampsi vs rdpower
aux1 <- rdsampsi(data=Z,tau=5)
aux2 <- rdpower(data=Z,tau=5,sampsi=c(aux1$sampsi.h.l,aux1$sampsi.h.r))
# rdsampsi without data
aux1 <- rdsampsi(data=Z,tau=5)
aux2 <- rdsampsi(tau = 5,
nsamples = c(aux1$N.l,aux1$Nh.l,aux1$N.r,aux1$Nh.r),
bias = c(aux1$bias.l,aux1$bias.r),
variance = c(aux1$var.l,aux1$var.r),
samph = c(aux1$samph.l,aux1$samph.r),
init.cond = aux1$init.cond)
# comparing sample sizes across designs
aux1 <- rdsampsi(data=Z,tau=5,p=0,h=20,plot=TRUE)
aux2 <- rdsampsi(data=Z,tau=5,p=1,h=20,plot=TRUE)
aux3 <- rdsampsi(data=Z,tau=5,p=0,plot=TRUE)
aux4 <- rdsampsi(data=Z,tau=5,p=1,plot=TRUE)
# rdsampsi with clustering
aux <- rdsampsi(data=Z,tau=5,cluster=clustvar)
#################################################################
## RDMDE
#################################################################
# rdmde with default options
aux <- rdmde(data=Z)
# rdmde for a power of 0.75
aux <- rdmde(data=Z,beta=0.75)
# rdmde manually setting the bandwidths
aux <- rdmde(data=Z,samph=c(12,13))
# rdmde manually setting the sample size inside the bandwidths
aux <- rdmde(data=Z,sampsi=c(350,320))
# rdmde without data
rdr <- rdrobust(Y,R)
bias.l <- rdr$bias[1]/(rdr$bws[1]^2)
bias.r <- rdr$bias[2]/(rdr$bws[2]^2)
VL <- rdr$V_rb_l
VR <- rdr$V_rb_r
N <- sum(rdr$N)
Vl.rb <- N*rdr$bws[1]*VL[1,1]
Vr.rb <- Vl.rb*1.1
aux <- rdmde(nsamples=c(600,350,700,320),bias=c(bias.l,bias.r),var=c(Vl.rb,Vr.rb),samph=17,init.cond=4)
# compare rdmde and rdpower
aux1 <- rdmde(data=Z)
mde <- aux1$mde
aux2 <- rdpower(data=Z,tau=mde)
# compare rdmde and rdsampsi
aux1 <- rdmde(data=Z)
mde <- aux1$mde
nhl <- aux1$Nh.l
nhr <- aux1$Nh.r
aux2 <- rdsampsi(data=Z,tau=mde,nratio=nhr/(nhl+nhr))
# rdmde with clustering
aux <- rdmde(data=Z,cluster=clustvar)