I first run the multilevel model for analysis, which has no issue.

mlm4 <- lmer(y ~ grp + t + (1|id), data=d1)
summary(mlm4)

Now, I wrote my own function to obtain some details about the results and diagnosis. However, the last parts with emmeans() do not work. The error message says, “Error in emmeans(p_m, ~p_t | p_g, at = list(p_t = c(-700, -365, 0, 365, : No variable named p_t in the reference grid”

This is my first time writing my own function in R. Please advise!

mlm_sum <- function(p_m, p_t, p_g) {
  sum <- summary(p_m)
  ci <- confint(p_m)
  
  # variance-covariance matric
  mat <- Matrix::bdiag(VarCorr(p_m))
  
  # diagnosis
  pred_res <- data.frame(predicted=predict(p_m), residual=residuals(p_m))
  plot1 <- ggplot(pred_res, aes(x=predicted, y=residual)) + geom_point() + geom_hline(yintercept=0, lty=3)
  plot2 <- ggplot(pred_res, aes(x=residual)) + geom_histogram(bins=20, color="black")
  plot3 <- ggplot(pred_res, aes(sample=residual)) + stat_qq() + stat_qq_line()

  # estimated marginal means
  emm1 <- emmeans(p_m, ~ p_t | p_g, 
          at = list(p_t = c(-700, -365, 0, 365, 700)),
          pbkrtest.limit = 43788)
  
  emm_plot1 <- emmip(p_m, ~ p_t | p_g,
        at = list(p_t = c(-700, -365, 0, 365, 700)),
        pbkrtest.limit = 43788
  )
  
  return(list(sum, ci, mat, plot1, plot2, plot3, emm1, emm_plot1))
}

mlm_sum(mlm4, dif_date, grp_poct)

working function returning model summary, 95% CI, diagnostic plots, and estiamted marginal means and its plot