一般化加法モデルのVIFの算出方法をメモしておきます。

加法モデルについてはこちらを参考に。

今回はサンプルデータとして、Rに組み込みのニューヨークの大気データを使用します。

data.airquality

aq<-airquality
head(aq,10)


   Ozone Solar.R Wind Temp Month Day
1     41     190  7.4   67     5   1
2     36     118  8.0   72     5   2
3     12     149 12.6   74     5   3
4     18     313 11.5   62     5   4
5     NA      NA 14.3   56     5   5
6     28      NA 14.9   66     5   6
7     23     299  8.6   65     5   7
8     19      99 13.8   59     5   8
9      8      19 20.1   61     5   9
10    NA     194  8.6   69     5  10

5月から9月まで毎日データがあるので、とりあえず日付を5月1日からの連番にします。
ところどころNAがありますが、ここではとりあえず削除します。
列名を簡単な表記に変更しておきます。

data.airquality

library(tidyverse)
aq<-aq %>%
  mutate(date = row_number()) %>% 
  na.omit()
aq=data.frame(ozone=aq$Ozone,solar=aq$Solar.R,wind=aq$Wind,temp=aq$Temp,date=aq$date)
plot(aq)

オゾンと関係ありそうなのは風速と気温でしょうか。
日射量と日付は微妙なのでとりあえずスプライン化してみます。
また気温と風速、気温と日付はそこそこ相関がありそうです。

準備ができたので早速一般化加法モデルを実行します。

gam

library(mgcv)
fit=gam(ozone~s(solar)+wind+temp+s(date),data=aq)
summary(fit)

Family: gaussian 
Link function: identity 

Formula:
ozone ~ s(solar) + wind + temp + s(date)

Parametric coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept) -67.4562    24.6202  -2.740  0.00724
wind         -3.1675     0.6360  -4.980 2.56e-06
temp          1.8130     0.2729   6.642 1.48e-09

Approximate significance of smooth terms:
           edf Ref.df     F p-value
s(solar) 3.717  4.618 2.521  0.0410
s(date)  1.000  1.000 6.377  0.0131

R-sq.(adj) =  0.624   Deviance explained = 64.7%
GCV = 447.77  Scale est. = 416.64    n = 111

lmやglmで使えるcarパッケージのvif関数をgamでもやってみます。

vif

library(car)
vif(fit)

            GVIF Df GVIF^(1/(2*Df))
wind  1344018355  0             Inf
temp  1344018355  0             Inf
solar 1344018355  0             Inf
date  1344018355  0             Inf

計算できません。
mgcv.helperなるパッケージがあるようですが、このバージョンではインストールできずのメッセージ。

さらに調べてみるとgithubにvifを計算するコードが公開されてました！
リンク

vif

samvif <- function(mod){
  # mod is an mgcv object
  # this function calculates the variance inflation factors for GAM as no one else has written code to do it properly
  
  # this is used to summarise how well the GAM performed
  
  mod.sum <- summary(mod)
  s2 <- mod$sig2 # estimate of standard deviation of residuals
  X <- mod$model # data used to fit the model
  n <- nrow(X) # how many observations were used in fitting?
  v <- -1 # omit the intercept term, it can't inflate variance
  varbeta <- mod.sum$p.table[v,2]^2 # variance in estimates
  varXj <- apply(X=X[,row.names(mod.sum$p.table)[v]],MARGIN=2, var) # variance of all the explanatory variables
  VIF <- varbeta/(s2/(n-1)*1/varXj) # the variance inflation factor, obtained by rearranging
  # var(beta_j) = s^2/(n-1) * 1/var(X_j) * VIF_j
  
  VIF.df <- data.frame(variable=names(VIF),
                       vif=VIF, 
                       row.names=NULL)
  
  return(VIF.df)
}
samvif(fit)

  variable      vif
1     wind 1.351802
2     temp 1.786334

計算されたのは風速と気温のみ。非線形の関数になっている日射量と日付は計算できないようです。
思ったより気温のvifは高くないですね。2以下なので問題なさそうです。

非線形の関数となっている説明変数同士はmgcvパッケージにあるconcurvity（共曲線性？）なる指標があるでようです。

concurvity

concurvity(fit)

              para  s(solar)   s(date)
worst    0.9969669 0.4109838 0.6976096
observed 0.9969669 0.2045333 0.3594092
estimate 0.9969669 0.2134965 0.3682248

worst、observed、estimateの3つの推定値が出てます。
こちらのページに英語の説明がありますがいまいちわからん。
concurvityは0から1までの値をとり、1に近づくほど多重共線性（曲線性）が強くなるようです。

どれくらいでまずくなるのかよくわかりませんが、
論文をいくつか見てみると、
0.8より下回るようにする（リンク）、
0.2より下回るようにする（リンク）、
0.5より大きければエラーが生じる（リンク）、

みたいなのがありいろいろ。

今回のケースでは日射量はworstでも0.41、日付に関してもworstでは0.7程度ですが、observedでは0.36でした。
まあ微妙なところですが少なくとも日射量に関しては問題なさそうです。

今後ももっと調べて更新します。