require("faraway")
## Loading required package: faraway
## Warning: package 'faraway' was built under R version 3.5.3
library("nlme")
library("ggplot2")
library("stats")
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
d<- as.data.frame(lapply(divusa[,2:7], diff, lag=1))
d$year <- tail(divusa$year,-1)
plot.ts(divusa[,2:7],main="Original")
plot.ts(d[,1:6],main="Differenced")
# Fit model 1
l_mod <- lm(divorce ~ .-year, data = d)
summary(l_mod)
##
## Call:
## lm(formula = divorce ~ . - year, data = d)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.79611 -0.33769 -0.03295 0.36598 1.85120
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.24215 0.11536 2.099 0.03942 *
## unemployed 0.09578 0.04436 2.159 0.03428 *
## femlab -0.05093 0.16738 -0.304 0.76183
## marriage 0.15735 0.02220 7.088 8.63e-10 ***
## birth -0.02825 0.02275 -1.241 0.21857
## military 0.04792 0.01788 2.680 0.00917 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.6971 on 70 degrees of freedom
## Multiple R-squared: 0.4889, Adjusted R-squared: 0.4524
## F-statistic: 13.39 on 5 and 70 DF, p-value: 3.679e-09
cor(residuals(l_mod)[-1],residuals(l_mod)[-length(residuals(l_mod))]) # Correlation between consecutive residuals
## [1] 0.4260145
# Fit model 2 (Generalized Least Squares)
gl_mod <- gls(divorce ~ .-year, correlation = corAR1(), method="ML", data = na.omit(d))
summary(gl_mod)
## Generalized least squares fit by maximum likelihood
## Model: divorce ~ . - year
## Data: na.omit(d)
## AIC BIC logLik
## 151.4896 170.1355 -67.74482
##
## Correlation Structure: AR(1)
## Formula: ~1
## Parameter estimate(s):
## Phi
## 0.508056
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 0.31418746 0.15899758 1.976052 0.0521
## unemployed 0.10163537 0.04030327 2.521765 0.0140
## femlab -0.20067958 0.14522541 -1.381849 0.1714
## marriage 0.14485763 0.02056861 7.042657 0.0000
## birth -0.00235827 0.02169264 -0.108713 0.9137
## military 0.04921917 0.01795131 2.741815 0.0078
##
## Correlation:
## (Intr) unmply femlab marrig birth
## unemployed -0.072
## femlab -0.446 0.211
## marriage -0.008 0.623 0.177
## birth 0.208 -0.141 -0.276 -0.035
## military 0.339 0.267 -0.630 0.500 0.349
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -3.14261680 -0.53905878 0.02802262 0.44735606 2.63066153
##
## Residual standard error: 0.6836979
## Degrees of freedom: 76 total; 70 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.2783656 0.508056 0.6827794
## attr(,"label")
## [1] "Correlation structure:"
##
## Residual standard error:
## lower est. upper
## 0.5540987 0.6836979 0.8436093