使用 geom_smooth 将 glm 拟合到分数

这篇文章有点相关这个帖子 https://stackoverflow.com/questions/62974766/removing-alpha-transparency-from-ggplot-legend-and-setting-x-axis-labels-to-be-f/62986181#62986181.

这里我有xy分组数据，其中y是分数：

library(dplyr)
library(ggplot2)
library(ggpmisc)

set.seed(1)
df1 <- data.frame(value = c(0.8,0.5,0.4,0.2,0.5,0.6,0.5,0.48,0.52),
                 age = rep(c("d2","d4","d45"),3),
                 group = c("A","A","A","B","B","B","C","C","C")) %>%
  dplyr::mutate(time = as.integer(age)) %>%
  dplyr::arrange(group,time) %>%
  dplyr::mutate(group_age=paste0(group,"_",age))
df1$group_age <- factor(df1$group_age,levels=unique(df1$group_age))

我想要实现的是绘制df1作为条形图，如下所示：

ggplot(df1,aes(x=group_age,y=value,fill=age)) +
  geom_bar(stat='identity')

但我想适合每一个group a binomial glm with a logit link function，它估计这些分数如何受到影响time.

假设每个观察值有 100 个age (time） 每一个group:

df2 <- do.call(rbind,lapply(1:nrow(df1),function(i){
  data.frame(age=df1$age[i],group=df1$group[i],time=df1$time[i],group_age=df1$group_age[i],value=c(rep(T,100*df1$value[i]),rep(F,100*(1-df1$value[i]))))
}))

然后glm对于每个group (e.g., group A) is:

glm(value ~ time, dplyr::filter(df2, group == "A"), family = binomial(link='logit'))

所以我想添加到估计上方的图中regression slopes对于每个group以及它们相应的p-values（类似于我为连续所做的事情df$value in 这个帖子 https://stackoverflow.com/questions/62974766/removing-alpha-transparency-from-ggplot-legend-and-setting-x-axis-labels-to-be-f/62986181#62986181).

我认为使用：

ggplot(df1,aes(x=group_age,y=value,fill=age)) +
  geom_bar(stat='identity') +
  geom_smooth(data=df2,mapping=aes(x=group_age,y=value,group=group),color="black",method='glm',method.args=list(family=binomial(link='logit')),size=1,se=T) +
  stat_poly_eq(aes(label=stat(p.value.label)),formula=my_formula,parse=T,npcx="center",npcy="bottom") +
  scale_x_log10(name="Age",labels=levels(df$age),breaks=1:length(levels(df$age))) +
  facet_wrap(~group) + theme_minimal()

可以工作，但我收到错误：

Error in Math.factor(x, base) : ‘log’ not meaningful for factors

知道如何做对吗？

我相信这会有所帮助：

library(tidyverse)
library(broom)

df2$value <- as.numeric(df2$value)

#Estimate coefs
dfmodel <- df2 %>% group_by(group) %>%
  do(fitmodel = glm(value ~ time, data = .,family = binomial(link='logit')))
#Extract coeffs
dfCoef = tidy(dfmodel, fitmodel)
#Create labels
dfCoef %>% filter(term=='(Intercept)') %>% mutate(Label=paste0(round(estimate,3),'(p=',round(p.value,3),')'),
                                                  group_age=paste0(group,'_','d4')) %>%
  select(c(group,Label,group_age)) -> Labels
#Values
df2 %>% group_by(group,group_age) %>% summarise(value=sum(value)) %>% ungroup() %>%
  group_by(group) %>% filter(value==max(value)) %>% select(-group_age) -> values
#Combine
Labels %>% left_join(values) -> Labels
Labels %>% mutate(age=NA) -> Labels
#Plot
ggplot(df2,aes(x=group_age,y=value,fill=age)) +
  geom_text(data=Labels,aes(x=group_age,y=value,label=Label),fontface='bold')+
  geom_bar(stat='identity')+
  facet_wrap(.~group,scales='free')