Defining contrast between different levels of treatments in emmeans

Unfortunately, I do not know how to achieve this with emmeans, but it is easy to do using the marginaleffects package. (Disclaimer: I am the author.)

First, we load the package and fit a model. Since you did not provide data, I am using the mtcars dataset, where gear is your treatment and hp is your fertilizer:

library(marginaleffects)

mod <- lm(mpg ~ factor(gear) + hp, mtcars)

Then, we create a prediction grid, with one row per treatment arm (gear), and different baseline values of fertilizer (hp) for each treatment:

dat <- data.frame(
    gear = c(3, 4, 5),
    hp = c(100, 110, 110))

dat
#>   gear  hp
#> 1    3 100
#> 2    4 110
#> 3    5 110

Now, we use the predictions() function to compute adjusted predictions for each of the groups:

predictions(mod, newdata = dat)
#>   rowid     type predicted std.error statistic       p.value
#> 1     1 response  21.19650 1.1705800  18.10769  2.771147e-73
#> 2     2 response  23.16282 0.9381255  24.69054 1.351597e-134
#> 3     3 response  27.10272 1.6979850  15.96170  2.362147e-57
#>   conf.low conf.high gear  hp
#> 1 18.79868  23.59433    3 100
#> 2 21.24116  25.08448    4 110
#> 3 23.62456  30.58089    5 110

Finally, we use the hypothesis argument to specify a custom contrast. This will use the delta method to compare the predictions in the first treatment level to the average of predictions in the other two treatment arms:

predictions(mod, newdata = dat,
    hypothesis = c(-1, .5, .5))
#>       type   term predicted std.error statistic     p.value conf.low
#> 1 response custom  3.936269   1.26463  3.112586 0.001854562  1.45764
#>   conf.high
#> 1  6.414898

The hypothesis argument is quite powerful. You can read more about it here: https://vincentarelbundock.github.io/marginaleffects/articles/hypothesis.html

As I read it, the stated objectives are to obtain predictions (not marginal means) for each combination of treatment and dose, and to compare each of those with the prediction for the control condition. If that is so, then two things stand out in the procedure displayed that are contrary to those goals:

1. Obtaining marginal. This averages over the doses (as is implied BTW by the term "marginal"), thus not providing estimates for the combinations of the two factors.
2. Using cld(). This is a way of doing all pairwise comparisons, as opposed to doing selected ones, i.e. each with control. (In addition, CLDs are at best a poor way of showing pairwise comparisons.)
3. I know I said two things, but I also question the model you fitted; and the meaning of dose2 is not explained in the OP. I suspect it creates a nesting structure that accounts for the fact that the control group can have no dose effect.

Here is something that might work. First, obtain the means for each condition:

M <- emmeans(ModelFert, ~ treatment * dose * dose2)
M   # displays the estimates

One of the estimates shown (I hope) is that of the control group that we want to contrast everything with. Suppose it is the 7th one. Then obtain the contrasts accordingly, changing 7 with the appropriate index:

contrast(M, "trt.vs.ctrl", ref = 7)

All this said, since you are interested in combinations of the factors, I think it is probably a mistake to have a model with only additive effects of treatment and dose. So I strongly recommend using the model with interaction:

ModelFert2 <- lmer (yield ~ treatment * dose + dose2 + 
          (1 | environment/block/cultivar) + 
          (0 + rate_mgha + rate2_mgha || treatment:environment),
          data = fert)

Then use analogous code M2 <- emmeans(ModelFert2, ...) and contrast(M2, ...).

Addendum

Now that I understand that dose is a quantitative predictor, and that dose2 equals dose^2, the plot thickens somewhat.

The second matter first: emmeans assumes that each predictor can be varied independently of the others. That is not the case here because you cannot have dose2 at any value other than dose^2. There are various remedies, but the easiest and most reliable one is to re-fit the model:

ModelFert3 <- lmer (yield ~ treatment * (dose + I(dose^2)) + 
              (1 | environment/block/cultivar),
              data = fert)

Yes, I put the interaction in the fixed-effects part of the model, where it belongs. I would never fit the model shown in the OP; I should have examined it more closely.

That done, the main issue is specifying the different dose levels and then getting the contrasts. Here, I suggest

EMM <- emmeans(ModelFert3, ~ treatment * dose,
               at = list(dose = c(0, 40, 80, 120, 160))
EMM <- EMM[-c(2:5, 6, 11, 16, 21)]
EMM    # List the means to make sure you have it right
contrast(EMM, "trt.vs.ctrl1")

The first statement creates a grid of 25 means for the combinations of 5 treatments and 5 doses. However, you don't want some of them because the first treatment occurs only with dose = 0, and dose = 0 occurs only with the first treatment. This leaves us with the 17 combinations we are actually interested in. Then the contrast call compares the first one with each of the others.

Finally, I think we got it. The code that work is:

emm1 = emmeans(ModelFert,~treatment*dose, at=list (dose= 100))
emm2 = emmeans(ModelFert,~treatment*dose, at=list (dose= 0))
emm3 = emm2[1] + emm1      # or rbind(emm2[1], emm1)
     
confint(emm3, adjutst= "none")

contrast(emm3, "trt.vs.ctrl")

The important thing in these cases is to create two sets of emmeans, and then merge them to finally do the contrast between treatments at different levels. I found the clue in a post made from Russ Lenth!