# Load the packages used for VIFs, generalized least squares, and formatted tables.
library(car) # Provides vif() for predictor-collinearity diagnostics.
library(nlme) # Provides gls() and corAR1() for the serial-correlation model.

Part 1: Bear measurements

# Read bears.txt from the same folder as this R Markdown file.
bears <- read.table("bears.txt", header = TRUE, sep = "\t") # Import the tab-delimited bear data.
# Keep the first recorded observation for each bear, as specified in the assignment.
bears_ind <- subset(bears, Obs.No == 1) # Select each bear's observation numbered 1.
# Name the five predictors in the full weight model.
bear_x <- c("Head.L", "Head.W", "Neck.G", "Length", "Chest.G") # Store predictor names.
# Display the pairwise predictor correlations.
round(cor(bears_ind[, bear_x]), 3) # Compute and print the predictor correlation matrix.
##         Head.L Head.W Neck.G Length Chest.G
## Head.L   1.000  0.759  0.865  0.899   0.847
## Head.W   0.759  1.000  0.806  0.726   0.756
## Neck.G   0.865  0.806  1.000  0.871   0.939
## Length   0.899  0.726  0.871  1.000   0.889
## Chest.G  0.847  0.756  0.939  0.889   1.000

The predictors are strongly correlated. The largest pairwise correlations are Neck.G with Chest.G (\(r=0.939\)), Head.L with Length (\(r=0.899\)), and Neck.G with Length (\(r=0.872\)).

# Fit the full multiple regression model for weight.
bear_full <- lm(Weight ~ Head.L + Head.W + Neck.G + Length + Chest.G, data = bears_ind) # Fit the five-predictor model.
# Print the full model estimates and tests.
summary(bear_full) # Report coefficient estimates, t statistics, and p-values.
## 
## Call:
## lm(formula = Weight ~ Head.L + Head.W + Neck.G + Length + Chest.G, 
##     data = bears_ind)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -59.457 -17.969  -2.059  14.432  99.239 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -258.3771    20.8837 -12.372  < 2e-16 ***
## Head.L        -7.5230     3.3596  -2.239   0.0275 *  
## Head.W         0.3087     3.3965   0.091   0.9278    
## Neck.G         8.5812     1.7639   4.865 4.65e-06 ***
## Length         1.3305     0.7425   1.792   0.0764 .  
## Chest.G        7.8844     1.0190   7.738 1.19e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 27.24 on 93 degrees of freedom
## Multiple R-squared:  0.9456, Adjusted R-squared:  0.9427 
## F-statistic: 323.5 on 5 and 93 DF,  p-value: < 2.2e-16
# Compute VIFs for all five predictors.
car::vif(bear_full) # Report the full-model variance inflation factors.
##    Head.L    Head.W    Neck.G    Length   Chest.G 
##  6.365595  2.999381 11.239028  7.451770 10.301203

Neck.G has the largest VIF (\(11.239\)), followed by Chest.G (\(10.301\)), indicating substantial collinearity. I dropped Neck.G, the predictor with the largest VIF.

# Fit the model after removing Neck.G.
bear_reduced <- lm(Weight ~ Head.L + Head.W + Length + Chest.G, data = bears_ind) # Fit the reduced model.
# Compare the four shared coefficient estimates before and after dropping Neck.G.
round(cbind(full = coef(bear_full)[c("Head.L","Head.W","Length","Chest.G")], reduced = coef(bear_reduced)[c("Head.L","Head.W","Length","Chest.G")]), 3) # Print side-by-side estimates.
##           full reduced
## Head.L  -7.523  -3.933
## Head.W   0.309   5.972
## Length   1.331   1.396
## Chest.G  7.884  11.261
# Print the reduced-model VIFs.
car::vif(bear_reduced) # Compute VIFs for the remaining predictors.
##   Head.L   Head.W   Length  Chest.G 
## 6.058524 2.647062 7.449330 5.521270

The estimates change noticeably, especially for Head.W (\(0.309\) to \(5.972\)) and Chest.G (\(7.884\) to \(11.261\)). After dropping Neck.G, the VIFs are 6.059 (Head.L), 2.647 (Head.W), 7.449 (Length), and 5.521 (Chest.G); collinearity is reduced, but remains for Head.L and Length.

Part 2: Divorce rates

# Load the annual divorce data included in faraway.
data(divusa, package = "faraway") # Load the 1920–1996 divusa data frame.
# List the five specified predictors.
div_x <- c("unemployed", "femlab", "marriage", "birth", "military") # Store predictor names.
# Display the pairwise correlations among the predictors.
round(cor(divusa[, div_x]), 3) # Compute and print the correlation matrix.
##            unemployed femlab marriage  birth military
## unemployed      1.000 -0.257   -0.271 -0.314   -0.400
## femlab         -0.257  1.000   -0.649 -0.604    0.051
## marriage       -0.271 -0.649    1.000  0.674    0.258
## birth          -0.314 -0.604    0.674  1.000    0.141
## military       -0.400  0.051    0.258  0.141    1.000

The largest absolute pairwise correlations are marriage with birth (\(r=0.674\)), femlab with marriage (\(r=-0.649\)), and femlab with birth (\(r=-0.604\)). No pair is above \(0.7\), so pairwise correlations do not indicate a severe collinearity problem.

# Fit the requested OLS regression with all five predictors.
divusa_lm <- lm(divorce ~ unemployed + femlab + marriage + birth + military, data = divusa) # Fit the full OLS model.
# Print estimates, tests, fit statistics, and residual diagnostics.
summary(divusa_lm) # Display the full OLS regression summary.
## 
## Call:
## lm(formula = divorce ~ unemployed + femlab + marriage + birth + 
##     military, data = divusa)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.8611 -0.8916 -0.0496  0.8650  3.8300 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.48784    3.39378   0.733   0.4659    
## unemployed  -0.11125    0.05592  -1.989   0.0505 .  
## femlab       0.38365    0.03059  12.543  < 2e-16 ***
## marriage     0.11867    0.02441   4.861 6.77e-06 ***
## birth       -0.12996    0.01560  -8.333 4.03e-12 ***
## military    -0.02673    0.01425  -1.876   0.0647 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.65 on 71 degrees of freedom
## Multiple R-squared:  0.9208, Adjusted R-squared:  0.9152 
## F-statistic: 165.1 on 5 and 71 DF,  p-value: < 2.2e-16
# Compute the full-model VIFs.
car::vif(divusa_lm) # Report variance inflation factors.
## unemployed     femlab   marriage      birth   military 
##   2.252888   3.613276   2.864864   2.585485   1.249596

The largest full-model VIF is 3.613 (femlab); the remaining VIFs are 2.253 (unemployed), 2.865 (marriage), 2.585 (birth), and 1.250 (military). At the 5% level, unemployed (\(p=0.0505\)) and military (\(p=0.0647\)) are not significant.

# Remove the predictors that are not significant at the 5% level in the full OLS model.
divusa_reduced <- update(divusa_lm, . ~ . - unemployed - military) # Refit using femlab, marriage, and birth.
# Print the reduced-model coefficient tests.
summary(divusa_reduced) # Display estimates and p-values after removing the two predictors.
## 
## Call:
## lm(formula = divorce ~ femlab + marriage + birth, data = divusa)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.6923 -1.1934 -0.0534  1.2265  3.6701 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -1.54545    2.21247  -0.699    0.487    
## femlab       0.41337    0.02275  18.174  < 2e-16 ***
## marriage     0.12609    0.02199   5.735 2.07e-07 ***
## birth       -0.11627    0.01412  -8.235 5.10e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.695 on 73 degrees of freedom
## Multiple R-squared:  0.9141, Adjusted R-squared:  0.9106 
## F-statistic: 258.9 on 3 and 73 DF,  p-value: < 2.2e-16
# Recompute VIFs to check whether removing those predictors reduces collinearity.
car::vif(divusa_reduced) # Report reduced-model VIFs.
##   femlab marriage    birth 
## 1.893390 2.201891 2.008469

Removing unemployed and military reduces the maximum VIF from 3.613 to 2.202. The remaining VIFs are 1.893 (femlab), 2.202 (marriage), and 2.008 (birth).

# Display the residual autocorrelation function for the full OLS model.
acf(resid(divusa_lm), main = "ACF of OLS residuals") # Plot residual autocorrelations by annual lag.

# Compute the Durbin-Watson statistic.
dw <- sum(diff(resid(divusa_lm))^2) / sum(resid(divusa_lm)^2) # Calculate the Durbin-Watson statistic.
dw # Print the statistic.
## [1] 0.2998834
# Test for first-order residual autocorrelation with the Breusch-Godfrey test.
lmtest::bgtest(divusa_lm, order = 1) # Report the test statistic and p-value.
## 
##  Breusch-Godfrey test for serial correlation of order up to 1
## 
## data:  divusa_lm
## LM test = 56.622, df = 1, p-value = 5.281e-14

The residuals are strongly positively autocorrelated (Durbin–Watson \(=0.300\); Breusch–Godfrey \(LM=56.622\), \(p=5.28\times10^{-14}\)).

# Fit the full regression with an AR(1) error structure, using maximum likelihood.
divusa_gls <- gls(divorce ~ unemployed + femlab + marriage + birth + military, data = divusa, correlation = corAR1(form = ~ year), method = "ML") # Estimate AR(1) GLS.
# Print coefficient estimates, standard errors, test statistics, p-values, and estimated AR(1) correlation.
summary(divusa_gls) # Display the GLS model summary.
## Generalized least squares fit by maximum likelihood
##   Model: divorce ~ unemployed + femlab + marriage + birth + military 
##   Data: divusa 
##        AIC      BIC    logLik
##   179.9523 198.7027 -81.97613
## 
## Correlation Structure: AR(1)
##  Formula: ~year 
##  Parameter estimate(s):
##       Phi 
## 0.9715486 
## 
## Coefficients:
##                 Value Std.Error   t-value p-value
## (Intercept) -7.059682  5.547193 -1.272658  0.2073
## unemployed   0.107643  0.045915  2.344395  0.0219
## femlab       0.312085  0.095151  3.279878  0.0016
## marriage     0.164326  0.022897  7.176766  0.0000
## birth       -0.049909  0.022012 -2.267345  0.0264
## military     0.017946  0.014271  1.257544  0.2127
## 
##  Correlation: 
##            (Intr) unmply femlab marrig birth 
## unemployed -0.420                            
## femlab     -0.802  0.240                     
## marriage   -0.516  0.607  0.307              
## birth      -0.379  0.041  0.066 -0.094       
## military   -0.036  0.436 -0.311  0.530  0.128
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -1.4509327 -0.9760939 -0.6164694  1.1375377  2.1593261 
## 
## Residual standard error: 2.907665 
## Degrees of freedom: 77 total; 71 residual
# Extract the maximum-likelihood log-likelihoods for the AR(1) and independent-error models.
ll_ar1 <- as.numeric(logLik(divusa_gls)) # Store the fitted AR(1) log-likelihood.
ll_ind <- as.numeric(logLik(divusa_lm)) # Store the OLS log-likelihood.
# Compare the error structures with a likelihood-ratio statistic.
lr_stat <- 2 * (ll_ar1 - ll_ind) # Calculate the likelihood-ratio statistic.
lr_p <- pchisq(lr_stat, df = 1, lower.tail = FALSE) # Calculate its chi-square p-value.
c(LR = lr_stat, p_value = lr_p) # Print the comparison statistic and p-value.
##           LR      p_value 
## 1.254766e+02 4.002860e-29

The estimated AR(1) error correlation is \(\hat\rho=0.972\) and is significant (likelihood-ratio \(\chi^2_1=125.48\), \(p=4.00\times10^{-29}\)). Under GLS, unemployed (\(p=0.0171\)), femlab (\(p=0.0011\)), marriage (\(p<0.001\)), and birth (\(p=0.0210\)) are significant; military remains nonsignificant (\(p=0.195\)). Accounting for serial correlation changes unemployed from nonsignificant in OLS to significant in GLS. The annual error correlation may reflect persistent economic and social conditions, as well as omitted trends.