4
$\begingroup$

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 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

$\endgroup$
6
  • $\begingroup$ See here: stats.stackexchange.com/questions/13859/… $\endgroup$
    – Stefan
    Jan 27 '18 at 14:47
  • 1
    $\begingroup$ 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. $\endgroup$
    – Stefan
    Jan 27 '18 at 14:57
  • $\begingroup$ 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. $\endgroup$
    – Jo Harris
    Jan 27 '18 at 16:09
  • $\begingroup$ 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 $\endgroup$
    – Stefan
    Jan 27 '18 at 17:35
  • $\begingroup$ 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? $\endgroup$
    – Stefan
    Jan 27 '18 at 17:38

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.