2
$\begingroup$

I am working through a spatial GAM example in one of Simon Wood's lectures (found here), and ran into some problems surrounding concurvity and penalization. I've read a number of other posts on concurvity in GAMs (here, here, and here), and a few on penalization and concurvity (here), but none of them appeared to address this particular problem, namely:

  • The model appears to fit well, and penalizes less-important terms out of the model, but...
  • Concurvity in the parameters causes very uncertain estimates of the smoothing parameters, possibly causing model to fail "behind the scenes" without notifying the user

Here's the model that's being fit:

library(gamair) #Load data
data(mack)
mack$log.net.area <- log(mack$net.area) #Log-transform net area

#Fit model of egg count, using s(lon,lat) to model spatial variation
gm <- gam(egg.count~s(lon,lat,k=100)+s(I(b.depth^.5))+
             s(c.dist)+s(temp.surf)+ 
            s(salinity)+s(temp.20m)+offset(log.net.area),
          data=mack,family=quasipoisson,method="REML",
          select=T) #Penalizes null space to remove terms (double-penalization)

The output from the model looks reasonable (output from gam.check looks fine too). The extra penalization has effectively removed s(temp.surf) and s(c.dist) from the model:

Family: quasipoisson 
Link function: log 

Formula:
egg.count ~ s(lon, lat, k = 100) + s(I(b.depth^0.5)) + s(c.dist) + 
    s(temp.surf) + s(salinity) + s(temp.20m) + offset(log.net.area)

Parametric coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   2.9910     0.1233   24.26   <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(lon,lat)        6.060e+01     99 3.400  < 2e-16 ***
s(I(b.depth^0.5)) 2.165e+00      9 0.908  0.00218 ** 
s(c.dist)         4.205e-01      9 0.062  0.04797 *  
s(temp.surf)      2.526e-04      9 0.000  0.43865    
s(salinity)       1.830e+00      9 0.913  0.00117 ** 
s(temp.20m)       5.173e+00      9 2.771 1.08e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =  0.828   Deviance explained = 90.2%
-REML = 515.08  Scale est. = 4.8428    n = 330

However, there are serious problems with concurvity; space (lat/lon smoother) is highly "correlated" with the other predictors (similar to collinearity in a GLM):

round(concurvity(gm,full=F)$worst)
                  para s(lon,lat) s(I(b.depth^0.5)) s(c.dist) s(temp.surf) s(salinity) s(temp.20m)
para                 1      0.000             0.000     0.000        0.000       0.000       0.000
s(lon,lat)           0      1.000             0.985     0.998        0.976       0.921       0.957
s(I(b.depth^0.5))    0      0.985             1.000     0.420        0.489       0.678       0.443
s(c.dist)            0      0.998             0.420     1.000        0.220       0.320       0.157
s(temp.surf)         0      0.976             0.489     0.220        1.000       0.682       1.000
s(salinity)          0      0.921             0.678     0.320        0.682       1.000       0.666
s(temp.20m)          0      0.957             0.443     0.157        1.000       0.666       1.000

On pg.3 of his lecture, Wood states that if the predictors are highly concurved, "smoothing parameter estimates may become highly correlated and variable, which degrades the performance of... confidence intervals and p-values", and that "Under ML or REML smoothness selection sp.vcov and gam.vcomp can help diagnose [the consequences of concurvity]".

After following his advice, I can see that there is some correlation between smoothing parameter estimates, but it doesn't look too bad:

    #Extract smoothing parameter covariance matrix and convert to correlation matrix
    round(cov2cor(sp.vcov(gm)),3) 
        [,1]   [,2]   [,3]   [,4]   [,5]   [,6]   [,7]   [,8]   [,9]  [,10]  [,11]  [,12]  [,13]
 [1,]  1.000  0.001 -0.068  0.000  0.003 -0.233  0.000 -0.001  0.163  0.000  0.163 -0.001  0.601
 [2,]  0.001  1.000 -0.002 -0.153 -0.089  0.002 -0.169 -0.180 -0.001 -0.147 -0.002 -0.197 -0.001
 [3,] -0.068 -0.002  1.000 -0.011  0.010  0.243 -0.001 -0.001  0.032 -0.001 -0.101 -0.002  0.028
 [4,]  0.000 -0.153 -0.011  1.000 -0.092  0.008 -0.175 -0.187 -0.002 -0.153 -0.001 -0.204  0.000
 [5,]  0.003 -0.089  0.010 -0.092  1.000 -0.065 -0.102 -0.108  0.005 -0.089  0.007 -0.118  0.008
 [6,] -0.233  0.002  0.243  0.008 -0.065  1.000  0.001  0.002  0.029  0.002 -0.086  0.003 -0.089
 [7,]  0.000 -0.169 -0.001 -0.175 -0.102  0.001  1.000 -0.206 -0.002 -0.169 -0.003 -0.226 -0.001
 [8,] -0.001 -0.180 -0.001 -0.187 -0.108  0.002 -0.206  1.000 -0.001 -0.180 -0.002 -0.241 -0.002
 [9,]  0.163 -0.001  0.032 -0.002  0.005  0.029 -0.002 -0.001  1.000 -0.005  0.028 -0.002  0.178
[10,]  0.000 -0.147 -0.001 -0.153 -0.089  0.002 -0.169 -0.180 -0.005  1.000 -0.001 -0.197 -0.002
[11,]  0.163 -0.002 -0.101 -0.001  0.007 -0.086 -0.003 -0.002  0.028 -0.001  1.000  0.000  0.245
[12,] -0.001 -0.197 -0.002 -0.204 -0.118  0.003 -0.226 -0.241 -0.002 -0.197  0.000  1.000 -0.002
[13,]  0.601 -0.001  0.028  0.000  0.008 -0.089 -0.001 -0.002  0.178 -0.002  0.245 -0.002  1.000

The variance components give a hint that something might be wrong; the confidence intervals are enormous:

#Extract variance components for smoothing parameter (1) and penalization term (2)
gam.vcomp(gm)
Standard deviations and 0.95 confidence intervals:

                       std.dev        lower        upper
s(lon,lat)1        5.222450031 3.972972e+00 6.864882e+00
s(lon,lat)2        0.005721280 9.588122e-75 3.413917e+69
s(I(b.depth^0.5))1 0.007041528 2.327681e-03 2.130151e-02
s(I(b.depth^0.5))2 0.006402433 5.447190e-77 7.525192e+71
s(c.dist)1         0.030582323 2.337054e-47 4.001954e+43
s(c.dist)2         0.213280517 6.130439e-03 7.420118e+00
s(temp.surf)1      0.001336795 3.747653e-84 4.768370e+77
s(temp.surf)2      0.002261587 2.343874e-88 2.182188e+82
s(salinity)1       2.911410660 7.610990e-01 1.113694e+01
s(salinity)2       0.004252481 8.852023e-75 2.042877e+69
s(temp.20m)1       1.321205095 4.819389e-01 3.622001e+00
s(temp.20m)2       0.011757738 5.143517e-94 2.687740e+89
scale              2.083335131 1.888021e+00 2.298854e+00

Rank: 13/13

Taking a closer look at the smoothing parameters and penalization terms, they also have enormous standard errors:

#Get smoothing terms
data.frame(terms=c(names(gm$sp),'scale'), #Names of terms
           log.sp=c(unname(log(gm$sp)),log(gmsig2)), #Log smoothing params
           log.se.sp=sqrt(diag(sp.vcov(gm)))) #Log smoothing param SE (sqrt of trace of Hessian)
                terms    log.sp   log.se.sp
1         s(lon,lat)1 -2.367212   0.3275753
2         s(lon,lat)2 11.795065 157.4245705
3  s(I(b.depth^0.5))1  1.884133   1.1278863
4  s(I(b.depth^0.5))2 11.570095 161.6735859
5          s(c.dist)1  8.955096 103.2894284
6          s(c.dist)2  4.558234   3.6113969
7       s(temp.surf)1 11.553715 173.2620884
8       s(temp.surf)2 13.651317 180.7825516
9        s(salinity)1  1.268129   1.3833522
10       s(salinity)2 12.388446 157.2485245
11       s(temp.20m)1 -2.400153   1.0490736
12       s(temp.20m)2 10.354428 190.9814056
13              scale  1.556048   0.1004503

My questions are:

  1. Does this indicate that the smoothing or penalization (or both) has failed? Is this model still interpretable, or are some of the terms inestimable?
  2. Since the SE estimates for the smoothing parameters are large, are the p-values from summary reliable? It seems like these variables should be discarded from the model, not just penalized out, since the estimates seem unreliable.
  3. How should I interpret variance components for penalty terms (e.g. s(temp.20m)2)? It doesn't seem like they're comparable to the variance components of the smoothing terms.
$\endgroup$

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.