-
Notifications
You must be signed in to change notification settings - Fork 0
/
post.R
201 lines (171 loc) · 8.62 KB
/
post.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
library(here)
library(tidyverse)
library(reshape2)
library(gtable)
library(grid)
library(gridExtra)
## BIAS
a <- read_delim("./results_test.txt",delim="\t");a
a <- a %>% filter(type=="bias"&!is.na(type));a
a$N<-a$n
nrow(a)
A<-c("PMT","NPT","PMF","NPF")
vv<-c(paste("ipw",A,sep=""),paste("reg",A,sep=""),
paste("aipw",A,sep=""),paste("tmle",A,sep=""))
res1 <- melt(a,measure.vars=vv)[,c("N","variable","value")]
head(res1);tail(res1)
res1$N<-as.factor(res1$N);res1$length<-str_length(res1$variable)
table(res1$length)
res1$estimator<-str_sub(res1$variable,1,res1$length-3)
res1$estimator<-factor(res1$estimator)
res1$estimator<-factor(res1$estimator,levels(res1$estimator)[c(3,2,1,4)])
res1$model_type<-str_sub(res1$variable,-3,-1)
res1$variable<-res1$length<-NULL
head(res1);tail(res1)
lngth<-nrow(res1[res1$N==100&res1$estimator=="reg"&res1$model_type=="PMT",])
summary(res1[res1$N==100&res1$estimator=="aipw"&res1$model_type=="PMT",]$value)
data.frame(res1 %>% group_by(N,estimator,model_type) %>% summarize(q=quantile(value,.9)))
res1 <- res1 %>% mutate(value=ifelse(abs(value)>60,0,value))
D <- res1 %>% group_by(N,estimator,model_type) %>% summarize(mv = mean(value),sd.mv=sd(value)/sqrt(lngth),M=n())
D <- D %>% ungroup(N,estimator,model_type)
D
facet_labs<-list(
"NPF"="Nonpar Misspec",
"NPT"="Nonpar Correct",
"PMF"="Par Misspec",
"PMT"="Par Correct"
)
type_labeller <- function(variable,value){
return(facet_labs[value])
}
ggplot(res1,aes(N,value,fill = estimator)) +
facet_grid(model_type ~ .,labeller = type_labeller) + theme_light() +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank()) +
labs(x = "Sample Size",y = "Bias",fill="Estimator\n") +
scale_fill_grey(labels = c("g Comp", "IPW", "AIPW","TMLE")) +
geom_boxplot(width=.5,position = position_dodge(.65)) +
geom_hline(yintercept = 0,color="red",linetype="dashed",size=.25)
ggplot(D,aes(as.factor(N),mv,group=estimator)) +
theme_light() + facet_grid(model_type ~ .,labeller = type_labeller) +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank()) +
labs(x = "Sample Size",y = "Bias",fill = "Estimator\n") +
geom_boxplot(width=.25,position = position_dodge(),outlier.colour = NA)
#geom_boxplot(aes(fill=estimator), width = 0.5,position = position_dodge(),stat="identity") +
#geom_bar(aes(fill=estimator), width = 0.5,position = position_dodge(),stat="identity") +
#geom_text(size = 5, position = position_dodge(width=.5)) +
#geom_errorbar(aes(ymin=mv-sd.mv, ymax=mv+sd.mv),width=.1,position = position_dodge(.5)) +
scale_fill_grey(labels = c("g Comp", "IPW", "AIPW","TMLE")) +
#coord_cartesian(ylim=c(-2,2)) +
geom_hline(yintercept = 0,color="black",
linetype="dashed",size=.25)
# ind_seq<-c("PMT","PMF","NPT","NPF")
# plotFunc<-function(ind){
# ggplot(D[D$model_type==ind,],aes(as.factor(N),mv,group=estimator,label=M)) +
# theme_light() + facet_grid(model_type ~ .,labeller = type_labeller) +
# theme(panel.grid.major = element_blank(),
# panel.grid.minor = element_blank()) +
# labs(x = "Sample Size",y = "Bias",fill = "Estimator\n") +
# geom_bar(aes(fill=estimator), width = 0.5,position = position_dodge(),stat="identity") +
# geom_text(size = 5, position = position_dodge(width=.5)) +
# #geom_errorbar(aes(ymin=mv-sd.mv, ymax=mv+sd.mv),width=.1,position = position_dodge(.5)) +
# scale_fill_grey(labels = c("g Comp", "IPW", "AIPW","TMLE")) +
# #coord_cartesian(ylim=c(-2,2)) +
# geom_hline(yintercept = 0,color="black",
# linetype="dashed",size=.25)
# }
# plot_object<-lapply(ind_seq,function(x) plotFunc(x))
#
# legend = gtable_filter(ggplotGrob(plot_object[[1]]), "guide-box")
# thm<-theme(legend.position="none",
# axis.title.x=element_blank(),
# axis.title.y=element_blank(),
# axis.ticks.x=element_blank(),
# axis.text.x=element_blank())
#
# grid.arrange(arrangeGrob(plot_object[[4]] + thm,
# plot_object[[3]] + thm,
# plot_object[[2]] + thm,
# plot_object[[1]] + theme(legend.position="none",axis.title.y=element_blank()),
# nrow = 4,
# left = textGrob("Bias", rot = 90, vjust = 1)),
# legend,
# widths=unit.c(unit(1, "npc") - legend$width, legend$width),
# ncol=2)
#
# D
# res1 %>% filter(N==100,estimator=="ipw",model_type=="NPF")
## MSE
a <- read_delim("./results_test.txt",delim="\t");a
a$N<-a$n
a <- a %>% filter(type=="rmse")
nrow(a)
a %>% group_by(N) %>% summarize(n=n())
A<-c("PMT","NPT","PMF","NPF")
vv<-c(paste("ipw",A,sep=""),paste("reg",A,sep=""),
paste("aipw",A,sep=""),paste("tmle",A,sep=""))
res1 <- melt(a,measure.vars=vv)[,c("N","variable","value")] ############# vv
head(res1);tail(res1)
res1$N<-as.factor(res1$N);res1$length<-str_length(res1$variable)
res1$estimator<-str_sub(res1$variable,1,res1$length-3)
res1$estimator<-factor(res1$estimator)
res1$estimator<-factor(res1$estimator,levels(res1$estimator)[c(3,2,1,4)])
res1$model_type<-str_sub(res1$variable,-3,-1)
res1$variable<-res1$length<-NULL
data.frame(res1 %>% group_by(N,estimator,model_type) %>% summarize(q=quantile(value,.9)))
res1 <- res1 %>% mutate(value=ifelse(value>1000,1000,value))
D <- res1 %>% group_by(N,estimator,model_type) %>% summarize(mv = mean(value))
D <- D %>% ungroup(N,estimator,model_type)
plotFunc<-function(ind){
ggplot(D[D$model_type==ind,],aes(as.factor(N),mv,group=estimator)) +
facet_grid(model_type ~ .,labeller = type_labeller) + theme_light() +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank()) +
labs(x = "Sample Size",y = "rMSE",shape = "Estimator\n") +
scale_y_log10() +
scale_shape(labels = c("g Comp", "IPW", "AIPW","TMLE")) +
geom_point(aes(shape=estimator),position=position_dodge(.5),size=2.5)
}
plot_object<-lapply(ind_seq,function(x) plotFunc(x))
legend = gtable_filter(ggplotGrob(plot_object[[1]]), "guide-box")
thm<-theme(legend.position="none",
axis.title.x=element_blank(),
axis.title.y=element_blank(),
axis.ticks.x=element_blank(),
axis.text.x=element_blank())
grid.arrange(arrangeGrob(plot_object[[4]] + thm,
plot_object[[3]] + thm,
plot_object[[2]] + thm,
plot_object[[1]] + theme(legend.position="none",axis.title.y=element_blank()),
nrow = 4,
left = textGrob("rMSE", rot = 90, vjust = 1)),
legend,
widths=unit.c(unit(1, "npc") - legend$width, legend$width),
ncol=2)
###### TABLE
a <- read_delim("./results_ipw.txt",delim="\t");a
a <- a %>% filter(!is.na(type));a
a$N<-factor(a$N,levels=c("100","200","600","1200"))
ggplot(subset(a,type=="bias"),aes(ipwPMT,group=N,color=N)) +
geom_density() + xlab("Bias")
ggplot(subset(a,type=="bias"),aes(ipwNPT,group=N,color=N)) +
geom_density() + xlab("Bias")
ggplot(subset(a,type=="bias"),aes(regNPT,group=N,color=N)) +
geom_density() + xlab("Bias")
a %>% filter(type=="bias") %>% group_by(N) %>% summarise(M=n(),ipwPT=mean(ipwPMT),ipwNPT=mean(ipwNPT),
ipwPF=mean(ipwPMF),ipwNPF=mean(ipwNPF),
regPT=mean(regPMT),regNPT=mean(regNPT),
regPF=mean(regPMF),regNPF=mean(regNPF))
a %>% filter(type=="rmse") %>% group_by(N) %>% summarise(M=n(),ipwPT=mean(ipwPMT),ipwNPT=mean(ipwNPT),
ipwPF=mean(ipwPMF),ipwNPF=mean(ipwNPF),
regPT=mean(regPMT),regNPT=mean(regNPT),
regPF=mean(regPMF),regNPF=mean(regNPF))
a %>% filter(type=="cov") %>% group_by(N) %>% summarise(M=n(),ipwPT=mean(ipwPMT),ipwNPT=mean(ipwNPT),
ipwPF=mean(ipwPMF),ipwNPF=mean(ipwNPF),
regPT=mean(regPMT),regNPT=mean(regNPT),
regPF=mean(regPMF),regNPF=mean(regNPF))
a %>% filter(type=="width") %>% group_by(N) %>% summarise(M=n(),ipwPT=mean(ipwPMT),ipwNPT=mean(ipwNPT),
ipwPF=mean(ipwPMF),ipwNPF=mean(ipwNPF),
regPT=mean(regPMT),regNPT=mean(regNPT),
regPF=mean(regPMF),regNPF=mean(regNPF))