I built a linear mixed effect model with nlme with Body length, Habitat and Sex as fixed effects. Body length is added for body size correction while the other two both have two levels. In the random part of the model Population was used as a random effect (with 8 levels). To account for the violation of heterogeneity I used the VarIdent function. The model is:
Model<- lme(A~ Body_length+Sex X Habitat, random= ~ 1|Population, method="REML", weights = varIdent(form = ~1|PopulationSex))*
As populations are clearly linked to habitats each population can only belong to one Habitat level. As a result of the connection between Habitat and Population (I think) the DF in case of Habitat is 6 while in case of Sex it is 500. To my knowledge low DF results in less chance to detect significant result and also increases the inaccuracy of the estimation. The model was followed by a contrast analysis. When visually inspecting the emmeans it seems like the Standard Errors estimated are very high and are almost the same for each Habitat X Sex. Also Sex and Habitat seems to have simiar effect on the response variable. I tried to check the data building an identical GLS model (without the random effect) and Habitat also became significant , had the same DF as Sex and the Standard Errors estimates became lower and differed between groups. Therefore the random effect -due to its few levels (as suggested by some authors)- seems to case this discrepancy. At the same time leaving out Population effect from the model would be incorrect as Populations clearly differ for each other and individuals within a population are more alike. So my questions would be: Is the random effect specified correctly? Is there a way to fix my model? If not could you please recommend other type models to be used with this kind of data?