1
$\begingroup$

I got a residuals vs fitted plot with a clear downward pattern from a linear regression model and am wondering what are some of the sensible things to do in this situation? enter image description here

Doing a boxcox transformation does not help much in correcting the pattern (bc=0.3 from boxcox()). A log transformation (with an offset ) of the responses does not help either. enter image description here

My response variable is vehicle miles traveled, so it is truncated at 0. I am aware that a truncated model such as tobit models likely works better for my data, but the downward pattern persists in the residuals vs fitted plot from the tobit model, so I am hoping to address the issue with the simplest model possible.

In term of independent variables, HHC_MSA is a categorical variable labeling different regions in the US, while HHSIZE is the number of persons in a household. I tried non-linear transformation of the HHSIZE variable with splines::ns, the downward pattern tilts up somewhat, but persists enter image description here.

Any suggestions and comments are appreciated!

UPDATE

Following @Glen_b 's suggestion, I added a residuals vs fitted plot from a hurdle model, however, the downward band persists there, similarly for a zero-inflated negative binomal model and a tobit model.

residuals vs fitted plot from hurdle model <code>mfit3 <- hurdle(as.integer(VMT_MILE)~HHC_MSA+HHSIZE, data=act2)</code>

Below are the code and data that replicate my results (please be aware of large data files - csv files more than 400MB):

library(readr)
library(dplyr)
library(AER)
library(MASS)

#Down data from [ORNL](http://nhts.ornl.gov/download.shtml#2009) and unzip
zip_url <- "http://nhts.ornl.gov/2009/download/Ascii.zip"
data_dir <- "data/"
data_file <- file.path(data_dir, "Ascii/DAYV2PUB.CSV")
zip_file <- file.path(data_dir, "NHTS2009.zip")

if (!dir.exists(data_dir)) dir.create(data_dir)

if (!file.exists(data_file)) {
  download.file(zip_url, destfile=zip_file, method="libcurl")
  unzip(zip_file, exdir = data_dir)
}
stopifnot(file.exists(data_file))

act <- read_csv(file.path(data_dir, "Ascii/DAYV2PUB.CSV"))
act2 <- act %>% 
  mutate(VMT_MILE=ifelse(VMT_MILE < 0, 0, VMT_MILE),
         HHC_MSA=ifelse(HHC_MSA %in% c("XXXX", "-1", "-9"), NA, HHC_MSA)) %>% 
  group_by(HOUSEID) %>% 
  summarize(VMT_MILE=sum(VMT_MILE, na.rm=0),
            HHSIZE=first(HHSIZE),
            HHVEHCNT=first(HHVEHCNT),
            HHC_MSA=first(HHC_MSA)) %>% 
  filter(HHVEHCNT<=12, VMT_MILE<=200)

summary(mfit0 <- lm(VMT_MILE~HHC_MSA+HHSIZE, data=act2))
plot(mfit0, which=1)

summary(mfit0_ns <- lm(VMT_MILE~HHC_MSA+ns(HHSIZE, 3), data=act2))
plot(mfit0_ns, which=1)

# adds 0.1 to VMT_MILE as boxcox requires positive responses
bct <- boxcox(mfit1 <- lm(VMT_MILE+0.1~HHC_MSA+HHSIZE, data=act2), plotit=F)
(bc <- bct$x[which.max(bct$y)])
summary(mfit1 <- lm(((VMT_MILE+0.1)^bc + 1)/bc~HHC_MSA+HHSIZE, data=act2))
plot(mfit1, which=1)

summary(mfit2 <- lm(log(VMT_MILE+0.1)~HHC_MSA+HHSIZE, data=act2))
plot(mfit2, which=1)

library(pscl)
summary(mfit3 <- hurdle(as.integer(VMT_MILE)~HHC_MSA+HHSIZE, data=act2))
residuals_ <- resid(mfit3, type="response")
fitted_ <- fitted(mfit3, type="response")
plot(fitted_, residuals_)
$\endgroup$
  • $\begingroup$ Note that the "top" of the band of residuals is parallel to the line at the bottom. The line at the bottom is cause by the "0" values. Is there a hard maximum on the response ("vehicle miles travelled") in your data set? $\endgroup$ – Glen_b Jul 9 '17 at 5:19
  • $\begingroup$ max(act$VMT_MILE)=3600, but I excluded observations whose VMT_MILE > 200 (99.9 percentile), as lm diagnostic plots identify higher values as potential outliers. $\endgroup$ – LmW. Jul 9 '17 at 5:33
  • 1
    $\begingroup$ That choice (with the 0's) accounts for the down-sloping band in your residuals. Rather than use an obviously inappropriate distributional model, one actually suitable for the kind of measurement you have (where responses are necessarily non-negative but with exact zeros possible) might be better. A zero-inflated gamma GLM / gamma hurdle model would seem to be an obvious thing to consider. An alternative if you really want a linear regression would be to build a logistic regression model for the "0" vs ">0" part and then for the ">0" observations, model log-miles against your predictors... $\endgroup$ – Glen_b Jul 9 '17 at 6:14
  • $\begingroup$ ... That might avoid the need to remove data (with consequent issues that may bring) $\endgroup$ – Glen_b Jul 9 '17 at 6:15
  • $\begingroup$ @Glen_b, thanks and I appreciate your suggestion. I did try hurdle and zero-inflated negative binomial models in the pscl package (after converting my responses to integer), as well as an AER::tobit model, but the downward band persists. I updated my question with a residuals vs fitted plot from the hurdle model. Any additional suggestions are highly appreciated! $\endgroup$ – LmW. Jul 9 '17 at 16:07

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.