Как создать цикл for в r, который назначает ggplots как объект без сохранения поверх них

Вариант map. Прокрутите имена столбцов в виде строки и передайте ее aes_string

library(tidyverse)
v1 <- c('tn', 'tp', 'tss')
out <- map(v1, ~ 
            ggplot()+
     geom_boxplot(
      data = dat, 
      mapping = aes_string(x = "site", y =.x) #or
      # mapping = aes(x, !! rlang::sym(.x))
      ))

Или вы можете использовать .data[[]] для извлечения желаемых имен столбцов

library(tidyverse)

# define plotting function
plot_gg <- function(dat, x_var, y_var) {

  dattn <-  ggplot() +
    geom_boxplot(
      data = dat, 
      mapping = aes(x = .data[[x_var]], y = .data[[y_var]])
    ) +
    labs(x = x_var, y = y_var)
  return(dattn)
}

# save plots in a list
plot_list <- c('tn', 'tp', 'tss') %>% 
  map(~ plot_gg(dat, 'site', .x))
plot_list
#> [[1]]

#> 
#> [[2]]

#> 
#> [[3]]

^{Создано 5 августа 2019 г. пакетом REPEX (v0.3.0)}

Подумайте о том, чтобы преобразовать ваши широкие данные в длинный формат, а затем построить график с facet_wrap или facet_grid без каких-либо сложных циклов, сопоставления или сохранения множества графиков для gridExtra::grid. Ниже демонстрируется случайная выборка данных.

Данные (предполагается, что приведенная ниже структура отражает фактические данные OP)

set.seed(852019)

### DATA BUILD
random_df <- data.frame(
  site = sample(c("sas", "stata", "spss", "python", "r", "julia"), 500, replace=TRUE),
  var2 = NA,
  var3 = NA,
  var4 = NA,
  var5 = NA,
  var6 = NA,
  var7 = NA,
  var8 = NA,
  var9 = NA,
  tn = rnorm(500),
  tss = rnorm(500),
  tp = rnorm(500)
)

head(random_df)

#    site var2 var3 var4 var5 var6 var7 var8 var9         tn         tss          tp
# 1 stata   NA   NA   NA   NA   NA   NA   NA   NA  2.0237416 -1.30919981 -1.71918905
# 2     r   NA   NA   NA   NA   NA   NA   NA   NA  0.6052126 -0.27231149  0.18739618
# 3     r   NA   NA   NA   NA   NA   NA   NA   NA  1.3270657 -0.70308896  0.04996251
# 4   sas   NA   NA   NA   NA   NA   NA   NA   NA -0.8690220  0.09934931 -0.12513907
# 5 julia   NA   NA   NA   NA   NA   NA   NA   NA -1.8871174  0.08761820 -0.45409606
# 6   sas   NA   NA   NA   NA   NA   NA   NA   NA  0.3205017 -0.61696052  0.32586570

Сюжет

# RESHAPE WIDE TO LONG
yvars <- c("tn", "tss", "tp")

long_df <- reshape(random_df, varying = yvars, v.names = "y_value",
                   times = yvars, timevar = "y_var",
                   idvar = c("site"), drop = c(2:9),
                   new.row.names = 1:1E4, direction = "long")

head(long_df)
#    site y_var    y_value
# 1 stata    tn  2.0237416
# 2     r    tn  0.6052126
# 3     r    tn  1.3270657
# 4   sas    tn -0.8690220
# 5 julia    tn -1.8871174
# 6   sas    tn  0.3205017

# BOXPLOT WITH FACET
ggplot() + geom_boxplot(data = long_df, mapping = aes(x = site, y = y_value)) + 
  facet_wrap(~y_var)

Если вам нужны отдельные графики, по-прежнему рассмотрите возможность использования длинного формата и используйте by вне столбца индикатора y_var для построения списка графиков. Затем постройте с gridExtra::grid.arrange:

plot_list <- by(long_df, long_df$y_var, function(sub) {
  ggplot() + geom_boxplot(data = sub, mapping = aes(x = site, y = y_value)) +
    ggtitle(sub$y_var[[1]])
}) 

do.call(grid.arrange, args=list(grobs=plot_list, nrow = 1))

Использование образцов данных OP

# RESHAPE WIDE TO LONG
yvars <- c("tn.1", "tss.1", "tp.1")

long_df <- reshape(df, varying=yvars, v.names="y_value",
                   times = yvars, timevar = "y_var",
                   idvar = c("site"), drop=c(2:9),
                   new.row.names = 1:1E4, direction = "long")
long_df$y_var <- gsub(".1", "", long_df$y_var)

# GRID ARRANGE PLOT
plot_list <- by(long_df, long_df$y_var, function(sub) {
  ggplot() + geom_boxplot(data = sub, mapping = aes(x = site, y = y_value)) +
    ggtitle(sub$y_var[[1]])
}) 

do.call(grid.arrange, args=list(grobs=plot_list, nrow = 1))