title Introduction to Statistical Modeling subtitle Outliers author Joris Vankerschaver ```{r include = FALSE} CWD <- read.table("./datasets/01-linear-regression/christ.csv", header = T, sep = ",", dec = ".") model <- lm(CWD.BASA ~ RIP.DENS, data = CWD) model3 <- lm(I(log(CWD.BASA)) ~ RIP.DENS + I(RIP.DENS^2), data = CWD) summary(model) confint(model) needles <- read.table("./datasets/01-linear-regression/needles.txt", header = T, sep = "\t", dec = ".") model_l5 <- lm(length ~ nitrogen * phosphor + potassium + phosphor * residu, data = needles) library(car) body.fat <- read.table("./datasets/01-linear-regression/bodyfatNKNW.txt", header = T, dec = ".") birnh <- read.table("./datasets/01-linear-regression/birnh.txt", header = T, sep = "\t", dec = ".") ``` ## Outliers / Influential observations - Dataset often contains extreme values for outcome $Y$ and/or predictors $X$ - These **can** influence regression line strongly (but don't have to) ## Influence of influential observations ```{r} #| out-width: 4.5in #| fig-width: 6 #| fig-height: 4.5 #| fig-align: center set.seed(12345) n <- 10 X <- sort(3*runif(n)) Y <- X + 0.3*rnorm(n) plot_data_outliers <- function(outlier = NULL) { X_outliers <- c(X, outlier[1]) Y_outliers <- c(Y, outlier[2]) m_regular <- lm(Y ~ X) m_outliers <- lm(Y_outliers ~ X_outliers) plot(X, Y, xlim = c(-1, 4), ylim = c(-2, 5), pch = 19, xlab = "", ylab = "") if (!is.null(outlier)) { # Draw the outlier points(outlier[1], outlier[2], pch = 19, col = "red") # Draw the regression line without the outlier abline(m_regular, lty = "dashed") } abline(m_outliers) } plot_data_outliers() ``` ## Influence of influential observations ```{r} #| out-width: 4.5in #| fig-width: 6 #| fig-height: 4.5 #| fig-align: center plot_data_outliers(c(3, 4)) ``` ## Influence of influential observations ```{r} #| out-width: 4.5in #| fig-width: 6 #| fig-height: 4.5 #| fig-align: center plot_data_outliers(c(0, 4)) ``` ## Influence of influential observations ```{r} #| out-width: 4.5in #| fig-width: 6 #| fig-height: 4.5 #| fig-align: center plot_data_outliers(c(4, -2)) ``` ## Tracking influential observations **Residuals**: - Indicate how far outcome deviates from regression line - Normally distributed with mean 0 and variance $\sigma^2 = MSE$. Hence, can be used to identify extreme outcomes: - 95% of residuals expected in interval $[-2\sigma, 2\sigma]$ - Observations where residual is much larger are probably outliers. ## Extreme outcomes in analysis of larches? ```{r} #| out-width: 4.5in #| fig-width: 6 #| fig-height: 4.5 #| fig-align: center qqnorm(model_l5$resid, pch = 19) ``` ## Tracking influential observations - **Scatterplots** of outcome in function of predictors can be used to identify extreme outcomes and predictors - When multiple predictors, these plots have serious shortcomings ## Multivariate outliers: $Y$ versus $X_1$ or $X_2$ ```{r} #| out-width: 4.5in #| fig-width: 6 #| fig-height: 4.5 #| fig-align: center set.seed(1234) n <- 30 x1 <- 3*runif(n) x2 <- x1 + 0.5*rnorm(n) y <- 2*x1 + x2 + 5*rnorm(n) x1_outlier <- 2.5 x2_outlier <- 0 y_outlier <- 5 range <- function(...) { c(floor(min(...)), ceiling(max(...))) } x1_lim <- range(x1, x1_outlier) x2_lim <- range(x2, x2_outlier) y_lim <- range(y, y_outlier) par(cex = 1.5, mfrow=c(1, 2)) plot(x1, y, xlim = x1_lim, ylim = y_lim) points(x1_outlier, y_outlier, pch = 15, col = "red") plot(x2, y, xlim = x2_lim, ylim = y_lim) points(x2_outlier, y_outlier, pch = 15, col = "red") ``` ## Multivariate outliers: $X_1$ versus $X_2$ ```{r} #| out-width: 4.5in #| fig-width: 6 #| fig-height: 4.5 #| fig-align: center plot(x1, x2, xlim = x1_lim, ylim = x2_lim) points(x1_outlier, x2_outlier, pch = 15, col = "red") ``` ## Leverage - Diagnostic measure to identify influential predictor-observations - Data point has high leverage if it has "extreme" predictor values (low or high) Mathematically: - Weighted distance between predictor for observation $i$ and mean predictor. - How much the $i$th observed value affects the $i$th fitted value: $$ h_{ii} = \frac{\partial \hat{y}_i}{\partial y_i}. $$ ## Interpretation of leverage - If high leverage for $i^{th}$ observation, then - it has predictor values that deviate strongly from the mean - it **possibly** has large influence on regression coefficients and predictions - Leverage is on average $p/n$ with $p$ number of unknown parameters - **Extreme leverage**: larger than $2p/n$ ## Leverage in analysis of larches ```{r} #| out-width: 4.5in #| fig-width: 6 #| fig-height: 4.5 #| fig-align: center p=7 n=26 lev5 <- hatvalues(model_l5) cutoff <- 2*p/n plot(lev5, pch = 20, xlab = "", ylab = "Leverage") abline(h=cutoff, lty = "dashed") text(1.5, 0.91, 1) text(4.5, 0.59, 4) ``` Note: larches model has $p = 7$, so $2p/n = 0.54$. ## Cook's distance - Diagnostic measure for influence of $i^{th}$ observation on all predictions / estimated coefficients. - Cook's distance for $i^{th}$ observation is obtained by comparing each prediction $\hat{Y}_j$ with prediction $\hat{Y}_{j(i)}$ that would be obtained **if $i^{th}$ observation was deleted**: $$ D_i=\frac{\sum_{j=1}^n(\hat{Y}_j-\hat{Y}_{j(i)})^2}{p \cdot\mathrm{MSE}} $$ ## Interpretation Cook's distance - If Cook's distance $D_i$ large, then $i^{th}$ observation has large influence on predictions and coefficients - **Extreme Cook's distance**: exceeds the 50th percentile (median) of the $F_{p,n-p}$-distribution Example: - In analysis of larches is $p=7, n=26$ and the 50\% percentile of $F_{p,n-p}$-distribution 0.94 - Cook's distance of first observation is 1.5 and corresponds to 77\% percentile - Conclusion: first observation has large influence on estimated regression coefficients ## Cook's distance in analysis of larches ```{r} #| out-width: 4.5in #| fig-width: 6 #| fig-height: 4.5 #| fig-align: center cd5 <- cooks.distance(model_l5) plot(cd5, type = "h", ylab = "Cook's distance", xlab = "") abline(h = qf(0.5, 7, 19), lty = "dashed") text(1.5, 1.5, 1) text(6.5, 0.22, 6) text(7.5, 0.24, 7) ``` ## Analysis of larches: residual plots ```{r} par(mfrow=c(2,2)) plot(model_l5) ``` ## DFBETAs On what coefficient(s) will first observation have large influence? - Diagnostic measure for influence of $i^{th}$ observation **on each regression coefficient separately** - DFBETAs for $i^{th}$ observation and $j^{th}$ coefficient is obtained by comparing $j^{th}$ coefficient $\hat{\beta}_j$ with coefficient $\hat{\beta}_{j(i)}$ from model **if $i^{th}$ observation would have been deleted** $$ \textrm{DFBETA}_{j(i)}=\frac{\hat{\beta}_{j}-\hat{\beta}_{j(i)}}{\textrm{SE}(\hat{\beta}_{j})} $$ ## Interpretation DFBETAs - Sign indicates if deleting observation $i$ causes an increase (DFBETA$<0$) or decrease (DFBETA$>0$) in each coefficient - **Extreme DFBETAs**: exceeds 1 in small to moderate datasets, and $2/\sqrt{n}$ in large datasets ## DFBETAs in analysis of larches ```{r} #| out-width: 4.5in #| fig-width: 6 #| fig-height: 4.5 #| fig-align: center dfb5 <- dfbetas(model_l5)[1,] plot(dfb5, pch = 20, ylab = "DFBETA for observation 1", xlab = "Parameter number") abline(h = c(-1, 1), lty = "dashed") ``` ## DFBETAs in analysis of larches First observation has **large influence on interaction between phosphorus and residual ash**: - current coefficient is -598.08 (SE 290.02); - DFBETA is 2.17; - after deletion of first observation, interaction between phosphorus and residual ash will be around $$ -598.08-2.17\times 290.02=-1227.42 $$ ## Histogram and QQ-plot of interaction ```{r} #| out-width: 4.5in #| fig-width: 6 #| fig-height: 4.5 #| fig-align: center inter5 <- needles$phosphor * needles$residu par(mfrow=c(1,2)) hist(inter5, xlab = "phosphor * residual ash", main = "Histogram of phosphor * residu") qqnorm(inter5) ``` ## Analysis of larches after deletion $1^{st}$ observation \small ```{r} needles2 <- needles[-1, ] model_l6 <- lm(length ~ nitrogen * phosphor + potassium + phosphor * residu, data = needles2) ``` ```{r echo=FALSE} summary(model_l6)$coefficients ``` ## Analysis of larches after deletion interaction \small ```{r} model_l7 <- lm(length ~ nitrogen * phosphor + potassium + residu, data = needles) ``` ```{r echo=FALSE} summary(model_l7)$coefficients ``` ## Analysis of larches: final model \small ```{r} model_l8 <- lm(length ~ nitrogen * phosphor + potassium, data = needles) ``` ```{r echo=FALSE} summary(model_l8)$coefficients ``` ## Final analysis: leverage ```{r} #| out-width: 4.5in #| fig-width: 6 #| fig-height: 4.5 #| fig-align: center p2=5 lev8 <- hatvalues(model_l8) cutoff2 <- 2*p2/n par(mfrow=c(1,2)) plot(lev5, pch = 20, xlab = "", ylab = "Leverage") abline(h=cutoff) text(3, 0.91, 1) text(6, 0.59, 4) plot(lev8, pch = 20, xlab = "", ylab = "Leverage", ylim = c(0.0,0.9)) abline(h=cutoff2) text(6, 0.59, 4) ``` ## Final analysis: Cook's distance ```{r} #| out-width: 4.5in #| fig-width: 6 #| fig-height: 4.5 #| fig-align: center cd8 <- cooks.distance(model_l8) plot(cd8, type = "h", ylab = "Cook's distance", xlab = "") text(4.5, 0.12, 4) text(6.5, 0.17, 6) text(7.5, 0.37, 7) ``` ## Final analysis: residual plots \centering ```{r echo=FALSE} par(mfrow=c(2,2)) plot(model_l8) ```