R not plotting observations with leverage one

Question

I am building a really simple linear model. I want to test if the frass I got over 3 days from butterfly larvae depend upon the food they ate (diet), the butterfly family (the mother line) and subsequent survival (called "survived", obviously larvae which may latter die are likely to show e.g. problems to eat at larval stage). I'm also interested in the two way interactions: diet:family and survived:family. The interaction diet:survived could not be included because there is only one individual in one of the two diet that died.

Model:

mod=lm(log(frass.weight)~diet*family+survived+family:survived,data=dat)
Anova(mod)   # all the variables are significant.
summary(mod)  #the R2adj is of 0.81
shapiro.test(resid(mod)) # p-value = 0.2389

I have not looked at the variance as I have a small sample size. Only 3 individuals of seven family have been recorded for their frass on both the 2 diets.

Problem:

All looks nice, except that when I plot the model plot(mod) I get the following warning:

"Warning messages: 1: not plotting observations with leverage one:  
 3, 20, 30, 35 "

Is it just a warning or I have a real issue that these points clearly influence the variance?

When I remove these points, the final model I get is simplier:

mod1=lm(log(frass.weight)~diet*family+survived,data=datout)

The residuals are good and the plot works now fine.

Therefore, is the warning about the leverage something to not really consider and my first model should be kept or not? Are my points real outliers?

Answer 1

Points with leverage equal 1.0 indicates something strange about your model, it could be groups with only one member. But it could also be some error in model formulation, so should in all cases be thoroughly checked out! The leverages are the diagonal elements of the hat matrix $H= X (X^T X)^{-1} X^T$ and those diagonals satisfy $0 \le h_i \le 1$. The fitted values can be written $\hat{Y} = HY$ explaining the name "hat matrix": it puts the hat on $Y$ (it was named so by W Tukey). Writing this out for observation number $i$ we get: $$ \hat{y}_i = \sum_{j=1}^n h_{ij} y_j = h_i y_i + \sum_{j\not = i} h_{ij} y_j. $$ When $h_i=1.0$ we can show (*) that the other $h_{ij}=0, j\not= i$. So when $h_i$=1, we can se: $$ \hat{y_i} = y_i. $$ That is, only its own observation influences its prediction. All the other observations are irrelevant for observation $i$! That is quite strange, and could well be caused by some error in model formulation. You should investigate your model (and data!) and try to understand why that happens.

Finally, let us prove our claim (*): $H$ is an idempotent matrix, that is, $HH=H$. Using that we get for the element $h_i$: $$ h_i =\sum_{j=1}^n h_{ij}^2 = h_i^2 + \sum_{j\not = i} h_{ij}^2 $$ which simplifies to $1=1+ \sum_{j\not = i} h_{ij}^2$ proving the claim.

Answer 2

This can happen rather easily. Just fit a model that fits lines for each of two or more groups. In any group where there are just two points, both points will have leverage 1.0:

dat <- data.frame(y=c(1,1.75,3.5, 4.25, 1.5, 4), x=c(1:4,1.5, 3.5),
     gps=rep(c("one" ,"two"),c(4,2)))
dat.lm <- lm(y ~ 0+gps/x, data=dat)
hatvalues(dat.lm)

Output is:

  1   2   3   4   5   6 
0.7 0.3 0.3 0.7 1.0 1.0 

Either of points 5 or 6 totally determines its own predicted value, which is the same as its observed value. It is not so very strange!

Incidentally, the numbers that are given when noting points with leverage 1.0 are row numbers, not the row labels that by default are used to identify the most extreme points in diagnostic plots. (If a data frame where the row labels are row numbers is subsetted, the existing row labels are transferred across, and in general are no longer the row numbers in the subsetted dataset.)