Using dplyr and pipes for logistic regression plotting I have successfully created a plot of a binomial glm using example data. https://sciences.ucf.edu/biology/d4lab/wp-content/uploads/sites/125/2018/11/parasites.txt

The predictors of the model include 3 predictors (one categorical, 2 continuous)

The code works fine but I have been wanting to try and incorporating more dplyr functions and pipes to streamline code. Ultimately, I want to make my block of code into a function that works with any model with the same type and number of predictors for a binomial glm. Are there better ways to carry out my code with more tidyverse/dplyr code?

#import parasites file

df&lt;-parasites
m1&lt;-glm(data=df, infected~age+weight+sex, family = "binomial")
summary(m1)
age_grid &lt;- round(seq(min(df$age), max(df$age), length.out = 15))
weight_grid &lt;- round(seq(min(df$weight), max(df$weight), length.out = 15))
newdat &lt;- expand.grid(weight =weight_grid,
                      age = age_grid, sex = c("female", "male")) 

pred &lt;- predict.glm(m1, newdata = newdat, type="link", se=TRUE)
ymin &lt;- m1$family$linkinv(pred$fit - 1.96 * pred$se.fit)
ymax &lt;- m1$family$linkinv(pred$fit + 1.96 * pred$se.fit)
fit &lt;- m1$family$linkinv(pred$fit) 
z &lt;- matrix(fit, length(age_grid))
ci.low &lt;- matrix(ymin, length(age_grid))
ci.up &lt;- matrix(ymax, length(age_grid))

x&lt;-data.frame(pred = fit,
              low = ymin,
              high = ymax,
              newdat) %&gt;% mutate(category=cut(age, breaks=c(0, 69, 138, 206), labels = 
                                                c("0-69", "70-139", "139-206")))

x$age&lt;-as.factor(x$age)

library(ggplot2)
finalgraph&lt;-ggplot(data=x)+
  geom_line(aes(x = weight, y = pred, color = age))+
  geom_ribbon(aes(x = weight, ymin = low, ymax = high, fill = age), alpha = 0.1) +
  facet_grid(category~sex) +theme(panel.grid.major = element_blank(),
                                  panel.grid.minor = element_blank())+
  ylab(expression(bold(y = "Infection Probability"))) + xlab(expression(bold("Weight"))) +
  theme(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"))+ 
  labs(linetype="Age (months)", colour="Age (months)", fill = "Age (months)")
finalgraph


Code notes:
Essentially I made a model, created a bunch of values from my predictors (age_grid, v_grid) and made all possible combinations of these values along with the categorical variable of sex using expand.grid.

Then I just used the predict.glm function to extract predicted values based off of expand.grid object. I also extracted std. errors and calculated confidence intervals (ci.up and ci. low). Then I used some dplyr functions to create a dataframe with all this information and also made a new column called category. Category breaks down one of my variables (age) into four distinct groups based of  f of breaks I decided on and labelled as decided as well.  Then I plotted all of this data using ggplot2.
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