myDF = data.frame(cbind(               # Create some data w/missing variables
  c(1, 3, 4, NA, 2, 5, NA, 7, 9, 3, 2, 5),
  c(7, 8, 2, 4, 3, 6, 9, 1, 0, 4, 2, NA),
  c(6, NA, 3, NA, 4, 6, 1, 9, NA, 2, NA, 8),
  c(5, 3, 2, NA, 2, 1, 8, 6, 7, 4, 7, 9)))

rownames(myDF) = paste(                # Paste together observation labels
  rep("Subject", 12),                  #  by concatenating "subject"
  c(seq(1:12)))                        #  with numbers 1-12

# Missing Values Grid ----------------------------------------------------------

missobs = rowSums(is.na(myDF))         # Num of missing values by obs
missvars = colSums(is.na(myDF))        # Num of missing values by col
myDF2 = myDF[order(missobs),           # Reorder data by num missing obs
  order(missvars)]                     #   & number of missing by variable

png(filename = "illustrations/fig-11-2-missing-values-grid.png",
  units = "in",                        # Set measurements in inches
  res = 1200,                          # Set resolution at 1200dpi
  width = 6,                           # Width at 6 inches
  height = 6)                          # Height at 6 inches

par(mai = c(.25, 1.25, .5, .25))       # Set margins for special labels

plot(x = 0, y = 0,                     # Set up a couple of dummy points
  type = "n",                          # Turn off display of points
  xlim = c(.5, ncol(myDF2) + .5),      # Set range for x axis
  ylim = c(.5, nrow(myDF2) + .2),      # Set range for y axis
  xaxt = "n", xlab = "",               # Turn off x axis values & labels
  yaxt = "n", ylab = "",               # Turn off y axis values & labels
  bty = "n")                           # Turn off surrounding box

axis(3,                                # Put x axis at top
  at = seq(1:ncol(myDF2)),             # Set location of labels
  labels = colnames(myDF2),            # Set labels
  las = 2,                             # Rotate text perpendicular
  lty = 0)                             # Turn off tick marks & axis line

axis(2,                                # Setup y axis
  at = seq(1:nrow(myDF2)),             # Set location of labels
  labels = rownames(myDF2),            # Set labels
  las = 1,                             # Rotate labels perpendicular to axis
  lty = 0)                             # Turn off tick marks & axis line

for(i in seq(1, nrow(myDF2))){         # For each row number
  for(j in seq(1, ncol(myDF2))){       #   and each column number
    polygon(                           # Define a box - x values
      x = c(j-.5, j+.5, j+.5, j-.5),  
      y = c(i-.5, i-.5, i+.5, i+.5),   # y values for box
      col = ifelse(is.na(myDF2[i,j]),  # Set box color depending on NA
        "white", "gray"))              # White for NA, otherwise gray
  }                                    # End for j  loop
}                                      # End for i loop
    
dev.off()                              # Output png file