https://www.andrewheiss.com/blog/2021/12/18/bayesian-propensity-scores-weights/
这个博文里的内容还挺多的,我们只关注其中关于频率分布直方图的实现代码。
nets_with_weights<-read.csv("nets_with_weights.csv")
isfahan <- MetBrewer::met.brewer("Isfahan1")
length(isfahan)
isfahan[1]
这里用到的配色包是 https://github.com/BlakeRMills/MetBrewer 这个用到的都是博物馆里的油画的配色,挺有意思的,大家可以试试
这里频率分布直方图用到的是geom_histogram()
函数,这里的代码多了一个weight
参数,暂时没有想明白这个参数起到什么作用
还遇到一个新函数colorspace::lighten()
操作颜色,看帮助文档是是颜色更亮。做一个散点图试试效果
library(ggplot2)
library(patchwork)
p1<-ggplot()+
geom_point(aes(x=1,y=1),size=50,color="darkgreen")
p2<-ggplot()+
geom_point(aes(x=1,y=1),size=50,
color=colorspace::lighten("darkgreen",0.9))
p1+p2
ggplot() +
geom_histogram(data = filter(nets_with_weights, net_num == 1),
bins = 50,
aes(x = propensity, weight = iptw),
fill = colorspace::lighten(isfahan[2], 0.35),
color="white")
ggplot() +
geom_histogram(data = filter(nets_with_weights, net_num == 1),
bins = 50,
aes(x = propensity, weight = iptw),
fill = colorspace::lighten(isfahan[2], 0.35),
color="white")+
geom_histogram(data = filter(nets_with_weights, net_num == 0),
bins = 50, aes(x = propensity, weight = iptw,
y = -..count..),
fill = colorspace::lighten(isfahan[6], 0.35),
color="white")
ggplot() +
geom_histogram(data = filter(nets_with_weights, net_num == 1),
bins = 50,
aes(x = propensity, weight = iptw),
fill = colorspace::lighten(isfahan[2], 0.35),
color="white")+
geom_histogram(data = filter(nets_with_weights, net_num == 0),
bins = 50, aes(x = propensity, weight = iptw,
y = -..count..),
fill = colorspace::lighten(isfahan[6], 0.35),
color="white")+
geom_histogram(data = filter(nets_with_weights, net_num == 1),
bins = 50, aes(x = propensity),
fill = isfahan[2],color="white") +
geom_histogram(data = filter(nets_with_weights, net_num == 0),
bins = 50, aes(x = propensity, y = -..count..),
fill = isfahan[6],
color="white")+
annotate(geom = "label",
x = 0.8, y = 70,
label = "Treated (actual)",
fill = isfahan[2],
color = "white", hjust = 1) +
annotate(geom = "label", x = 0.8,
y = 90, label = "Treated (IPTW pseudo-population)",
fill = colorspace::lighten(isfahan[2], 0.35),
color = "white", hjust = 1) +
annotate(geom = "label", x = 0.8, y = -60,
label = "Untreated (actual)",
fill = isfahan[6],
color = "white", hjust = 1) +
annotate(geom = "label",
x = 0.8, y = -80,
label = "Untreated (IPTW pseudo-population)",
fill = colorspace::lighten(isfahan[6], 0.35),
color = "white", hjust = 1)
ggplot() +
geom_histogram(data = filter(nets_with_weights, net_num == 1),
bins = 50,
aes(x = propensity, weight = iptw),
fill = colorspace::lighten(isfahan[2], 0.35),
color="white")+
geom_histogram(data = filter(nets_with_weights, net_num == 0),
bins = 50, aes(x = propensity, weight = iptw,
y = -..count..),
fill = colorspace::lighten(isfahan[6], 0.35),
color="white")+
geom_histogram(data = filter(nets_with_weights, net_num == 1),
bins = 50, aes(x = propensity),
fill = isfahan[2],color="white") +
geom_histogram(data = filter(nets_with_weights, net_num == 0),
bins = 50, aes(x = propensity, y = -..count..),
fill = isfahan[6],
color="white")+
annotate(geom = "label",
x = 0.8, y = 70,
label = "Treated (actual)",
fill = isfahan[2],
color = "white", hjust = 1) +
annotate(geom = "label", x = 0.8,
y = 90, label = "Treated (IPTW pseudo-population)",
fill = colorspace::lighten(isfahan[2], 0.35),
color = "white", hjust = 1) +
annotate(geom = "label", x = 0.8, y = -60,
label = "Untreated (actual)",
fill = isfahan[6],
color = "white", hjust = 1) +
annotate(geom = "label",
x = 0.8, y = -80,
label = "Untreated (IPTW pseudo-population)",
fill = colorspace::lighten(isfahan[6], 0.35),
color = "white", hjust = 1) +
geom_hline(yintercept = 0, color = "white", size = 0.25) +
scale_y_continuous(label = abs) +
coord_cartesian(xlim = c(0.1, 0.8), ylim = c(-80, 100)) +
labs(x = "Propensity", y = "Count")+
theme_minimal() +
theme(panel.grid.minor = element_blank(),
plot.background = element_rect(fill = "white", color = NA),
plot.title = element_text(face = "bold"),
axis.title = element_text(face = "bold"),
strip.text = element_text(face = "bold", size = rel(0.8), hjust = 0),
strip.background = element_rect(fill = "grey80", color = NA),
legend.title = element_text(face = "bold"))
示例数据和代码大家可以自己到推文开头提到的链接去下载