Reputation: 1823
I would like to visualized the output of multiple comparisons for GLM in ggplot. In my example:
#Simulation data set
set.seed(123)
d<-NULL
N<-54
d$Production <- rgamma(N,10)
d$Feature <- ifelse(d$Production >7 & d$Production<10, c("A"),ifelse(d$Production>11,
c("B"), c("C")))
#d$Temp<-rnorm(N,20,5)
d$Temp<-rep(1:3,18)
d<-as.data.frame(d)
#
# Gamma glm model
mG<- glm(Production~ Feature + Temp, family= Gamma, data = d)
summary(mG)
anova(mG,test="Chi")
#If for example a have interaction between Temp and Feature
d$BV <- interaction(d$Feature, d$Temp)
# GLM fit
mI <- glm(Production ~ -1 + BV, data = d, family = "Gamma")
#summary(glht(mI, linfct = mcp(BV = "Tukey")))
cld(glht(mI, linfct = mcp(BV = "Tukey")))
# Prediction values
pred.data = data.frame(
Feature<-d$Feature,
Temp<-d$Temp
)
pred.data$Production = predict(mG, newdata=pred.data, type="response")
# Select means
library(dplyr)
d2<-d %>%
group_by(Feature, Temp) %>%
summarize(Production = mean(Production, na.rm = TRUE))
#Final plot
library("ggplot2")
ggplot(d2, aes(Temp, Production, colour = Feature)) +
geom_point() +
stat_smooth(method = "glm", formula = y ~ x, family = Gamma)
#
Now I've like to visualized in ggplot the letters of cld(glht(mI, linfct = mcp(BV = "Tukey"))) for an output:
This is possible? Thanks
Upvotes: 0
Views: 587
Reputation: 2399
Try this:
library(tidyverse)
res <- cld(glht(mI, linfct = mcp(BV = "Tukey")))
d2 <-
res$mcletters$Letters %>%
tibble(
reg = names(.),
lett = .
) %>%
mutate(
Feature = substr(reg, 1, 1),
Temp = substr(reg, 3, 3) %>% as.integer()
) %>%
left_join(d2) %>%
mutate(lett_pos = if_else(Feature == 'C', Production - .5, Production + .5))
ggplot(d2, aes(Temp, Production, colour = Feature)) +
geom_point() +
stat_smooth(method = "glm", formula = y ~ x) +
geom_text(aes(y = lett_pos, label = lett), color = 'black')
Upvotes: 2