推文内容来自于链接
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]
image.png
这里用到的配色包是 https://github.com/BlakeRMills/MetBrewer 这个用到的都是博物馆里的油画的配色，挺有意思的，大家可以试试
使用ggplot2作图
这里频率分布直方图用到的是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
image.png
频率分布直方图
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")
image.png
如果要倒过来加一个负号就可以了
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")
image.png
添加文本注释
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)
image.png
对细节的一些调整
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"))
image.png
欢迎大家关注我的公众号
小明的数据分析笔记本
小明的数据分析笔记本 公众号 主要分享：1、R语言和python做数据分析和数据可视化的简单小例子；2、园艺植物相关转录组学、基因组学、群体遗传学文献阅读笔记；3、生物信息学入门学习资料及自己的学习笔记！
