diff --git a/R/fitMod_helpers.R b/R/fitMod_helpers.R
index 33891cb..fb72d56 100644
--- a/R/fitMod_helpers.R
+++ b/R/fitMod_helpers.R
@@ -453,7 +453,7 @@ plotFunc <- function(x, CI = FALSE, level = 0.95,
     if(plotData %in% c("means", "meansCI")){
       ## produce estimates for ANOVA type model
       data$doseFac <- as.factor(data[[doseNam]])
-      form <- as.formula(paste(respNam, "~ doseFac +", addCovars[2]))
+      form <- as.formula(paste0("`", respNam, "` ~ doseFac + ", deparse(addCovars[[2]])))
       fit <- lm(form, data=data)
       ## build design matrix for prediction
       dose <- sort(unique(data[[doseNam]]))
@@ -548,7 +548,7 @@ plotFunc <- function(x, CI = FALSE, level = 0.95,
   ylim <- c(rng[1] - 0.05 * dff, rng[2] + 0.05 * dff)
   
   ## produce plot
-  form <- as.formula(paste(respNam, "~", doseNam, "|model", sep=""))
+  form <- as.formula(paste0("`", respNam, "` ~ `", doseNam, "` | model"))
   print(
     lattice::xyplot(form, groups = plotdf$group, data = plotdf, pList=pList, ...,
            ylim = ylim, panel = function(x, y, ..., pList){
diff --git a/tests/testthat/test-fitMod.R b/tests/testthat/test-fitMod.R
index 606f29c..0c54853 100644
--- a/tests/testthat/test-fitMod.R
+++ b/tests/testthat/test-fitMod.R
@@ -286,3 +286,16 @@ test_that("fitMod complains if `resp` is a row-vector", {
                         type = "general", bnds = defBnds(max(doses))$emax),
                  "resp_row is not a numeric but a matrix, converting with as.numeric()")
 })
+
+test_that("plotting works when dose/resp are supplied as expressions", {
+  data(IBScovars)
+  anovaMod <- lm(resp~factor(dose)+gender, data=IBScovars)
+  drFit <- coef(anovaMod)[2:5]
+  vCov <- vcov(anovaMod)[2:5,2:5]
+  gfit <- fitMod(sort(unique(IBScovars$dose))[-1], drFit, S=vCov, model = "emax",
+                 type = "general", placAdj = TRUE, bnds = c(0.01, 2))
+  pdf(NULL)
+  on.exit(dev.off())
+  expect_no_error(plot(gfit, CI=TRUE, plotData = "meansCI"))
+  expect_no_error(plot(gfit, CI=TRUE, plotData = "means"))
+})