# 基于CRAN安装R包，检测没有则安装 Installing R packages based on CRAN and installing them if they are not detectedp_list = c("ggplot2", "latex2exp", "readxl", "patchwork", "dplyr", "ggdag",           "dagitty","broom","ggsci","brms" )for(p in p_list){if (!requireNamespace(p)){install.packages(p)}    library(p, character.only = TRUE, quietly = TRUE, warn.conflicts = FALSE)}
# 加载R包 Loading R packagessuppressWarnings(suppressMessages(library(ggplot2)))suppressWarnings(suppressMessages(library(latex2exp)))suppressWarnings(suppressMessages(library(readxl)))suppressWarnings(suppressMessages(library(patchwork)))suppressWarnings(suppressMessages(library(dplyr)))suppressWarnings(suppressMessages(library(ggdag)))suppressWarnings(suppressMessages(library(dagitty)))suppressWarnings(suppressMessages(library(broom)))suppressWarnings(suppressMessages(library(ggsci)))suppressWarnings(suppressMessages(library(brms)))

# 读取数据# Load datafig3a <- read_excel("data/data2.xlsx", sheet = "Fig3a", skip = 1)
# 设置颜色# Set colormypal <- pal_npg("nrc", alpha = 0.7)(10)#mypal <- c("#4DBBD5", "#00A087", "#3C5488")mypal <- c("#d5cabd", "#318eab", "#845ec2")
# 绘图# plotp_Frequency <- ggplot(data=fig3a, aes(x=h2)) +  geom_histogram(aes(fill=type), bins = 100, alpha = 0.7, color = "white") +  scale_fill_manual(values = mypal,                    labels = c(TeX(r"(\textit{h}${^2}$ (       )"),                                 TeX(r"(\textit{h}${^2}$ ()"),                                 TeX(r"(\textit{h}${^2}$ (    )")),                    name = "") +  scale_x_continuous(expand = expansion(mult = c(0,0)),                     breaks = seq(0, 1, 0.2),                     limits = c(0, 1)) +  scale_y_continuous(expand = expansion(mult = c(0,0))) +  labs(y = "Counts",       x = TeX(r"(\textit{h}${^2}$)")) +  geom_vline(xintercept = c(0.27, 0.37, 0.62), linetype = "dashed", color = mypal) +  theme_bw() +  theme(panel.grid = element_blank(),        legend.position = c(0.8, 0.8),        legend.background = element_blank(),        legend.box.background = element_rect(color = "black"),        axis.text = element_text(size = 12),        axis.title = element_text(size = 14, face = "bold"),        plot.margin = margin(10, 10, 10, 10)) +  geom_text(data = data.frame(x = c(0.27, 0.37, 0.62), y = 80, label = c("0.27", "0.37", "0.62")),            aes(x = x, y = y, label = label),            hjust = -0.1, size = 4, color = "black")
# 保存为PDF# Save as PDFggsave(filename = "results/Frequency_distribution_histogram1.pdf", plot = p_Frequency, width = 10, height = 6)

# 载入数据# Load datadata <- read.table("data/data_hisogram.txt", header = TRUE, sep = "\t", comment.char = "")
# 绘制上下两轴直方图# Draw a histogram of the upper and lower axes=p_Frequency2 <- ggplot(data,aes(x = value, fill = variable))+   geom_histogram(data = filter(data, variable == "group1"),                 aes(y = ..count..),                 binwidth = 100    #aes(fill=type),#直方图函数(histgram)                 #bins = 100    )+  geom_histogram(data = filter(data, variable == "group2"),                 aes(y = ..count..*(-1)),                 binwidth = 100    #aes(fill=type),#直方图函数(histgram)                 #bins = 100    )+  scale_x_continuous(limits = c(0,1100),                                          breaks = c(0,250,500,750,1000),                                          expand = c(0,0))+    scale_y_continuous(limits = c(-3.5,5.5),                                          breaks = seq(-4,6,1),                               labels = abs(seq(-4,6,1)))+  theme_classic(base_size = 15)+    theme(panel.border = element_rect(size = 1,fill='transparent'),                legend.position = 'none',          axis.text = element_text(colour = 'black'))+    geom_vline(xintercept = 700,linetype=2,cex=1)+    scale_fill_manual(values = c('grey60','#2AC643'))#p_Frequency2# 保存为PDF# Save as PDFggsave(filename = "results/Frequency_upper_lower1.pdf", width = 10, height = 6)

# 绘制上下多组频率分布直方图# Draw multiple groups of frequency distribution histograms above and belownets <- read.csv("data/mosquito_nets.csv")mosquito_dag <- dagify(  malaria_risk ~ net + income + health + temperature,  net ~ income + health + temperature,  health ~ income,  exposure = "net",  outcome = "malaria_risk",  coords = list(x = c(malaria_risk = 7, net = 3, income = 4, health = 5, temperature = 6),                y = c(malaria_risk = 2, net = 2, income = 3, health = 1, temperature = 3)),  labels = c(malaria_risk = "Risk of malaria", net = "Mosquito net", income = "Income",             health = "Health", temperature = "Nighttime temperatures",             resistance = "Insecticide resistance"))adjustmentSets(mosquito_dag)#> { health, income, temperature }
# 创建使用混杂因素预测治疗状况的模型# Create model that predicts treatment status using confoundersmodel_treatment_freq <- glm(net ~ income + temperature + health,                            data = nets,                            family = binomial(link = "logit"))
# 计算治疗权重的倾向得分和逆概率# Calculate propensity scores and inverse probability of treatment weightsnets_with_weights <- augment(model_treatment_freq, nets,                             type.predict = "response")nets_with_weights <- as.data.frame(nets_with_weights)colnames(nets_with_weights)[colnames(nets_with_weights) == ".fitted"] <- "propensity"nets_with_weights = nets_with_weights %>%  mutate(iptw = (net_num / propensity) + ((1 - net_num) / (1 - propensity)))
# 评估治疗对结果的影响# Estimates the effect of treatment on outcomemodel_outcome_freq <- lm(malaria_risk ~ net,                         data = nets_with_weights,                         weights = iptw)
# `net` 的系数应为 -10ish# Coefficient for `net` should be -10ishtidy(model_outcome_freq)# 自定义颜色方案 (Customize color palette)isfahan_colors <- c("#4DBBD5B2", "#3C5488B2")
# 绘图# Plotp_Frequency3 <- ggplot() +   geom_histogram(data = filter(nets_with_weights, net_num == 1),                  bins = 50, aes(x = propensity, weight = iptw),                  fill = colorspace::lighten(isfahan_colors[1], 0.35)) +  geom_histogram(data = filter(nets_with_weights, net_num == 0),                  bins = 50, aes(x = propensity, weight = iptw, y = -after_stat(count)),                 fill = colorspace::lighten(isfahan_colors[2], 0.35)) +  geom_histogram(data = filter(nets_with_weights, net_num == 1),                  bins = 50, aes(x = propensity),                  fill = isfahan_colors[1]) +   geom_histogram(data = filter(nets_with_weights, net_num == 0),                  bins = 50, aes(x = propensity, y = -after_stat(count)),                 fill = isfahan_colors[2]) +  annotate(geom = "label", x = 0.8, y = 70, label = "Treated (actual)",            fill = isfahan_colors[1], color = "white", hjust = 1, size = 4) +  annotate(geom = "label", x = 0.8, y = 90, label = "Treated (IPTW pseudo-population)",            fill = colorspace::lighten(isfahan_colors[1], 0.35), color = "white", hjust = 1, size = 4) +  annotate(geom = "label", x = 0.8, y = -60, label = "Untreated (actual)",            fill = isfahan_colors[2], color = "white", hjust = 1, size = 4) +  annotate(geom = "label", x = 0.8, y = -80, label = "Untreated (IPTW pseudo-population)",            fill = colorspace::lighten(isfahan_colors[2], 0.35), color = "white", hjust = 1, size = 4) +  geom_hline(yintercept = 0, color = "white", linewidth = 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_classic(base_size = 15)+    theme(panel.border = element_rect(size = 1,fill='transparent'),                legend.position = 'none',          axis.text = element_text(colour = 'black'))
# 保存为PDF# Save as PDFggsave(filename = "results/Frequency_upper_lower2.pdf", width = 10, height = 6)
# 组合图library(cowplot)width = 89height = 59p0 = plot_grid(p_Frequency2, p_Frequency3, labels = c("A", "B"), ncol = 2)ggsave("results/Frequency_distribution_histogram2.pdf", p0, width = width * 3, height = height * 2, units = "mm")