require(astsa)
## Loading required package: astsa
require(knitr)
## Loading required package: knitr
mort = cbind(Mortality=cmort, Temperature=tempr, Particulate=part)
plot.ts(mort)
pairs(mort)
trend=time(cmort)
temp=tempr - mean(tempr)
temp2 = temp^2
temp3 = temp^3
temp4 = temp^4
part2 = part^2
part3 = part^3
fit1 = lm( cmort ~ trend )
fit2 = lm( cmort ~ trend + temp )
fit3 = lm( cmort ~ trend + temp + part)
fit4 = lm( cmort ~ trend + temp + part + temp2)
fit5 = lm( cmort ~ trend + temp + part + temp2 + temp3 )
fit6 = lm( cmort ~ trend + temp + part + temp2 + temp3 + temp4 )
fit7 = lm( cmort ~ trend + temp + part + temp2 + temp3 + temp4 + part2 )
fit8 = lm( cmort ~ trend + temp + part + temp2 + temp3 + temp4 + part2 + part3)
AIC(fit1, fit2, fit3, fit4, fit5, fit6, fit7, fit8)
## df AIC
## fit1 3 3665.899
## fit2 4 3544.889
## fit3 5 3390.986
## fit4 6 3332.282
## fit5 7 3327.014
## fit6 8 3329.012
## fit7 9 3329.865
## fit8 10 3330.168
BIC(fit1, fit2, fit3, fit4, fit5, fit6, fit7, fit8)
## df BIC
## fit1 3 3678.591
## fit2 4 3561.811
## fit3 5 3412.138
## fit4 6 3357.664
## fit5 7 3356.627
## fit6 8 3362.856
## fit7 9 3367.939
## fit8 10 3372.473
Both AIC and BIC suggest that the preferred model is fit5, which includes squared and cubed temperature.
```