2015-09-06 122 views
1

我在尋找將中的正態分佈擬合添加到分組直方圖中最優雅的方式。我知道這個問題之前已經被問過很多次了,但是沒有一個建議的選項,比如this onethis one讓我覺得非常優雅,至少沒有,除非stat_function可以用於每個特定的數據子部分。R:在ggplot2中添加正態擬合到分組直方圖

將正態分佈擬合疊加到非分組直方圖上的一種相對優雅的方法是使用geom_smoothmethod="nls"(除了事實之外,它不是自啓動函數,起始值必須指定):

library(ggplot2) 
myhist = data.frame(size = 10:27, counts = c(1L, 3L, 5L, 6L, 9L, 14L, 13L, 23L, 31L, 40L, 42L, 22L, 14L, 7L, 4L, 2L, 2L, 1L)) 
ggplot(data=myhist, aes(x=size, y=counts)) + geom_point() + 
    geom_smooth(method="nls", formula = y ~ N * dnorm(x, m, s), se=F, 
       start=list(m=20, s=5, N=300)) 

enter image description here

我想知道是否雖然這種方法也可以用來添加正態分佈適合於分組直方圖作爲

library(devtools) 
install_github("tomwenseleers/easyGgplot2",type="source") 
library("easyGgplot2") # load weight data 
ggplot(weight,aes(x = weight)) + 
+  geom_histogram(aes(y = ..count.., colour=sex, fill=sex),alpha=0.5,position="identity") 

enter image description here

如果有可能定義的任何包我也想知道一個+ stat_distrfit()+ stat_normfit()爲ggplot2任何機會(與分組的可能性)? (我真的找不到任何東西,但這似乎是一個普通的任務,所以我只是想知道)

原因我希望代碼儘可能短是因爲這是一門課程,我想保持儘可能簡單...

PS geom_density不適合我的目標,我也想繪製計數/頻率而不是密度。我也想讓他們在同一個面板上,並避免使用facet_wrap

+0

看看[這篇文章](http://www.stackoverflow。COM /問題/ 25075428#25091231)。 – jlhoward

回答

2

喜歡這個?

## simulate your dataset - could not get easyGplot2 to load.... 
set.seed(1)  # for reproducible example 
weight <- data.frame(sex=c("Female","Male"), weight=rnorm(1000,mean=c(65,67),sd=1)) 

library(ggplot2) 
library(MASS)  # for fitdistr(...) 
get.params <- function(z) with(fitdistr(z,"normal"),estimate[1:2]) 
df <- aggregate(weight~sex, weight, get.params) 
df <- data.frame(sex=df[,1],df[,2]) 
x <- with(weight, seq(min(weight),max(weight),len=100)) 
gg <- data.frame(weight=rep(x,nrow(df)),df) 
gg$y <- with(gg,dnorm(x,mean,sd)) 
gg$y <- gg$y * aggregate(weight~sex, weight,length)$weight * diff(range(weight$weight))/30 

ggplot(weight,aes(x = weight, colour=sex)) + 
    geom_histogram(aes(y = ..count.., fill=sex), alpha=0.5,position="identity") + 
    geom_line(data=gg, aes(y=y)) 

我想 「優雅」 是在旁觀者的眼睛。使用stat_function(...)的問題是args=...列表無法使用aes(...)進行映射,因爲註釋中的帖子解釋了該列表。因此,您必須創建一個輔助數據框架(本例中爲gg),該數據具有適合的分佈的x值和y值,並使用geom_line(...)

上面的代碼在MASS包中使用fitdistr(...)來計算根據正態性假設(如果有意義,您可以使用不同的分佈)按性別分組的數據平均值和sd的最大似然估計值。然後通過將weight中的範圍除以100個增量創建x軸,並計算dnorm(x,...)的適當均值和sd。由於結果是密度,所以我們必須調整:

gg$y <- gg$y * aggregate(weight~sex, weight,length)$weight * diff(range(weight$weight))/30 

因爲您要將其與計數數據進行映射。請注意,這假定您使用geom_histogram中的默認分箱(將x中的範圍分成30等分增量)。最後,我們使用gg作爲圖層特定數據集添加對geom_line(...)的調用。

+0

非常感謝你 - 這是我一直在尋找的!仍然有點令人驚訝的是,stat_function()不能被映射 - 我真的沒有看到任何遲早不應該實現的內在原因......我將嘗試將其包裝在ggplot2.normhist()中,函數在我的easyGgplot2 fork來保存我的學生一些代碼... :-) –