HHernandez_DATA621_Discussion11

LMR Exercise 8.2

require("faraway")

data("divusa")
summary(divusa)

##       year         divorce        unemployed         femlab     
##  Min.   :1920   Min.   : 6.10   Min.   : 1.200   Min.   :22.70  
##  1st Qu.:1939   1st Qu.: 8.70   1st Qu.: 4.200   1st Qu.:27.47  
##  Median :1958   Median :10.60   Median : 5.600   Median :37.10  
##  Mean   :1958   Mean   :13.27   Mean   : 7.173   Mean   :38.58  
##  3rd Qu.:1977   3rd Qu.:20.30   3rd Qu.: 7.500   3rd Qu.:47.80  
##  Max.   :1996   Max.   :22.80   Max.   :24.900   Max.   :59.30  
##     marriage          birth           military     
##  Min.   : 49.70   Min.   : 65.30   Min.   : 1.940  
##  1st Qu.: 61.90   1st Qu.: 68.90   1st Qu.: 3.469  
##  Median : 74.10   Median : 85.90   Median : 9.102  
##  Mean   : 72.97   Mean   : 88.89   Mean   :12.365  
##  3rd Qu.: 80.00   3rd Qu.:107.30   3rd Qu.:14.266  
##  Max.   :118.10   Max.   :122.90   Max.   :86.641

# Fit model 1
l_mod <- lm(divorce ~ .-year, data = divusa)
summary(l_mod)

## 
## Call:
## lm(formula = divorce ~ . - year, 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

cor(residuals(l_mod)[-1],residuals(l_mod)[-length(residuals(l_mod))]) # Correlation between consecutive residuals

## [1] 0.8469792

plot(divorce~year, divusa)

library(nlme)
# Fit model 2 (Generalized Least Squares)
gl_mod <- gls(divorce ~ .-year, correlation = corAR1(), method="ML", data = na.omit(divusa))
summary(gl_mod)

## Generalized least squares fit by maximum likelihood
##   Model: divorce ~ . - year 
##   Data: na.omit(divusa) 
##        AIC      BIC    logLik
##   179.9523 198.7027 -81.97613
## 
## Correlation Structure: AR(1)
##  Formula: ~1 
##  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

plot(ACF(gl_mod, form =~Year), alpha=0.05)

intervals(gl_mod,which="var-cov")

## Approximate 95% confidence intervals
## 
##  Correlation structure:
##         lower      est.     upper
## Phi 0.6528156 0.9715486 0.9980191
## attr(,"label")
## [1] "Correlation structure:"
## 
##  Residual standard error:
##      lower       est.      upper 
##  0.7974459  2.9076645 10.6019895