En este punto, todavía no ha descrito la estructura de heterocedasticidad dentro del grupo en su modelo. Pruebe weights=varPower() como se muestra en el ejemplo de ?gls . Eso elimina la heterocedasticidad en su caso.
Compare:
m1 <- gls(salary ~ age*sex) plot(m1)
m2 <- gls(salary ~ age*sex, weights=varPower()) plot(m2)
También si busca en el Capítulo 5.2.1 (página 208) en Modelos de efectos mixtos en S y S -Además de Pinheiro y Bates 2000, hay mucha información sobre las funciones de variación en nlme. Esta respuesta también puede ser útil: Modelado de regresión con varianza desigual .