Я хотел бы создать график в ggplot2, который бы совмещал два разных уровня в измененной модели пуассона GLM.Но мое представление доверительного интервала (IC - 95%) не работает.
#Packages
library(ggplot2)
library(dplyr)
library(tidyverse)
#Artificial data set
Consumption <- c(501, 502, 503, 504, 26, 27, 50, 56, 63, 60, 72, 93, 78, 43, 59, 70, 53, 80)
Gender <- gl(n = 3, k = 6, length = 3*6, labels = c("Male", "Female","hermaphrodite"), ordered = FALSE)
Income <- c(5010, 5020, 5030, 5040, 260, 270, 550, 560, 680, 690, 720, 550, 560, 680, 690, 720,500,512)
df3 <- data.frame(Consumption, Gender, Income)
df3
# GLM Regression
fm1 <- glm(Consumption~Gender+Income, data=df3, family=poisson)
summary(fm1)
# ANOVA
anova(fm1,test="Chi")
#Comparing Gender
sort(tapply(df3$Consumption,df3$Gender,mean))
#Pairwise comparison - Stepwise like
Gender2<-df3$Gender
levels(Gender2)
levels(Gender2)[2]<-"Fem_Her"
levels(Gender2)[3]<-"Fem_Her"
levels(Gender2)
fm2<-glm(Consumption~Gender2+Income, data=df3, family=poisson)
anova(fm1,fm2,test="Chi")
# 0.7824 Female/Hermaphrodite are equal
#Genders are different than I have one parameter for male and another for Female/Hermaphrodite
pred <- predict(fm2, type="response", se.fit = TRUE)
df3 = cbind(df3, pred = pred$fit)
df3 = cbind(df3, se = pred$se.fit)
df3 = cbind(df3, ucl=df3$pred + 1.96*df3$se)
df3 = cbind(df3, lcl=df3$pred - 1.96*df3$se)
df3 = cbind(df3, Gender2)
df<-df3 %>%
dplyr::group_by(Income, Gender2) %>%
dplyr::summarize(Consumption = mean(Consumption, na.rm = TRUE))
df<-as.data.frame(df)
#Plot
df3 %>%
tidyr::gather(type, value, Consumption) %>%
ggplot(mapping=aes(x=type, y=value, color = Gender2)) +
geom_smooth(mapping=aes(ymin = lcl, ymax = ucl), stat = "identity") +
geom_point(df,mapping=aes(x=Income, y=Consumption, color = Gender2)) +
geom_line(mapping=aes(x=Income, y=pred))
#
Мой окончательный вывод не имеет доверительного интервала:
Есть идеи, пожалуйста?