```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```



## Species diversity on the Galapagos Islands

```{r}
#install.packages("faraway")
library(faraway)
data(gala)
galaNew <- gala[, -2]
```

## Model Selection

```{r}
library(tidyverse)
library(caret)
library(leaps)
models <- regsubsets(Species ~ ., data = galaNew, nvmax = 5)
summary(models)

res.sum <- summary(models)

criteria <- data.frame(
  Adj.R2 = res.sum$adjr2,
  Cp = res.sum$cp,
  BIC = res.sum$bic)

criteria

plot(2:6, criteria$Cp, las = 1, xlab = "p", ylab = "Cp",
     pch = 16, col = "gray", ylim = c(1, max(criteria$Cp)))
abline(0, 1)

plot(2:6, criteria$Adj.R2, las = 1, xlab = "p", ylab = "", pch = 16, col = "gray",
     main = expression(R['adj']^2))
points(5, criteria$Adj.R2[4], col = "blue", pch = 16)

plot(2:6, criteria$BIC, las = 1, xlab = "p", ylab = "", pch = 16, col = "gray", main = "BIC")
points(3, criteria$BIC[2], col = "blue", pch = 16)


full <- lm(Species ~ ., data = galaNew)
step(full)
step(full, direction = "backward")
```

## Model Diagnostics

```{r}
mod <- lm(Species ~ Elevation + Adjacent, data = galaNew)
plot(mod$fitted.values, galaNew$Species, pch = 16, col = "blue")
abline(0, 1, col = "red")

plot(mod$fitted.values, mod$residuals, pch = 16, col = "blue")
abline(h = 0, col = "red")

par(las = 1)
hist(mod$residuals, 12, prob = T,
     col = "lightblue", border = "gray")
xg <- seq(-200, 200, 1)
yg <- dnorm(xg, 0, 60.86)
lines(xg, yg)


plot(qnorm(1:30 / 31, 0, 60.86), sort(mod$residuals), pch = 16,
     col = "gray", xlab = "Normal Quantiles", ylab = "Residuals")
abline(0, 1)
```

## Diagnostics in Multiple Linear Regression

### Leverage

```{r}
full <- lm(Species ~ ., data = galaNew)
step_gala <- step(full)

X <- model.matrix(step_gala)
H <- X %*% solve((t(X) %*% X)) %*% t(X)
lev <- hat(X)
high_lev <- which(lev >= 2 * 3 / 30)
attach(gala)

par(las = 1)
plot(Elevation, Adjacent,
     cex = sqrt(5 * lev),
     col = "blue", ylim = c(0, 5000))
points(Elevation[high_lev],
       Adjacent[high_lev], col = "red",
       pch = 16,
       cex = sqrt(5 *lev[high_lev]))
```

### Studentized Residuals

```{r}
gs <- summary(step_gala)
gs$sig

studRes <- gs$res / (gs$sig * sqrt(1 - lev))

par(las = 1)
plot(studRes, pch = 16,
     cex = 0.8, col = "blue",
     ylab = expression(r[i]), main = "Studentized Residuals", xlab = "")
abline(h = 0, lty = 2, col = "gray")
```


### Jackknife Residuals

```{r}
jack <- rstudent(step_gala)

par(las = 1)
plot(jack, pch = 16,
     cex = 0.8, col = "blue", main =" Jacknife Residuals ", xlab = "",
     ylab = "")
abline(h = 0, lty = 2, col = "gray")
```




### Identifying Influential Observations: DFFITS
```{r, message=F}
library(olsrr)
ols_plot_dffits(step_gala)
```


### Residual Plot

```{r}
par(las = 1)
plot(step_gala$fitted.values,
     step_gala$residuals, 
     pch = 16, cex = 0.8, col = "blue", main =" Residuals ",
     xlab = expression(hat(Y)), ylab = expression(e))
abline(h = 0, lty = 2, col = "gray")

sqrt_fit <- lm(sqrt(Species) ~ Elevation + Adjacent)

par(las = 1)
plot(sqrt_fit$fitted.values,
     sqrt_fit$residuals, 
     pch = 16, cex = 0.8, col = "blue", main =" Residuals ",
     xlab = expression(hat(Y)), ylab = expression(e))
abline(h = 0, lty = 2, col = "gray")
```



