4
$\begingroup$

I am currently building a GAM model to describe house prices. The dataset is a collection of roughly 200K house sale prices $\{P_{i}\}$ together with a vector of house characteristics $\{\textbf{x}_{i}\}$. My strategy, which seems to be widely used in the literature on this subject, is to predict the prices using a hedonic regression model of the form:

\begin{align} P = s(x_{1}) + s(x_{2}) + … \end{align}

My approach seems to work rather well as it yields “reasonably looking” splines of the different covariates.

I am however faced with the issue that my database of house sales contains a lot of garbage data points that I suspect might “pollute” the quality of the fitted GAM model. I therefore have the following question: What methods exist for the removal of such data points when fitting GAM-models? I realize that the question is rather broad, so I would be happy to simply hear people’s thoughts on the subject without necessarily needing to go into specific details.

Let me also be clear, that I have intentionally avoided the use of the word “outlier” in my question, as I am not sure whether an outlier filtration mechanism is what I am looking for. Indeed, outlier filtration would simply remove points far from the mass of the data points, but this is not what I am looking for. I would like an approach that identifies data points that exhibit an abnormal influence on the resulting fit.

Having posed the main question, let me elaborate a bit on what my own thoughts have been so far. I can generally come up with two approaches for dealing with junk data.

  1. As a first approach one might consider a strictly data-driven approach in which one looks at the total data set of house prices together with characteristics $\{P_{i}, \textbf{x}_{i}\}$ and removes outliers in this N-dimensional space based on e.g. the Mahalanobis distance measure. Based on my discussion of outliers above, I am however not too fond of this approach.

  2. As an alternative approach I considered taking inspiration from the concept of Cook’s distance for linear models. Here only data points with both high leverage and influence would be removed from the fit, and indeed this idea has already been brought up once in the following post: Can I use the Cook´s Distance to find outliers in a GAM? As I understand from the answers however, it is not clear whether the cook’s distance formula used for linear models is applicable to GAM-models, which has therefore so far led to avoid digging too much into this approach.

Looking forward to hear your thoughts.

$\endgroup$
2
  • $\begingroup$ Do you mean unusual points among the features or values that are not predicted well by the model? Something else? $\endgroup$
    – Dave
    Commented Apr 3 at 14:00
  • $\begingroup$ I mean points that on their own have a high impact on the fitted function. For a linear model those with high cook's distance. $\endgroup$ Commented Apr 3 at 16:05

1 Answer 1

1
$\begingroup$

I suspect you could use Cook's distance, but as I'm not so familiar with these I'll demonstrate a way to do this using {brms}. You can make use of Pareto-smoothed approximate leave-one-out cross-validation (LOOCV) as a useful metric for identifying points with excess leverage (see some details and links to references in this help page from the Stan team). In the below example I simulate a smooth relationship and then add some excess noise to a few of the points to see if this diagnostic can identify them. It works reasonably well and could probably give you useful information about which points might have unnecessarily:

# Simulate a smooth relationship
library(mgcv)
#> Loading required package: nlme
#> This is mgcv 1.8-42. For overview type 'help("mgcv-package")'.
set.seed(555)
x <- rnorm(100)
y <- 0.2 * sin(x) + -0.6 * exp(x) + 0.8 * cos(x) -
  0.15 * x^3 + rnorm(100, sd = 0.3)

# Randomly add major noise to some points
n_garbages <- 6
garbages <- sample(1:length(y), n_garbages, replace = FALSE)
y[garbages] <- y[garbages] + 
  runif(n_garbages, 3.5, 5) * sample(c(-1, 1), n_garbages, 
                                     replace = TRUE)
# Plot the points
plot(x, y, pch = 16, cex = 0.9, 
     col = '#00000066', bty = 'l')

# Fit the GAM in brms (suppressing messages for nicer output)
library(brms)
dat <- data.frame(x, y)
mod <- brm(y ~ s(x),
           data = dat,
           backend = 'cmdstanr', 
           silent = 2, 
           refresh = 0)
#> Warning: 9 of 4000 (0.0%) transitions ended with a divergence.
#> See https://mc-stan.org/misc/warnings for details.

# Visualise the smooth on the outcome scale using 
# marginaleffects capability
library(marginaleffects); library(ggplot2)
plot_predictions(mod, condition = 'x', points = 0.5,
                 rug = TRUE, type = 'prediction') +
  theme_classic()

# Use approximate leave-one-out cross-validation to find
# influential points
(modloo <- loo(mod, k_threshold = 0.5))
#> Warning: Found 4 observations with a pareto_k > 0.5 in model 'mod'. We
#> recommend to set 'moment_match = TRUE' in order to perform moment matching for
#> problematic observations.
#> 
#> Computed from 4000 by 100 log-likelihood matrix
#> 
#>          Estimate   SE
#> elpd_loo   -154.4 19.2
#> p_loo        14.7  5.1
#> looic       308.9 38.4
#> ------
#> Monte Carlo SE of elpd_loo is NA.
#> 
#> Pareto k diagnostic values:
#>                          Count Pct.    Min. n_eff
#> (-Inf, 0.5]   (good)     96    96.0%   373       
#>  (0.5, 0.7]   (ok)        1     1.0%   145       
#>    (0.7, 1]   (bad)       3     3.0%   49        
#>    (1, Inf)   (very bad)  0     0.0%   <NA>      
#> See help('pareto-k-diagnostic') for details.

# Influential points identified; plot each point's 
# pareto_k estimate and label those with pareto_k > 0.5
plot(modloo, label_points = TRUE)

# Which of these were true "outliers"?
sort(garbages)
#> [1]  4 24 52 55 72 81
which(modloo$diagnostics$pareto_k > 0.5) %in% garbages
#> [1]  TRUE FALSE  TRUE  TRUE

# Highlight those points identified as influential in a plot
plot_predictions(mod, condition = 'x', points = 0.5,
                 rug = TRUE, type = 'prediction') +
  geom_point(data = data.frame(x = x, y = y, 
                        outlier = modloo$diagnostics$pareto_k > 0.5),
             aes(x, y, col = outlier)) +
  scale_color_discrete(type = c('grey30', 'firebrick')) +
  theme_classic() +
  theme(legend.position = 'none')

Created on 2024-04-04 with reprex v2.1.0

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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