Homework 8

Part 1: Swiss Fertility Data

swiss_data <- read.table("swiss.txt", header = TRUE)  # load data, header row present
str(swiss_data)  # check structure of data

## 'data.frame':    47 obs. of  6 variables:
##  $ Fertility       : num  64.1 75.5 54.3 79.4 76.9 35 72.7 65 67.6 83.8 ...
##  $ Agriculture     : num  62 85.9 15.2 64.9 43.5 1.2 16.7 75.9 18.7 70.2 ...
##  $ Examination     : int  21 3 31 7 17 37 22 9 25 16 ...
##  $ Education       : int  12 2 20 3 15 53 13 9 7 7 ...
##  $ Catholic        : num  8.52 99.71 2.15 98.22 5.16 ...
##  $ Infant.Mortality: num  16.5 15.1 10.8 20.2 20.6 18 18.9 17.8 19.5 23.6 ...

# fit model predicting Fertility from the other 5 variables
swiss_lm <- lm(Fertility ~ Agriculture + Examination + Education + Catholic + Infant.Mortality,
               data = swiss_data)
summary(swiss_lm)

## 
## Call:
## lm(formula = Fertility ~ Agriculture + Examination + Education + 
##     Catholic + Infant.Mortality, data = swiss_data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -15.2743  -5.2617   0.5032   4.1198  15.3213 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      66.91518   10.70604   6.250 1.91e-07 ***
## Agriculture      -0.17211    0.07030  -2.448  0.01873 *  
## Examination      -0.25801    0.25388  -1.016  0.31546    
## Education        -0.87094    0.18303  -4.758 2.43e-05 ***
## Catholic          0.10412    0.03526   2.953  0.00519 ** 
## Infant.Mortality  1.07705    0.38172   2.822  0.00734 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.165 on 41 degrees of freedom
## Multiple R-squared:  0.7067, Adjusted R-squared:  0.671 
## F-statistic: 19.76 on 5 and 41 DF,  p-value: 5.594e-10

Leverage

n <- nrow(swiss_data)       # number of observations
p <- length(coef(swiss_lm)) # number of parameters (including intercept)
h <- hatvalues(swiss_lm)    # leverage values (diagonal of hat matrix)

plot(h, type = "h", xlab = "Observation Number", ylab = "Leverage",
     main = "Leverage vs Observation Number")
points(h, pch = 16)
abline(h = 2 * p / n, col = "blue", lty = 2)   # 2p/n threshold
abline(h = 3 * p / n, col = "red", lty = 2)    # 3p/n threshold
legend("topleft", legend = c("2p/n", "3p/n"), col = c("blue", "red"), lty = 2)

which(h > 2 * p / n)   # observations exceeding 2p/n

##    La Vallee V. De Geneve 
##            3            6

which(h > 3 * p / n)   # observations exceeding 3p/n

## V. De Geneve 
##            6

2p/n = 0.2553 and 3p/n = 0.383. Two provinces exceed 2p/n: La Vallee (h = 0.351) and V. De Geneve (h = 0.456). Only V. De Geneve exceeds 3p/n.

Studentized Deleted Residuals

stud_del_res <- rstudent(swiss_lm)  # studentized deleted (externally studentized) residuals

plot(stud_del_res, xlab = "Observation Number", ylab = "Studentized Deleted Residual",
     main = "Studentized Deleted Residuals vs Observation Number")
abline(h = c(-2, 2), col = "red", lty = 2)

order(abs(stud_del_res), decreasing = TRUE)[1:3]           # indices of 3 most extreme

## [1] 32 25 39

rownames(swiss_data)[order(abs(stud_del_res), decreasing = TRUE)[1:3]]  # corresponding provinces

## [1] "Sierre"      "Rive Gauche" "Porrentruy"

stud_del_res[order(abs(stud_del_res), decreasing = TRUE)[1:3]]

##      Sierre Rive Gauche  Porrentruy 
##    2.445227   -2.394471   -2.367218

The two high-leverage points do not have large studentized residuals: La Vallee = 0.61, V. De Geneve = 0.04, both well within ±2. The three most extreme provinces by studentized deleted residual are Sierre (2.45), Rive Gauche (-2.39), and Porrentruy (-2.37).

Cook’s Distance

cooksd <- cooks.distance(swiss_lm)  # Cook's distance for each observation

plot(cooksd, type = "h", xlab = "Observation Number", ylab = "Cook's Distance",
     main = "Cook's Distance vs Observation Number")
points(cooksd, pch = 16)
abline(h = 4 / n, col = "red", lty = 2)   # common rule-of-thumb threshold

order(cooksd, decreasing = TRUE)[1:5]

## [1] 39 32 47 25 28

rownames(swiss_data)[order(cooksd, decreasing = TRUE)[1:5]]

## [1] "Porrentruy"  "Sierre"      "Neuchatel"   "Rive Gauche" "Rive Droite"

cooksd[order(cooksd, decreasing = TRUE)[1:5]]

##  Porrentruy      Sierre   Neuchatel Rive Gauche Rive Droite 
##   0.2077343   0.1476237   0.1249500   0.1123648   0.1024707

Cook’s distance for the two high-leverage points is small: La Vallee = 0.034, V. De Geneve = 0.0002. The high studentized-residual points have larger Cook’s distance: Sierre = 0.148, Rive Gauche = 0.112, Porrentruy = 0.208 (the largest in the data set). The high-leverage points do not influence the fitted model more than other points — their Cook’s distance is low because their residuals are small. The points with extreme studentized residuals do influence the model more than most other points, with Porrentruy the single most influential observation overall.

Part 2: Cheddar Cheese Data

library(faraway)   # load package (install.packages("faraway") first if needed)
data(cheddar)       # load the cheddar dataset into the environment
str(cheddar)        # preliminary look at structure

## 'data.frame':    30 obs. of  4 variables:
##  $ taste : num  12.3 20.9 39 47.9 5.6 25.9 37.3 21.9 18.1 21 ...
##  $ Acetic: num  4.54 5.16 5.37 5.76 4.66 ...
##  $ H2S   : num  3.13 5.04 5.44 7.5 3.81 ...
##  $ Lactic: num  0.86 1.53 1.57 1.81 0.99 1.09 1.29 1.78 1.29 1.58 ...

summary(cheddar)    # preliminary summary statistics

##      taste           Acetic           H2S             Lactic     
##  Min.   : 0.70   Min.   :4.477   Min.   : 2.996   Min.   :0.860  
##  1st Qu.:13.55   1st Qu.:5.237   1st Qu.: 3.978   1st Qu.:1.250  
##  Median :20.95   Median :5.425   Median : 5.329   Median :1.450  
##  Mean   :24.53   Mean   :5.498   Mean   : 5.942   Mean   :1.442  
##  3rd Qu.:36.70   3rd Qu.:5.883   3rd Qu.: 7.575   3rd Qu.:1.667  
##  Max.   :57.20   Max.   :6.458   Max.   :10.199   Max.   :2.010

pairs(cheddar)        # scatterplot matrix of all variables

cor(cheddar)           # correlation matrix

##            taste    Acetic       H2S    Lactic
## taste  1.0000000 0.5495393 0.7557523 0.7042362
## Acetic 0.5495393 1.0000000 0.6179559 0.6037826
## H2S    0.7557523 0.6179559 1.0000000 0.6448123
## Lactic 0.7042362 0.6037826 0.6448123 1.0000000

Taste ranges from 0.7 to 57.2 with mean 24.5. All three predictors are positively correlated with taste: H2S has the strongest association (r = 0.756), followed by Lactic (r = 0.704) and Acetic (r = 0.550). The predictors are also correlated with each other (r between 0.60 and 0.65), suggesting some multicollinearity.

cheddar_lm <- lm(taste ~ Acetic + H2S + Lactic, data = cheddar)  # fit full model
summary(cheddar_lm)

## 
## Call:
## lm(formula = taste ~ Acetic + H2S + Lactic, data = cheddar)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -17.390  -6.612  -1.009   4.908  25.449 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept) -28.8768    19.7354  -1.463  0.15540   
## Acetic        0.3277     4.4598   0.073  0.94198   
## H2S           3.9118     1.2484   3.133  0.00425 **
## Lactic       19.6705     8.6291   2.280  0.03108 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 10.13 on 26 degrees of freedom
## Multiple R-squared:  0.6518, Adjusted R-squared:  0.6116 
## F-statistic: 16.22 on 3 and 26 DF,  p-value: 3.81e-06

The model is significant overall (F = 16.22, p < 0.001) and explains 65.2% of the variance in taste (adjusted R² = 0.612). H2S (p = 0.004) and Lactic (p = 0.031) are significant predictors; Acetic (p = 0.942) is not, likely due to overlap with H2S and Lactic.

vif(cheddar_lm)   # variance inflation factors to check multicollinearity

##   Acetic      H2S   Lactic 
## 1.831589 1.992200 1.937912

VIFs are all below 2, so multicollinearity is not severe despite the pairwise correlations among predictors.

par(mfrow = c(2, 2))
plot(cheddar_lm)   # standard diagnostic plots: residuals vs fitted, QQ, scale-location, leverage

par(mfrow = c(1, 1))

shapiro.test(residuals(cheddar_lm))        # test normality of residuals

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(cheddar_lm)
## W = 0.98021, p-value = 0.8312

h_cheddar <- hatvalues(cheddar_lm)
p_cheddar <- length(coef(cheddar_lm))
n_cheddar <- nrow(cheddar)
which(h_cheddar > 2 * p_cheddar / n_cheddar)  # check for high-leverage points

## named integer(0)

cooks.distance(cheddar_lm)[order(cooks.distance(cheddar_lm), decreasing = TRUE)[1:3]]  # largest Cook's D

##        15        12        30 
## 0.1714653 0.1231599 0.1017165