一般化加法モデルのVIFの算出方法をメモしておきます。
加法モデルについてはこちらを参考に。
今回はサンプルデータとして、Rに組み込みのニューヨークの大気データを使用します。
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がありますが、ここではとりあえず削除します。
列名を簡単な表記に変更しておきます。
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)
オゾンと関係ありそうなのは風速と気温でしょうか。
日射量と日付は微妙なのでとりあえずスプライン化してみます。
また気温と風速、気温と日付はそこそこ相関がありそうです。
準備ができたので早速一般化加法モデルを実行します。
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でもやってみます。
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を計算するコードが公開されてました！
リンク
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(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でした。
まあ微妙なところですが少なくとも日射量に関しては問題なさそうです。
今後ももっと調べて更新します。