Using dplyr and pipes for logistic regression plotting
Here is what I would do leaning fully into the tidyverse style to get a nice "pipelined" set of steps that are easy to wrap up into functions.

Loading tidyverse just gets some extra tools from purrr here, but I find it much more productive for data manipulation in general, compared to just using dplyr alone. That said if you can only use dplyr you can replace the purrr::map() for lapply() and purrr::keep() for Filter(), but you lose a little pipe readability.

First time posting on the sub-Exchange, so feedback is welcome. 

df &lt;- read.csv("https://sciences.ucf.edu/biology/d4lab/wp-content/uploads/sites/125/2018/11/parasites.txt", header = T)

m1 &lt;- glm(data=df, infected ~ age + weight + sex, family = "binomial") # add spaces to variables separated by arithmetic operators
link_func &lt;- m1$family$linkinv # maybe this could become a generic function

library(tidyverse)

# anonymous functions are quick and easy to type, my preference if only one input arg
newdat_func &lt;- . %&gt;% # meant to start with df
  select(weight, age) %&gt;% # keep only column of interest
  map(~ round(seq(min(.), max(.), length.out = 15))) %&gt;% # don't repeat yourself and call the same operation on both columns in one line
  c(list(sex = c("female", "male"))) %&gt;% # prep a 3-element list for expand.grid to process
  expand.grid()

newdat2 &lt;- newdat_func(df)

# fall back to traditional function format for multiple inputs
x_func &lt;- function(model, newdata, link_func) {
  predict.glm(model, newdata = newdata, type="link", se=TRUE) %&gt;% # obviously this only works on glm objects, you could add checks to be defensive
    keep(~ length(.) == nrow(newdata)) %&gt;% # drop the third element that is length 1
    bind_cols() %&gt;% # build data frame with a column from each list element
    mutate(low = fit - 1.96 * se.fit,
           high = fit + 1.96 * se.fit) %&gt;%
    mutate_all(funs(link_func)) %&gt;% # again don't repeat yourself
    bind_cols(newdata) %&gt;% # bolt back on simulated predictors
    mutate(category = cut(age,
                          breaks = c(0, 69, 138, 206),
                          labels = c("0-69", "70-139", "139-206")),
           age = as.factor(age))
}

x2 &lt;- x_func(m1, newdat2, link_func)

ggplot(data = x2, aes(x = weight)) + # always use spaces around '+' and '=', do ggplot(data = data) +
  geom_line(aes(y = fit, color = age)) +
  geom_ribbon(aes(ymin = low, ymax = high, fill = age), alpha = 0.1) + # okay is all on one line (&lt;80 chars)
  facet_grid(category ~ sex) +
  labs(x = expression(bold("Weight")), # if a function goes beyond 1 line, split its args one per row
       y = expression(bold(y = "Infection Probability")),
       linetype = "Age (months)",
       colour = "Age (months)",
       fill = "Age (months)") +
  theme(panel.grid.major = element_blank(), # split args again
        panel.grid.minor = element_blank(),
        legend.position = "right",
        strip.text.x = element_text(face = "bold", size=12),
        strip.text.y = element_text(size=10),
        axis.text.y = element_text(size=10, face = "bold"),
        axis.text.x = element_text(size=10),
        axis.title = element_text(size=12), 
        legend.text = element_text(size=10),
        legend.title = element_text(size=12, face="bold"))


Minor tidy-style adjustments everywhere are adding spaces around ~/=/+ signs and only one argument per line for multiy line calls like theme() and labs(). See more here https://style.tidyverse.org/

Obviously I went the last inch and wrapped the processing steps into functions. But I developed the sequence as an open pipe chain, adding a step and printing the result to console as I progressed. The speed of that iterative/dev workflow is why I love leveraging pipes, but I think it also makes the code easier to read. Now instead of multiple intermediate variables and repeated patterns your have two code chunks/functions that handle the two distinct phases of this model plotting problem