last line of an R function prevent the plotting

Question

I'm really wondering why when I keep the last line (sigma) in my function foo1 my plot call stops working BUT when I delete the last line in foo2 the plot call works fine?!

Note: My understanding is that plot should show anywhere in the function like:

foo3 <- function(x = 1:3){; plot(x); return(x); }; foo3()

I need to keep the last line but I also need to plot, is this fixable?

library(lme4)
library(emmeans)

h <- read.csv('https://raw.githubusercontent.com/hkil/m/master/h.csv')
h$year <- as.factor(h$year)
m <- lmer(scale~ year*group + (1|stid), data = h)

foo1 <- function(fit, plot = T){

  vc <- VarCorr(fit)
  f  <- as.formula(bquote(pairwise ~ .(terms(fit)[[3]])))
  ems <- emmeans(fit, f, infer = c(T, T))

  if(plot) plot(ems)

  ## Why having this line prevents plotting?
  sigma <- sqrt(sum(as.numeric(c(attr(vc[[1]], "stddev"), attr(vc, "sc")))^2)) 
}

###### EXAMPLE: Will NOT plot !!! 
foo1(m, plot = T)

But simply delete the last line of foo1 and plot works fine!

foo2 <- function(fit, plot = T){

  f <- as.formula(bquote(pairwise ~ .(terms(fit)[[3]])))   
  ems <- emmeans(fit, f, infer = c(T, T))

  if(plot) plot(ems)
}

###### EXAMPLE: NOW plots fine !!!
foo2(m, plot = T)          

Answer 1

The underlying cause for the unexpected behavior is that emmeans uses ggplot2 as plotting method. You may examine the code of function emmeans:::.plot.srg. This is why you can store the plot in an object (p) and print it:

foo1 <- function(fit, plot=TRUE) {
  vc <- VarCorr(fit)
  f <- as.formula(bquote(pairwise ~ .(terms(fit)[[3]])))
  ems <- emmeans(fit, f, infer=c(TRUE, TRUE))
  sigma <- sqrt(sum(as.numeric(c(attr(vc[[1]], "stddev"), attr(vc, "sc")))^2))
  p <- plot(ems)
  if(plot) print(p)
  return(sigma)
}

foo1(m, plot=TRUE)  ## plots
# [1] 126.6106

Intuitively we expect the behavior like in your foo3 <- function(x = 1:3){; plot(x); return(x); }; foo3() example:

foo2 <- function(x, plot=TRUE) {
  y <- seq_len(x)^2
  if(plot) plot(seq_len(x), y)
  return(y)
}
foo2(20)  ## plots
# [1]   1   4   9  16  25  36  49  64  81 100 121 144 169
# [14] 196 225 256 289 324 361 400

But this works differently, when one uses ggplot. When we don't store the "ggplot" object, a following line will "overwrite" the plot.

library(ggplot2)
foo3a <- function(x, plot=TRUE) {
  y <- seq_len(x)^2
  if(plot) ggplot(mapping=aes(seq_len(x), y)) + geom_point()
  return(y)  ## try to comment/un-comment this line
}
foo3a(20)  ## won't plot, just output of y (depending on `return`)

Therefore we have to actually print the "ggplot" object (I show this by storing the plot in p in function foo3b), if it's not the last line of the function.

foo3b <- function(x, plot=TRUE) {
  y <- seq_len(x)^2
  p <- ggplot(mapping=aes(seq_len(x), y)) + geom_point()
  if(plot) print(p)
  return(y)  ## try to comment/un-comment this line (works in both cases)
}
foo3b(20)  ## plots
# [1]   1   4   9  16  25  36  49  64  81 100 121 144 169
# [14] 196 225 256 289 324 361 400

Note that the use of return also is useful.

Answer 2

Well, that makes sense because in foo2 plot is the last command in the function and it returns that whereas that is not the case in foo1.

A simple way if you want to plot and then do some additional calculation is to print it in the function.

foo1 <- function(fit, plot = T){

   vc <- VarCorr(fit)
   f <- as.formula(bquote(pairwise ~ .(terms(fit)[[3]])))

   ems <- emmeans(fit, f, infer = c(T, T))
   if(plot) print(plot(ems))
   sigma <- sqrt(sum(as.numeric(c(attr(vc[[1]], "stddev"), attr(vc, "sc")))^2)) 
}

then if you run foo1

x <- foo1(m, plot = T)
x
#[1] 126.61

it returns sigma value from the function and since we are using print it plots it as well.