I'm attempting to calculate the standard error (SE) of several experimental groups for the purpose of making a plot. However, the data do not satisfy homogeneity of variances- the difference in variance of fitness between treatments is quite large and cannot be solved by the transformations I've tried.
My model is pretty simple: Fitness ~ History * Treatment.
In R, I've tried using emmeans with my model as an lm, and that gives the exact same SE for each group as expected, since it assumes homogeneity of variances. I've read that the gls function of the package nlme should solve this issue here by allowing heterogeneity of variances, but even running emmeans on nlme::gls() gives extremely similar SE among groups (below).
> SM2 <- gls(Seed_mass_mg~History*Treatment, data1, na.action = na.omit)
emmeans(SM2, ~History * Treatment)
History Treatment emmean SE df lower.CL upper.CL
ancestral drought 35.0 5.93 909 23.3 46.6
dry drought 56.3 6.29 909 44.0 68.7
wet drought 39.1 6.12 909 27.1 51.1
ancestral watered 102.9 6.02 909 91.1 114.8
dry watered 131.0 6.38 909 118.5 143.6
wet watered 140.2 5.97 909 128.4 151.9
Degrees-of-freedom method: df.error
Confidence level used: 0.95
However, when I calculate the SE by formula I get quite different SE:
History Treatment Seed_mass_mg_SE
Ancestral Watered 7.008392
Ancestral Drought 1.60024
Drought Watered 8.693766
Drought Drought 2.740732
Watered Watered 9.229806
Watered Drought 2.234901
Can anyone help me understand what I'm misunderstanding about SE here?
A gls() call without a weights argument is just like a lm() call, because it defaults to a homoscedastic model. In particular, I suggest adding weights = varIdent(form = ~ 1 | History*Treatment). See the documentation for varIdent, and take another look at that FAQs vignette in emmeans.