I checked my linear regression model (WMAN = Species, WDNE = sea surface temp) and found auto-correlation so instead, I am trying generalized least squares with the following script;
library(nlme) modelwa <- gls(WMAN ~WDNE, data=dat, correlation = corAR1(form=~MONTH), na.action=na.omit) summary(modelwa)
I compared both models;
> library(MuMIn) > model.sel(modelw,modelwa) Model selection table (Intrc) WDNE class na.action correlation df logLik AICc delta modelwa 31.50 0.1874 gls na.omit crAR1(MONTH) 4 -610.461 1229.2 0.00 modelw 11.31 0.7974 lm na.excl 3 -658.741 1323.7 94.44 weight modelwa 1 modelw 0 Abbreviations: na.action: na.excl = ‘na.exclude’ correlation: crAR1(MONTH) = ‘corAR1(~MONTH)’ Models ranked by AICc(x)
I believe the results suggest I should use gls as the AIC is lower.
My problem is, I have been reporting F-value/R²/p-value, but the output from the gls does not have these?
I would be very grateful if someone could assist me in interpreting these results?
> summary(modelwa) Generalized least squares fit by REML Model: WMAN ~ WDNE Data: mp2017.dat AIC BIC logLik 1228.923 1240.661 -610.4614 Correlation Structure: ARMA(1,0) Formula: ~MONTH Parameter estimate(s): Phi1 0.4809973 Coefficients: Value Std.Error t-value p-value (Intercept) 31.496911 8.052339 3.911524 0.0001 WDNE 0.187419 0.091495 2.048401 0.0424 Correlation: (Intr) WDNE -0.339 Standardized residuals: Min Q1 Med Q3 Max -2.023362 -1.606329 -1.210127 1.427247 3.567186 Residual standard error: 18.85341 Degrees of freedom: 141 total; 139 residual
I have now overcome the problem of autocorrelation so I can use
Add lag1 of residual as an X variable to the original model. This can be done using the
slide() function in DataCombine package.
library(DataCombine) econ_data <- data.frame(economics, resid_mod1=lmMod$residuals) econ_data_1 <- slide(econ_data, Var="resid_mod1", NewVar = "lag1", slideBy = -1) econ_data_2 <- na.omit(econ_data_1) lmMod2 <- lm(pce ~ pop + lag1, data=econ_data_2)
This script can be found here