R: добавить нормальные соответствия сгруппированным гистограммам в ggplot2

Я ищу наиболее элегантный способ наложить совпадения нормального распределения на сгруппированные гистограммы в ggplot2. Я знаю, что этот вопрос задавался много раз раньше, но ни один из предложенных вариантов, например этот или этот мне показался очень элегантным, по крайней мере, если только stat_function нельзя заставить работать с каждым конкретным подразделом данных.

Один относительно элегантный способ наложить соответствие нормального распределения на несгруппированную гистограмму, с которым я столкнулся, заключался в использовании geom_smooth и method="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)) 

введите описание изображения здесь

Однако мне было интересно, можно ли использовать этот подход для добавления подходов нормального распределения к сгруппированным гистограммам, как в

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")

введите описание изображения здесь

Мне также было интересно, есть ли какие-нибудь пакеты, которые могут случайно определить + stat_distrfit() или + stat_normfit() для ggplot2 (с возможностью группировки)? (На самом деле я ничего не мог найти, но это могло показаться достаточно распространенной задачей, поэтому мне просто было интересно)

Причина, по которой я хочу, чтобы код был как можно короче, заключается в том, что это для курса, и что я хочу, чтобы все было как можно проще ...

PS geom_density не подходит для моей цели, и я также хотел бы построить график количества / частоты в противоположность плотности. Я также хотел бы, чтобы они были на одной панели и не использовали facet_wrap


person Tom Wenseleers    schedule 06.09.2015    source источник
comment
Прочтите этот пост.   -  person jlhoward    schedule 06.09.2015


Ответы (1)


Нравится?

## 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(...), как объясняется в сообщении в комментариях. Поэтому вам нужно создать вспомогательный data.frame (gg в этом примере), который имеет значения x и y для подобранных распределений, и использовать geom_line(...).

В приведенном выше коде используется fitdistr(...) в пакете MASS для вычисления оценок максимального правдоподобия среднего и SD ваших данных, сгруппированных по полу, на основе предположения о нормальности (вы можете использовать другое распределение, если это имеет смысл). Затем он создает ось x, разделяя диапазон в weight на 100 приращений, и вычисляет dnorm(x,...) для соответствующего среднего и стандартного отклонения. Поскольку результатом является плотность, мы должны скорректировать ее, используя:

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

потому что вы хотите сопоставить это с данными подсчета. Обратите внимание, что здесь предполагается, что вы используете биннинг по умолчанию в geom_histogram (который делит диапазон по x на 30 равных приращений). Наконец, мы добавляем вызов geom_line(...), используя gg в качестве набора данных для конкретного слоя.

person jlhoward    schedule 06.09.2015
comment
Большое спасибо за это - да, это было то, что я искал! Тем не менее немного удивительно, что stat_function () не может быть отображена - я действительно не вижу никакой внутренней причины, по которой это невозможно реализовать рано или поздно ... Я попытаюсь обернуть это в ggplot2.normhist () в моей вилке easyGgplot2, чтобы сэкономить моим ученикам немного кода ... :-) - person Tom Wenseleers; 06.09.2015