Why do I get this error when trying to obtain back-transformed estimated marginal means using emmeans? I am trying to examine the interaction between temperature (factor: 6 levels) and species (factor: 2 levels) in determining the duration of the first larval stage in insects using the following model and the package emmeans:
model3<-lm(sqrt(Duration.1st.days)~Temperature*Species, data = first.1)

Analysis of Variance Table

Response: sqrt(Duration.1st.days)
                     Df  Sum Sq Mean Sq  F value    Pr(>F)
Temperature           5 19.2397  3.8479 116.9385 < 2.2e-16
Species               1  2.4146  2.4146  73.3783 < 2.2e-16
Temperature:Species   5  0.7496  0.1499   4.5563 0.0004338
Residuals           625 20.5661  0.0329                   
                       
Temperature         ***
Species             ***
Temperature:Species ***
Residuals              
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

I am using emmeans in order to compare temperature within species and species within each level of temperature however I would like to get the estimated marginal means and their standard errors on the scale of Duration.1st.days rather than sqrt(Duration.1st days). I have tried the following code but get an error when I run the first line of code. When I run the second the EM means are still on the sqrt scale:
Within.species<-emmeans(model3, type = "response", specs = pairwise ~ Temperature|Species)
emmeans:::cld.emmGrid(Within.species$emmeans, Letter = letters, alpha = 0.05)

Note: Use 'contrast(regrid(object), ...)' to obtain contrasts of back-transformed estimates

Would anyone be able to tell me what this error means and how I might rectify it?
 A: If you have estimates of $\sqrt{\mu_1}$ and $\sqrt{\mu_2}$, the difference of those estimates is $\sqrt{\mu_1} - \sqrt{\mu_2}$. There is no meaningful way to back-transform that quantity, so for the pairwise differences, the software leaves it on the difference-of-square-root scale, and gives you a message suggesting something else you could do.
I suggest never using pairwise ~ ... in emmeans() because you are asking for two different sets of results, and they are worth thinking about separately. Here are some suggestions to consider:
EMM <- emmeans(model3, specs = ~ Temperature|Species)

EMM    # show means on square-root scale
pairs(EMM)    # pairwise differences on square-root scale
confint(EMM, type = "response")    # means on response scale

# If you want comparisons on the response scale, the message suggests...
regrid(EMM)    # convert to response scale and create a whole new grid of means
               # Display will be the same as above, but internal structure differs
pairs(regrid(EMM))   # pairwise differences on response scale
cld(regrid(EMM))     # compact letter display

Read the vignette on transformations for more discussion.
For lots of reasons, I don't think much of CLDs; the main reason is that they display what you can't show rather than what you can. I'd suggest this vignette as well.
