在 R 中使用 ggplot 直方图而不是 hist 函数

Posted

技术标签:

【中文标题】在 R 中使用 ggplot 直方图而不是 hist 函数【英文标题】:Using ggplot histogram instead of hist function in R 【发布时间】:2022-01-10 13:27:51 【问题描述】:

我在 R 中使用名为 BetaMixture 的包来拟合数据向量的 beta 分布的混合。将输出提供给hist,该hist 使用混合模型组件生成良好的直方图:

# Install and load the libraries

#install.packages("BetaModels")
library(BetaModels)

# Create a vector, fit mixture models and plot the histogram

vec <- c(rbeta(700, 5, 2), rbeta(300, 1, 10))
model <- BetaMixture(vec,2)
h <- hist(model, breaks = 35)

到目前为止一切顺利。现在我如何在ggplot中得到这个?我检查了h 对象,但这与model 对象没有什么不同。它们完全相同。我什至不知道这个hist 是如何为这个班级工作的。除了@datavec,它从model 中提取了什么来生成此图?

【问题讨论】:

【参考方案1】:

您可以使用getMethod("hist", "BetaMixture") 获取BetaMixed 对象的hist 函数。 您可以在下面找到此函数到“ggplot2world”的简单翻译。

myhist <- function (x, ...) 
    .local <- function (x, mixcols = 1:7, breaks=25, ...) 
    
        df1 <- data.frame(x=x@datavec)
        p <- ggplot(data=df1, aes(x=x)) + 
             geom_histogram(aes(y=..density..), bins=breaks, alpha=0.5, fill="gray50", color="black")
        while (length(mixcols) < ncol(x@mle)) mixcols <- c(mixcols, 
            mixcols)
        xv <- seq(0, 1, length = 502)[1:501]
        for (J in 1:ncol(x@mle)) 
            y <- x@phi[J] * dbeta(xv, x@mle[1, J], x@mle[2, J])
            df2 <- data.frame(xv, y)
            p <- p + geom_line(data=df2, aes(xv, y), size=1, col=mixcols[J])
        
        p <- p + theme_bw()
        invisible(p)
    
    .local(x, ...)


library(ggplot2)
# Now p is a ggplot2 object.
p <- myhist(model, breaks=35)
print(p)

【讨论】:

这很漂亮,谢谢! 另外,您使用 getMethod 读取本地函数的解释令人大开眼界。谢谢!【参考方案2】:

BetaMixture 返回的对象是一个 S4 类对象,有 2 个感兴趣的插槽。

Slot Z 返回属于每个分布的每个数据点的概率矩阵。 所以在前 6 行中,所有点都属于第二个分布。

head(model@Z)
#             [,1]      [,2]
#[1,] 1.354527e-04 0.9998645
#[2,] 4.463074e-03 0.9955369
#[3,] 1.551999e-03 0.9984480
#[4,] 1.642579e-03 0.9983574
#[5,] 1.437047e-09 1.0000000
#[6,] 9.911427e-04 0.9990089

mle返回参数的最大似然估计。

现在使用这些值来创建向量的 data.frame 和参数的 data.frame。

df1 <- data.frame(vec)
df1$component <- factor(apply(model@Z, 1, which.max))
colors <- as.integer(levels(df1$component))

params <- as.data.frame(t(model@mle))
names(params) <- c("shape1", "shape2")

绘制数据。

library(ggplot2)

g <- ggplot(df1, aes(x = vec, group = component)) +
  geom_histogram(aes(y = ..density..),
                 bins = 35, fill = "grey", color = "grey40")

for(i in 1:nrow(params))
  sh1 <- params$shape1[i]
  sh2 <- params$shape2[i]
  g <- g + stat_function(
    fun = dbeta,
    args = list(shape1 = sh1, shape2 = sh2),
    color = colors[i]
  )

suppressWarnings(print(g + theme_bw()))

【讨论】:

非常感谢您的详细解释,尤其是 Z 插槽。我在想它们是实际的组件值而不是概率。不知道有专门为这个包开发的本地hist方法。这是一个非常好的包,我应该说!此外,您使用概率的解决方案非常好。 但是,对于一般的拟合问题,您还需要使用权重参数,该参数由模型中的@phi 给出。否则,组件会爆炸。 @user2167741 是的,没错,我考虑过使用权重。但在这种情况下,它可以在没有它们的情况下工作,所以我选择不复杂化。

以上是关于在 R 中使用 ggplot 直方图而不是 hist 函数的主要内容,如果未能解决你的问题,请参考以下文章

使用ggplot将直方图转换为r中的小提琴图

在R中标记对数刻度显示

我无法使用 ggplot 在 r 中填充直方图

如何根据 R 中的直方图桶对数据进行分类?

用 ggplot2 (R) 覆盖 KDE 和填充直方图

R语言使用ggplot2绘制带有边缘直方图的散点图实战