问使用geom_smooth将glm拟合为分数
EN

Stack Overflow用户

提问于 2020-07-19 23:25:47

回答 2查看 169关注 0票数 2

这篇文章有点和this post有关。

这里我有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拟合一个带有logit link function的binomial glm，它估计这些片段如何受到time的影响。

假设我在每个age (time)中每个group都有100个观察值

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]))))
}))

则每个group (例如，group A)的glm为：

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

因此，我想在上面的图中添加每个group的估计regression slopes以及它们相应的p-value(类似于我在this post中对连续df$value所做的操作)。

我认为使用：

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

你知道怎么做才是对的吗？

ggplot2

bar-chart

glm

平台电商解决方案

为客户提供电商客服、电商图片系统识别及搜索的快速搭建方案

回答 2

Stack Overflow用户

发布于 2020-07-20 02:06:25

我相信这可能会有所帮助：

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

票数 3

Stack Overflow用户

发布于 2020-07-20 19:16:49

多亏了Pedro Aphalo，这几乎是一个完整的解决方案：

生成带有分数的data.frame (这里使用time作为整数，方法是删除age中的"d"，而不是使用time作为age的levels )：

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(gsub("d","",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扩大到每个group中的每个age 100个观察值，但将value指定为integer而不是binary

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(1,100*df1$value[i]),rep(0,100*(1-df1$value[i]))))
}))

现在使用geom_smooth和stat_fit_tidy绘制它

ggplot(df1,aes(x=time,y=value,group=group,fill=age)) +
  geom_bar(stat='identity') +
  geom_smooth(data=df2,mapping=aes(x=time,y=value,group=group),color="black",method='glm',method.args=list(family=binomial(link='logit'))) +
  stat_fit_tidy(data=df2,mapping=aes(x=time,y=value,group=group,label=sprintf("P = %.3g",stat(x_p.value))),method='glm',method.args=list(formula=y~x,family=binomial(link='logit')),parse=T,label.x="center",label.y="top") +
  scale_x_log10(name="Age",labels=levels(df2$age),breaks=unique(df2$time)) +
  facet_wrap(~group) + theme_minimal()