Plotting fitted glm output using ggplot2 for interactions

I am trying to plot my model output for two different models using ggplot2. When I visualize my plots they look exactly the same despite that they are fitted.

The only difference between the two models is the addition of a variable "age" to the second model (and all two-way interactions that include "age").

The first model is as follows:

ab_bi_f= glm(LRS_Bin ~ as.factor(DispersalFate) + MastN + as.factor(DispersalFate)*MastN, data = female, family = binomial(link = logit))

The second model is as follows:

rel_bi_f = glm(LRS_Bin ~ as.factor(DispersalFate) + Age + MastN + as.factor(DispersalFate)*MastN + Age*MastN +as.factor(DispersalFate)*Age, data = female, family = binomial(link = logit))

The code I use to fit and visualize the first models is:

fitted1<- function (fit) {
    ggplot(data=fit$ab_bi_f, aes(x=MastN, y=LRS_Bin, colour=factor(DispersalFate))) + 
         stat_smooth() + 
         ggtitle("Binary absolute LRS model - female") + 
         labs(x="Total masts in female's lifespan", y="Mean fitted binary absolute LRS", colour="Dispersal type")+
         scale_x_continuous(breaks=c(0,1,2))+
         theme_classic()+annotate("text", x = 0, y = 0.89, label = "a",size=6)+
         coord_cartesian(ylim=c(0.1, 0.9))+
         scale_y_continuous(breaks=c(0.1,0.3,0.5,0.7,0.9))+
         scale_color_manual(values=c("#000000", "#999999"))+
         theme(legend.position=c(0.3,0.85),
         plot.title = element_text(size = 10),
         legend.title=element_text(size=10),
         axis.title=element_text(size=8), 
         legend.key = element_rect(size = 0.1),
         legend.key.size = unit(0.5, "cm"),
         legend.direction="horizontal")
}
plot1<-fitted1(glm(MastN~LRS_Bin))

And, the code I used to visualize the second model is:

fitted2<- function (fit) {
    ggplot(data=fit$rel_bi_f, aes(x=MastN, y=LRS_Bin, colour=factor(DispersalFate))) + 
         stat_smooth() + 
         ggtitle("Binary absolute LRS model - female") + 
         labs(x="Total masts in female's lifespan", y="Mean fitted binary absolute LRS", colour="Dispersal type")+
         scale_x_continuous(breaks=c(0,1,2))+
         theme_classic()+annotate("text", x = 0, y = 0.89, label = "a",size=6)+
         coord_cartesian(ylim=c(0.1, 0.9))+
         scale_y_continuous(breaks=c(0.1,0.3,0.5,0.7,0.9))+
         scale_color_manual(values=c("#000000", "#999999"))+
         theme(legend.position=c(0.3,0.85),
         plot.title = element_text(size = 10),
         legend.title=element_text(size=10),
         axis.title=element_text(size=8), 
         legend.key = element_rect(size = 0.1),
         legend.key.size = unit(0.5, "cm"),
         legend.direction="horizontal")
}
plot2<-fitted2(glm(MastN~LRS_Bin))

This is the resulting plot:

The two plots look exactly the same!

I can't figure out if I am making a mistake in my code or if the two models (despite the second model being different from the first) result in the same output...

My data and code can be found here.

Solution

There are some errors in the code you are using. I am not sure how you are getting the output you have attached. However, if you just want to visualize the fitted values. You can do it as follows:-

ab_bi_f= glm(LRS_Bin ~ as.factor(DispersalFate) + MastN + as.factor(DispersalFate)*MastN, data = female, family = binomial(link = logit))

rel_bi_f = glm(LRS_Bin ~ as.factor(DispersalFate) + Age + MastN + as.factor(DispersalFate)*MastN + Age*MastN +as.factor(DispersalFate)*Age, data = female, family = binomial(link = logit))

female$fit1 <- predict(ab_bi_f, type = 'response')
female$fit2 <- predict(rel_bi_f, type = 'response')

plot1 <- ggplot(data=female, aes(x=MastN, y=fit1, colour=factor(DispersalFate))) + 
        geom_point() + 
        stat_smooth() + 
        ggtitle("Binary absolute LRS model - female") + 
        labs(x="Total masts in female's lifespan", y="Mean fitted binary absolute LRS", colour="Dispersal type")+
        scale_x_continuous(breaks=c(0,1,2))+
        theme_classic()+annotate("text", x = 0, y = 0.89, label = "a",size=6)+
        coord_cartesian(ylim=c(0.1, 0.9))+
        scale_y_continuous(breaks=c(0.1,0.3,0.5,0.7,0.9))+
        scale_color_manual(values=c("#000000", "#999999"))+
        theme(legend.position=c(0.3,0.85),
        plot.title = element_text(size = 10),
        legend.title=element_text(size=10),
        axis.title=element_text(size=8), 
        legend.key = element_rect(size = 0.1),
        legend.key.size = unit(0.5, "cm"),
        legend.direction="horizontal")

plot2 <- ggplot(data=female, aes(x=MastN, y=fit2, colour=factor(DispersalFate))) + 
        geom_point() + 
        stat_smooth() + 
        ggtitle("Binary absolute LRS model - female") + 
        labs(x="Total masts in female's lifespan", y="Mean fitted binary absolute LRS", colour="Dispersal type")+
        scale_x_continuous(breaks=c(0,1,2))+
        theme_classic()+annotate("text", x = 0, y = 0.89, label = "a",size=6)+
        coord_cartesian(ylim=c(0.1, 0.9))+
        scale_y_continuous(breaks=c(0.1,0.3,0.5,0.7,0.9))+
        scale_color_manual(values=c("#000000", "#999999"))+
        theme(legend.position=c(0.3,0.85),
        plot.title = element_text(size = 10),
        legend.title=element_text(size=10),
        axis.title=element_text(size=8), 
        legend.key = element_rect(size = 0.1),
        legend.key.size = unit(0.5, "cm"),
        legend.direction="horizontal")

grid.arrange(plot1, plot2)