More specifically, I have two predictor (age with three levels and sex with two levels) treated as factors and one response (salary). I want to test if age, sex and their interaction are significant related to salary. After conducting type III two-way anova (unbalanced observations),

age <-as.factor(c(2,3,2,4,2,3,4,3,4,4,3,3,2,2,2,4,4,3,3,3,2,2,2,3,4,4,2))
sex <- as.factor(c(1,0,0,0,0,1,1,1,0,1,0,1,1,0,1,1,0,1,1,0,0,1,0,0,1,0,0))
mod<- lm(salary~age*sex)

Anova(mod, type = "III")

I found that the graph of "residuals vs fitted value" has heteroscedasticity with respect to its variance. So I want to use gls() function in R to make the variance more constant (homoscedasticity). What I did by coding in R was:

mod_gls = gls(salary ~ age*sex)

However, this plot gave me a plot called "standardized residual vs fitted values but still have almost same pattern as the first plot.

So I was wondering if I did gls() wrong and what is the appropriate way?

  • $\begingroup$ Please add a reproducible example for people to work with. $\endgroup$ Feb 1, 2017 at 0:41
  • $\begingroup$ At this point you haven't described the within-group heteroscedasticity structure in your model yet. Check ?gls and have a look at the weights argument. If you need more than that, please provide a workable example as @gung pointed out. $\endgroup$
    – Stefan
    Feb 1, 2017 at 0:46
  • $\begingroup$ @Stefan Thanks for your suggestion. I did take a look at weights argument and found it needs to specify varFunc and etc,. However, due to my limited ability, I cannot figure out how to do that specifically in my case. I just want to make my variance more constant. Also, I have added my original dataset and right now it might be an appropriate reproducible example? $\endgroup$
    – Bratt Swan
    Feb 1, 2017 at 1:49

1 Answer 1


At this point you haven't described the within-group heteroscedasticity structure in your model yet. Try weights=varPower() as shown in the example in ?gls. That gets rid of the heteroscedasticity in your case.


m1 <- gls(salary ~ age*sex)

enter image description here

m2 <- gls(salary ~ age*sex, weights=varPower())

enter image description here

Also if you look in Chapter 5.2.1 (page 208) in Mixed Effects Models in S and S-Plus by Pinheiro and Bates 2000, there is a lot of information on the Variance Functions in nlme. This answer may also be helpful: Regression modelling with unequal variance .


Not the answer you're looking for? Browse other questions tagged or ask your own question.