autoplot() in R for lm: why do I get a "Constant Leverage: Residuals vs Factor Levels" instead of a "Residuals vs Leverage" plot?

Question

I am doing an ANCOVA in R with one continuous variable (DENSITY) and one factor (SEASON). When I check for the model assumptions I get a plot named: "Constant Leverage: Residuals vs Factor Levels" instead of the "Residuals vs Leverage" plot.

limp.mod <- lm(EGGS~DENSITY*SEASON, data=limp)
autoplot(limp.mod,smooth.colour = NA, which=5)

Image:What I get

Image:What I want

How can I get the "Residuals vs Leverage" plot? Why does the exact same code in my textbook give another autoplot() output?

Thanks in advance for your help!

Answer 1

Without a reproducible example, I will first create a model based on the built-in data set iris.

df1 <- iris[1:100, 3:5]
df1$Species <- droplevels(df1$Species)
str(df1)
#'data.frame':   100 obs. of  3 variables:
# $ Petal.Length: num  1.4 1.4 1.3 1.5 1.4 1.7 1.4 1.5 1.4 1.5 ...
# $ Petal.Width : num  0.2 0.2 0.2 0.2 0.2 0.4 0.3 0.2 0.2 0.1 ...
# $ Species     : Factor w/ 2 levels "setosa","versicolor": 1 1 1 1 1 1 1 1 1 1 ...

fit <- lm(Petal.Length ~ Petal.Width*Species, df1)

As for the plot, autoplot is a generic function. Package ggfortify includes methods for objects of class "lm", among others.

From help("autoplot.lm"):

which If a subset of the plots is required, specify a subset of the numbers 1:6.

The default is which = c(1, 2, 3, 5). Trying all 6 values for the argument, we see that the wanted graph is not one of them. So a custom graph needs to be built.

The residuals and the leverage values can be obtained from stats::resid and from stats::hatvalues, respectively.

library(ggplot2)

dflev <- data.frame(Leverage = hatvalues(fit), y = resid(fit))

ggplot(dflev, aes(Leverage, y)) +
  geom_point() +
  geom_hline(yintercept = 0, linetype = "dashed") +
  ggtitle("Residuals vs Leverage") +
  lims(y = c(-1, 1)) +
  ylab("") +
  theme(plot.title = element_text(hjust = 0.5, face = "bold"))