使用带有 ggplot 的两个不同向量创建重叠直方图

Posted 2023-02-16

【中文标题】使用带有 ggplot 的两个不同向量创建重叠直方图

【英文标题】：Creating an overlap histogram using two different vectors with ggplot

【发布时间】：2021-11-17 20:35:32

【问题描述】：

我们的分析师对我们的数据进行了倾向得分分析。基本上，他使用国家、年龄和生物起始年份来“平衡”我们数据集中的女性和男性人口。他在两组（女性和男性）之间进行了重叠评估，并查看了线性化倾向得分以查看是否存在“良好”重叠。

数据集：

structure(list(gender = c(0, 1, 0, 0, 1, 1, 1, 1, 0, 0, 1, 0, 
1, 0, 1, 1, 1, 0, 0, 1), country = structure(c(1L, 2L, 2L, 3L, 
1L, 1L, 1L, 1L, 1L, 1L, 3L, 3L, 3L, 1L, 1L, 1L, 1L, 1L, 1L, 1L
), .Label = c("CH", "CZ", "DK", "IS", "NL", "NO", "PT", "RO", 
"SE", "SF", "SI", "TR", "UK"), class = "factor"), age = c(39, 
37, 54, 33, 30, 62, 30, 48, 34, 40, 39, 41, 29, 31, 37, 27, 22, 
23, 21, 31), bio_drug_name = structure(c(1L, 1L, 4L, 3L, 1L, 
3L, 4L, 3L, 1L, 4L, 3L, 5L, 4L, 4L, 1L, 5L, 1L, 3L, 4L, 2L), .Label = c("adalimumab", 
"certolizumab", "etanercept", "golimumab", "infliximab"), class = "factor"), 
    bio_drug_start_year = c(2007, 2011, 2012, 2012, 2012, 2004, 
    2012, 2012, 2012, 2012, 2012, 2012, 2016, 2015, 2013, 2015, 
    2013, 2013, 2014, 2013), asdas_crp_cii_6month = c(1, 1, 0, 
    1, 0, 0, 1, 0, 1, 1, 0, 1, 1, 1, 0, 1, 1, 0, 1, 0), bio_drug_start_year_centered = c(-8, 
    -4, -3, -3, -3, -11, -3, -3, -3, -3, -3, -3, 1, 0, -2, 0, 
    -2, -2, -1, -2), age_std = structure(c(-0.211016383746095, 
    -0.375088510873223, 1.01952456970737, -0.70323276512748, 
    -0.949340955818173, 1.67581307821588, -0.949340955818173, 
    0.527308188325984, -0.621196701563916, -0.12898032018253, 
    -0.211016383746095, -0.046944256618966, -1.03137701938174, 
    -0.867304892254609, -0.375088510873223, -1.19544914650887, 
    -1.60562946432669, -1.52359340076312, -1.68766552789025, 
    -0.867304892254609), .Dim = c(20L, 1L)), ID = 1:20), na.action = structure(c(`111395` = 169L, 
`769107` = 2619L, `844107` = 2624L, `164325` = 2681L, `1011013` = 2728L, 
`114174` = 2763L, `116484` = 2778L, `231118` = 3058L), class = "omit"), row.names = c("463", 
"7729", "7756", "8306", "8324", "128", "8440", "8450", "8663", 
"8809", "8840", "8857", "9020", "9033", "9101", "9324", "9377", 
"9523", "9702", "9718"), class = "data.frame")

用于创建 PS 模型和计算男性和女性的线性化 PS 分数的代码

psmod = glm( gender ~ country + age_std + bio_drug_start_year_centered, family = 'binomial', data = dat)
psmod = step(psmod, scope = list(lower = ~country + age_std + bio_drug_start_year_centered, 
                                 upper = ~(country + age_std + bio_drug_start_year_centered)^2+
                                   poly(dat$age_std,degree=3)[,2] + poly(dat$age_std,degree=3)[,3] +
                                   poly(dat$bio_drug_start_year_centered,degree=3)[,2] +
                                   poly(dat$bio_drug_start_year_centered,degree=3)[,3]
),
direction='forward' )
summary(psmod)

# Predict ps-score
ps = predict(psmod, type= 'response')
lps = log(ps/(1-ps))

# Overlap assessment
par(mfrow=c(2,1))
min.lps = min(lps)
max.lps = max(lps)
hist(lps[dat$gender==0], breaks=50,main='male', xlab='Linearized ps-score', xlim=c(min.lps,max.lps))
hist(lps[dat$gender==1], breaks=50,main='female', xlab='Linearized ps-score', xlim=c(min.lps,max.lps))

这是image的输出

虽然这对他来说很好，但对于科学期刊来说还不够。我想使用 ggplot 创建一个漂亮的直方图并显示男性和女性之间的重叠。 here 上有一些很好的例子但是，由于线性化 PS 分数的长度不同，我不知道如何将其转换为数据集，然后在 ggplot 上使用。

【问题讨论】：

请注意您的模型没有运行，因为国家/地区没有超过一个级别

已更改 - 现在应该运行自己测试了

您也可以使用cobalt 包中的bal.plot() 来自动执行此操作。

【参考方案1】：

我无法运行提供的大部分代码，但如果问题是您要填充直方图的两个变量具有不同数量的值，那么类似以下的方法应该可以工作：

library(tidyverse)

score_a <- rnorm(n = 50, mean = 0, sd = 1)
score_b <- rnorm(n = 75, mean = 2, sd = 0.75)

# Basic plot:
ggplot() +
  # Add one histogram:
  geom_histogram(aes(score_a), color = "black", fill = "red", alpha = 0.7) +
  # Add second, which has a different number of values
  geom_histogram(aes(score_b), color = "black", fill = "blue", alpha = 0.7) +
  # Black and white theme
  theme_bw()
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

编辑：如果您想更好地控制 x 轴并根据您的值的最小值/最大值进行设置，它可能类似于以下示例。请注意，这里我使用了 round() 函数，因为我在示例中使用了这些值，但如果不需要四舍五入，您可以省略此函数和 labels = 或 breaks = seq(from = min_x, to = max_x, by = 0.5)。

# Labeling the x-axis based on the min/max might look like this:

# Define axis breaks & labels:
min_x <- min(c(score_a, score_b))
max_x <- max(c(score_a, score_b))

ggplot() +
  # Add one histogram:
  geom_histogram(aes(score_a), color = "black", fill = "red", alpha = 0.7) +
  # Add second, which has a different number of values
  geom_histogram(aes(score_b), color = "black", fill = "blue", alpha = 0.7) +
  # Black and white theme
  theme_bw() +
  scale_x_continuous(
    breaks = round(x = seq(from = min_x, to = max_x, by = 0.5),
                   digits = 1),
    labels = round(x = seq(from = min_x, to = max_x, by = 0.5),
                   digits = 1))
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

^{由reprex package (v2.0.0) 于 2021-09-24 创建}

【讨论】：

这是一个很好的解决方案。在他的直方图中，他使用了基于最小和最大 LPS 的 x 轴。是否也可以将其合并到您的代码中？

@CoinAtlas 查看编辑

很棒的作品，现在完美了！ @cactusoxbird

【参考方案2】：

这是一种解决方案。将lps绑定到原始数​​据，然后使用ggplot2绘制性别并填充颜色。

dat2 <- cbind(dat, lps)

library(ggplot2)
library(dplyr)
dat2 <- mutate(dat2, gender = as.character(gender)) 

ggplot(dat2)+
  geom_histogram(aes(x= lps, fill = gender), bins = 10)

【讨论】：

感谢您的努力。然而，这个答案对我来说并不理想，因为与男性组相比，女性组（一半）的观察结果较少，所以这个输出看起来好像有重叠但不成比例。在此处查看您的代码给我的输出：ibb.co/cbDNTWq。这是我从 cactusoxbird 得到的输出，我更喜欢这里：ibb.co/QjzzcB6