1
\$\begingroup\$

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<-parasites
m1<-glm(data=df, infected~age+weight+sex, family = "binomial")
summary(m1)
age_grid <- round(seq(min(df$age), max(df$age), length.out = 15))
weight_grid <- round(seq(min(df$weight), max(df$weight), length.out = 15))
newdat <- expand.grid(weight =weight_grid,
                      age = age_grid, sex = c("female", "male")) 

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

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

x$age<-as.factor(x$age)

library(ggplot2)
finalgraph<-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.

\$\endgroup\$
2
\$\begingroup\$

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 <- read.csv("https://sciences.ucf.edu/biology/d4lab/wp-content/uploads/sites/125/2018/11/parasites.txt", header = T)

m1 <- glm(data=df, infected ~ age + weight + sex, family = "binomial") # add spaces to variables separated by arithmetic operators
link_func <- 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 <- . %>% # meant to start with df
  select(weight, age) %>% # keep only column of interest
  map(~ round(seq(min(.), max(.), length.out = 15))) %>% # don't repeat yourself and call the same operation on both columns in one line
  c(list(sex = c("female", "male"))) %>% # prep a 3-element list for expand.grid to process
  expand.grid()

newdat2 <- newdat_func(df)

# fall back to traditional function format for multiple inputs
x_func <- function(model, newdata, link_func) {
  predict.glm(model, newdata = newdata, type="link", se=TRUE) %>% # obviously this only works on glm objects, you could add checks to be defensive
    keep(~ length(.) == nrow(newdata)) %>% # drop the third element that is length 1
    bind_cols() %>% # build data frame with a column from each list element
    mutate(low = fit - 1.96 * se.fit,
           high = fit + 1.96 * se.fit) %>%
    mutate_all(funs(link_func)) %>% # again don't repeat yourself
    bind_cols(newdata) %>% # 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 <- 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 (<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

\$\endgroup\$
  • \$\begingroup\$ Hey, welcome to Code Review! Here we usually look for reviews of the OP's code and not just alternative implementations. While you did write some comments on what is different about your approach, it might help if you could add some more on why/how this is better (is it more readable, more modular or maintainable, faster, more memory efficient, just plain best practices or more modern?). \$\endgroup\$ – Graipher Nov 25 '18 at 16:16
  • \$\begingroup\$ Thanks @Nate! the code worked very well. I do have a question though. For some of your function calls you use "." in the arguments section, e.g. length(.) or keep(~length(.) == nrow(newdat). What exactly does this do? \$\endgroup\$ – Leo Ohyama Nov 26 '18 at 23:16
  • \$\begingroup\$ Yea its from the purrr formula syntax, specified by the ~. See the help/man page for ?map for more, but the . is the reference to the element being passed into the function (each element in the list being mapped over). So the statement ~length(.) is the same as function(x) length(x) its just shorthand. You can always pass the function()... syntax to map(.f = ) and for more complicated functions the ~ shorthand may fail, but its faster for more simple cases. \$\endgroup\$ – Nate Nov 27 '18 at 14:30

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.