Table of Contents

• Load R packages
• The dataset
• Expanding binomial to binary with pmap_dfr()
• Aside: pmap functions with more columns
• Comparing analysis results
  • Binomial model
  • Binary model
• Expanding binomial to binary via nesting
• Just the code, please

Load R packages

I’m using purrr for looping through rows with pmap_dfr(). I also load dplyr and tidyr for a nest()/unnest() approach.

library(purrr) # 0.3.2
library(tidyr) # 1.0.0
library(dplyr) # 0.8.3

The dataset

I created a dataset with a total of 8 plots, 4 plots in each of two groups. The data has been summarized up to the plot level. The number of trials (total) per plot varied. The number of successes observed is in num_dead.

dat = structure(list(plot = structure(1:8, .Label = c("plot1", "plot2", 
"plot3", "plot4", "plot5", "plot6", "plot7", "plot8"), class = "factor"), 
    group = structure(c(1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L), .Label = c("g1", 
    "g2"), class = "factor"), num_dead = c(4L, 6L, 6L, 5L, 1L, 4L, 
    3L, 2L), total = c(5L, 7L, 9L, 7L, 8L, 10L, 10L, 7L)), class = "data.frame", row.names = c(NA, 
-8L))

dat
#    plot group num_dead total
# 1 plot1    g1        4     5
# 2 plot2    g1        6     7
# 3 plot3    g1        6     9
# 4 plot4    g1        5     7
# 5 plot5    g2        1     8
# 6 plot6    g2        4    10
# 7 plot7    g2        3    10
# 8 plot8    g2        2     7

Expanding binomial to binary with pmap_dfr()

To make the binomial data into binary data, I need to make a vector with a \(1\) for every “success” listed in num_dead and a \(0\) for every “failure” (the total number of trials minus the number of successes). Since I want to do a rowwise operation I’ll use one of the pmap functions. I want my output to be a data.frame so I use pmap_dfr().

I use an anonymous function within pmap_dfr() for creating the output I want from each row. I purposely make the names of the function arguments match the column names. You can either match on position or on names in pmap functions, and I tend to go for name matching. You can use formula coding with the tilde in pmap variants, but I find the code more difficult to understand when I have more than three or so columns.

Within the function I make a column for the response variable, repeating \(1\) num_dead times and \(0\) total - num_dead times for each row of the original data. I’m taking advantage of recycling in data.frame() to keep the plot and group columns in the output, as well.

binary_dat = pmap_dfr(dat, 
                      function(group, plot, num_dead, total) {
                           data.frame(plot = plot,
                                      group = group,
                                      dead = c( rep(1, num_dead),
                                                rep(0, total - num_dead) ) )
                      }
)

Here are the first 6 rows of this new dataset. You can see for the first plot, plot1, there are five rows of \(1\) and one row of \(0\).

head(binary_dat)
#    plot group dead
# 1 plot1    g1    1
# 2 plot1    g1    1
# 3 plot1    g1    1
# 4 plot1    g1    1
# 5 plot1    g1    0
# 6 plot2    g1    1

This matches the information in the first row of the original dataset.

dat[1, ]
#    plot group num_dead total
# 1 plot1    g1        4     5

Aside: pmap functions with more columns

My anonymous function in pmap_dfr() works fine in its current form as long as every column is included as a function argument. If I had extra columns that I didn’t want to remove and wasn’t using in the function, however, I would get an error.

To bypass this problem you can add dots, ..., to the anonymous function to refer to all other columns not being used.

function(group, plot, num_dead, total, ...)

Comparing analysis results

While I definitely learned that binomial data can be analyzed in binary format and returns identical results in a GLM class, for some reason I often have to re-convince myself this is true. ???? This is clear when I do an analysis with each dataset and compare results.

fit = glm( cbind(num_dead, total - num_dead) ~ group, 
           data = dat,
           family = binomial)
summary(fit)$coefficients
#              Estimate Std. Error   z value     Pr(>|z|)
# (Intercept)  1.098612  0.4364358  2.517237 0.0118279240
# groupg2     -2.014903  0.5748706 -3.504968 0.0004566621

fit_binary = glm( dead ~ group, 
                  data = binary_dat,
                  family = binomial)
summary(fit_binary)$coefficients
#              Estimate Std. Error   z value     Pr(>|z|)
# (Intercept)  1.098612  0.4364354  2.517239 0.0118278514
# groupg2     -2.014903  0.5748701 -3.504971 0.0004566575

Expanding binomial to binary via nesting

Doing the expansion with nesting plus purrr::map() inside mutate() is another option, although this seems less straightforward to me for this particular case. It does keep the other variables in the dataset, though, without having to manually include them in the output data.frame like I did above.

binary_dat2 = dat %>%
     nest(data = c(num_dead, total) ) %>%
     mutate(dead = map(data, ~c( rep(1, .x$num_dead),
                                 rep(0, .x$total - .x$num_dead) ) ) ) %>%
     select(-data) %>%
     unnest(dead)
head(binary_dat2)
# # A tibble: 6 x 3
#   plot  group  dead
#   <fct> <fct> <dbl>
# 1 plot1 g1        1
# 2 plot1 g1        1
# 3 plot1 g1        1
# 4 plot1 g1        1
# 5 plot1 g1        0
# 6 plot2 g1        1

Just the code, please

Here’s the code without all the discussion. Copy and paste the code below or you can download an R script of uncommented code from here.

library(purrr) # 0.3.2
library(tidyr) # 1.0.0
library(dplyr) # 0.8.3

dat = structure(list(plot = structure(1:8, .Label = c("plot1", "plot2", 
"plot3", "plot4", "plot5", "plot6", "plot7", "plot8"), class = "factor"), 
    group = structure(c(1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L), .Label = c("g1", 
    "g2"), class = "factor"), num_dead = c(4L, 6L, 6L, 5L, 1L, 4L, 
    3L, 2L), total = c(5L, 7L, 9L, 7L, 8L, 10L, 10L, 7L)), class = "data.frame", row.names = c(NA, 
-8L))

dat

binary_dat = pmap_dfr(dat, 
                      function(group, plot, num_dead, total) {
                           data.frame(plot = plot,
                                      group = group,
                                      dead = c( rep(1, num_dead),
                                                rep(0, total - num_dead) ) )
                      }
)

head(binary_dat)
dat[1, ]

function(group, plot, num_dead, total, ...)
     
fit = glm( cbind(num_dead, total - num_dead) ~ group, 
           data = dat,
           family = binomial)
summary(fit)$coefficients

fit_binary = glm( dead ~ group, 
                  data = binary_dat,
                  family = binomial)
summary(fit_binary)$coefficients

binary_dat2 = dat %>%
     nest(data = c(num_dead, total) ) %>%
     mutate(dead = map(data, ~c( rep(1, .x$num_dead),
                                 rep(0, .x$total - .x$num_dead) ) ) ) %>%
     select(-data) %>%
     unnest(dead)
head(binary_dat2)