3
$\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
$\endgroup$
  • $\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
0
$\begingroup$

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$
  • $\begingroup$ Do you need help adding this to your original question? $\endgroup$ – EngrStudent Oct 24 '18 at 18:08

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.