GAMM with spatial auto-correlation in R

Question

From what I understand, GAMM has a random error or a spatial autocorrelation error structure.

I am trying to run a GAMM model with a spatial autocorrelation error structure such as the corExp (see https://stat.ethz.ch/R-manual/R-devel/library/nlme/html/corClasses.html).

But I am getting confused modelling GAMM from GAM.

Below is my code for GAM (Generalised additive model).

library(mgcv)
mod.gam <- gam(Chave~s(Band_3) + s(Band_7) + s(Band_8) + s(BaCo_3_2) + s(BaCo_5_2) + 
  s(BaCo_5_3) + s(BaCo_5_4) + s(BaCo_8_2) + s(BaCo_8A_2),data = data)
summary(mod.gam)
fits = predict(mod.gam, newdata=data, type='response')
plot(data$Chave, fits, col='blue', ylab = "Predicted AGB KNN", xlab = "Estimated AGB")
plot(mod.gam)

How do I change this to GAMM? When I try to use GAMM (as mentioned below), I do not see much of a change:

model <- gamm(Chave~s(Band_3) + s(Band_7) + s(Band_8) + s(BaCo_3_2) +
 s(BaCo_5_2) + s(BaCo_5_3) + s(BaCo_5_4) + s(BaCo_8_2) + s(BaCo_8A_2),
 data = data, correlation=corExp(1,form = ~ Latitude + Longitude)))
summary(model$gam)
plot(model$gam,pages=1)

From what I understand the spatial auto-correlation function should have latitude and longitude as input. Am I doing it correctly? And how is the output of GAMM different from the output of GAM?

Output of different models:

model.gamm <- gamm(Chave~s(Band_3) + s(Band_7) + s(Band_8) + s(BaCo_3_2) + s(BaCo_5_2) + 
       s(BaCo_5_3) + s(BaCo_5_4) + s(BaCo_8_2) + s(BaCo_8A_2), data = data, 
       correlation=corSpher(form = ~ Latitude + Longitude))
summary(model.gam$lme)

[Output]:
Linear mixed-effects model fit by maximum likelihood
 Data: strip.offset(mf) 
       AIC      BIC    logLik
  1176.797 1232.729 -567.3985
Correlation Structure: Spherical spatial correlation
 Formula: ~Latitude + Longitude | g/g.0/g.1/g.2/g.3/g.4/g.5/g.6/g.7 
 Parameter estimate(s):
      range 
0.003845795 
Fixed effects: y ~ X - 1 
                      Value Std.Error DF   t-value p-value
X(Intercept)      125.44208   5.23505 96 23.961953  0.0000
Xs(Band_3)Fx1      95.48991  38.42171 96  2.485311  0.0147
Xs(Band_7)Fx1     -17.69103  40.27462 96 -0.439260  0.6615
Xs(Band_8)Fx1     -84.45737  60.83131 96 -1.388386  0.1682
Xs(BaCo_3_2)Fx1  -183.67952 179.22609 96 -1.024848  0.3080
Xs(BaCo_5_2)Fx1   111.84971 266.21261 96  0.420152  0.6753
Xs(BaCo_5_3)Fx1   -46.55365 131.32239 96 -0.354499  0.7237
Xs(BaCo_5_4)Fx1   -10.57644   9.59420 96 -1.102378  0.2731
Xs(BaCo_8_2)Fx1    76.30654  52.36900 96  1.457094  0.1484
Xs(BaCo_8A_2)Fx1   35.37155  39.34172 96  0.899085  0.3709

#########################

model.gam <- gamm(Chave~s(Band_3) + s(Band_7) + s(Band_8) + s(BaCo_3_2) + 
s(BaCo_5_2) + s(BaCo_5_3) + s(BaCo_5_4) + s(BaCo_8_2) + s(BaCo_8A_2), data = data)
summary(model.gamm$lme)

[Output]:
Linear mixed-effects model fit by maximum likelihood
 Data: strip.offset(mf) 
       AIC      BIC   logLik
  1177.106 1230.375 -568.553
Fixed effects: y ~ X - 1 
                      Value Std.Error DF   t-value p-value
X(Intercept)      123.71814   4.94795 96 25.003936  0.0000
Xs(Band_3)Fx1     102.63429  38.76940 96  2.647301  0.0095
Xs(Band_7)Fx1     -26.65351  42.00784 96 -0.634489  0.5273
Xs(Band_8)Fx1     -81.96547  61.57940 96 -1.331053  0.1863
Xs(BaCo_3_2)Fx1  -312.47086 191.52784 96 -1.631464  0.1061
Xs(BaCo_5_2)Fx1   283.77313 280.17983 96  1.012825  0.3137
Xs(BaCo_5_3)Fx1  -126.82534 139.08953 96 -0.911825  0.3641
Xs(BaCo_5_4)Fx1   -11.82504   9.69256 96 -1.220012  0.2254
Xs(BaCo_8_2)Fx1    78.95786  52.86484 96  1.493580  0.1386
Xs(BaCo_8A_2)Fx1   38.03149  40.01424 96  0.950449  0.3443

##############################

mod.gam.gam <- gam(Chave~s(Band_3) + s(Band_7) + s(Band_8) + s(BaCo_3_2) + s(BaCo_5_2) + 
+   s(BaCo_5_3) + s(BaCo_5_4) + s(BaCo_8_2) + s(BaCo_8A_2),data = data)
summary(mod.gam.gam)

[Output]:
Family: gaussian 
Link function: identity 

Formula:
Chave ~ s(Band_3) + s(Band_7) + s(Band_8) + s(BaCo_3_2) + s(BaCo_5_2) + 
    s(BaCo_5_3) + s(BaCo_5_4) + s(BaCo_8_2) + s(BaCo_8A_2)

Parametric coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  123.718      4.155   29.78   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
               edf Ref.df     F p-value  
s(Band_3)    1.000  1.000 1.125  0.2922  
s(Band_7)    5.854  6.945 2.419  0.0268 *
s(Band_8)    4.130  5.210 2.989  0.0143 *
s(BaCo_3_2)  4.724  5.694 3.128  0.0113 *
s(BaCo_5_2)  2.628  3.431 1.449  0.2275  
s(BaCo_5_3)  2.614  3.262 2.713  0.0530 .
s(BaCo_5_4)  4.441  5.454 1.045  0.4176  
s(BaCo_8_2)  1.000  1.000 1.144  0.2881  
s(BaCo_8A_2) 3.373  4.250 1.077  0.3394  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =  0.445   Deviance explained = 60.2%
GCV = 2578.1  Scale est. = 1829.9    n = 106

Answer 1

You are, I think , calling corExp() incorrectly. You use:

corExp(1, form = ~ Latitude + Longitude)

which is fixing the value of the correlation parameter in the exponential correlation function to be a fixed value of 1 rather than be estimated from the data, which would be achieved by instead using

corExp(form = ~ Latitude + Longitude)

You may also not see much of a change in the model if the residual spatial correlation from the GAM (gam()) was orthogonal or not non-linearly correlated with the smooth covariates: the correlation functions enter the model through the covariance structure of the response, not through the mean structure.

If you want to look at residuals from the gamm() that include the correlation function effect, then you need to use:

resid(model$lme, type = "normalized")

Also, if you want to compare a model with and without the structure, you'll need to use anova() on the $lme components of the models, hence you'll need to refit the gam() model using gamm(), just leave off the correlation part. Then the anova(mod1$lme, mod2$lme) will give you a generalised likelihood ratio test of the parameter in the exponential correlation function.