Part 1: Bears Data

Setup

# Read in the bears dataset
bears <- read.table("bears.txt", header = TRUE)

# Keep only the first observation for each bear so observations are independent
bears_indep <- bears[bears$Obs.No == 1, ]
nrow(bears_indep)
## [1] 99

Model 1: Weight ~ Length

m1 <- lm(Weight ~ Length, data = bears_indep)
summary(m1)
## 
## Call:
## lm(formula = Weight ~ Length, data = bears_indep)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -100.88  -36.19  -12.22   30.29  183.40 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -408.0766    33.8779  -12.05   <2e-16 ***
## Length         9.8490     0.5517   17.85   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 55.26 on 97 degrees of freedom
## Multiple R-squared:  0.7666, Adjusted R-squared:  0.7642 
## F-statistic: 318.7 on 1 and 97 DF,  p-value: < 2.2e-16
# Residuals vs fitted values
plot(fitted(m1), resid(m1),
     xlab = "Fitted Values", ylab = "Residuals",
     main = "Residuals vs Fitted: Weight ~ Length")
abline(h = 0, col = "red", lty = 2)

Linearity: There is a slight curved pattern in the residuals (negative in the low and high fitted ranges, positive in the middle), suggesting a mild violation of linearity. Adding a quadratic term (Length^2) or transforming the predictor could improve the fit.

Equivariance: The spread of residuals increases as fitted values increase, so there’s evidence against constant variance. Transforming the response (e.g., log(Weight)) should help stabilize the variance.

Model 2: log(Weight) ~ Length

m2 <- lm(log(Weight) ~ Length, data = bears_indep)
summary(m2)
## 
## Call:
## lm(formula = log(Weight) ~ Length, data = bears_indep)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.65091 -0.15842 -0.01098  0.11902  0.62023 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 1.394436   0.139293   10.01   <2e-16 ***
## Length      0.060317   0.002269   26.59   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2272 on 97 degrees of freedom
## Multiple R-squared:  0.8793, Adjusted R-squared:  0.8781 
## F-statistic: 706.9 on 1 and 97 DF,  p-value: < 2.2e-16
# Residuals vs fitted values
plot(fitted(m2), resid(m2),
     xlab = "Fitted Values", ylab = "Residuals",
     main = "Residuals vs Fitted: log(Weight) ~ Length")
abline(h = 0, col = "red", lty = 2)

Linearity: The residuals are scattered more randomly around zero with no obvious curved pattern, so linearity looks reasonable here.

Equivariance: The spread of residuals is roughly constant across the range of fitted values. The log transformation fixed the funnel shape seen in Model 1.

Q-Q Plot and Shapiro-Wilk Test (log model)

qqnorm(resid(m2), main = "Q-Q Plot: log(Weight) ~ Length")
qqline(resid(m2), col = "red")

shapiro.test(resid(m2))
## 
##  Shapiro-Wilk normality test
## 
## data:  resid(m2)
## W = 0.98805, p-value = 0.519

The Q-Q plot shows the residuals falling close to the reference line with no major departures in the tails. The Shapiro-Wilk test gives \(W = 0.988\), \(p = 0.519\). Since \(p > 0.05\), I fail to reject normality — the residuals from the log-transformed model are consistent with a normal distribution.

Part 2: Teengamb Data (faraway package)

Setup and Model

library(faraway)
data(teengamb)

# Fit full model: gamble as response, all other variables as predictors
gamb_model <- lm(gamble ~ sex + status + income + verbal, data = teengamb)
summary(gamb_model)
## 
## Call:
## lm(formula = gamble ~ sex + status + income + verbal, data = teengamb)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -51.082 -11.320  -1.451   9.452  94.252 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  22.55565   17.19680   1.312   0.1968    
## sex         -22.11833    8.21111  -2.694   0.0101 *  
## status        0.05223    0.28111   0.186   0.8535    
## income        4.96198    1.02539   4.839 1.79e-05 ***
## verbal       -2.95949    2.17215  -1.362   0.1803    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 22.69 on 42 degrees of freedom
## Multiple R-squared:  0.5267, Adjusted R-squared:  0.4816 
## F-statistic: 11.69 on 4 and 42 DF,  p-value: 1.815e-06

Normality Check

# Histogram of residuals
hist(resid(gamb_model), breaks = 10,
     main = "Histogram of Residuals", xlab = "Residuals")

# Q-Q plot
qqnorm(resid(gamb_model), main = "Q-Q Plot of Residuals")
qqline(resid(gamb_model), col = "red")

# Shapiro-Wilk test
shapiro.test(resid(gamb_model))
## 
##  Shapiro-Wilk normality test
## 
## data:  resid(gamb_model)
## W = 0.86839, p-value = 8.16e-05

The histogram is right-skewed rather than symmetric, and the Q-Q plot shows the residuals deviating from the reference line in the upper tail (large positive residuals). The Shapiro-Wilk test confirms this: \(p < 0.05\), so I reject the null hypothesis of normality. The normality assumption is violated for this model.

Equivariance and Nonlinearity Check

plot(fitted(gamb_model), resid(gamb_model),
     xlab = "Fitted Values", ylab = "Residuals",
     main = "Residuals vs Fitted: Teengamb Model")
abline(h = 0, col = "red", lty = 2)

The residuals fan out as the fitted values increase, showing a clear violation of the equivariance (constant variance) assumption. There’s also one point (observation 24) with a much larger residual than the rest, standing out as a potential outlier. There’s no strong evidence of nonlinearity — the issue is primarily non-constant variance and the one extreme observation, not curvature.