# Generalized least squares results interpretation

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

• – Stefan Jan 27 '18 at 14:47
• According to your code you model Species ~ Temperature? If so, what type of variable is species? Also I wouldn't compare lm() vs gls(), instead I would compare a gls() model with and without the correlation structure. – Stefan Jan 27 '18 at 14:57
• Thank you for your answer it is helpful - I am comparing an east and a west location so species are recorded as the percentage seen in the east/west each month. Actual numbers cannot be used as there are no true zeros as the data was collected via 'citizen science' thus, a day with no sighting does not mean the animal is not there. I am trying to model the percentage at each location with various environmental variables such as SST and Chl-a to see if there is any link. – Jo Harris Jan 27 '18 at 16:09
• Maybe can you provide some made data that represents the data you are analyzing? See here for how: adv-r.had.co.nz/Reproducibility.html – Stefan Jan 27 '18 at 17:35
• Depending on the distribution of your percentages, a beta regression model may be more appropriate. How many years are you analyzing? Is there a temporal trend in the monthly data when you plot it? Do you analyze east and west separately? – Stefan Jan 27 '18 at 17:38

I have now overcome the problem of autocorrelation so I can use lm().

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

• Do you need help adding this to your original question? – EngrStudent Oct 24 '18 at 18:08