「非自明」な応答のデータ

　今まで私たちは、線形回帰、非線形回帰、一般化線形モデル、と様々な応答を示すデータに関して、非線形な式やリンク関数などを通じて、うまいことフィットできるように工夫を凝らしてきた。では、こんなデータの解析はいかがだろうか？

------------------------------------------------------

library(plotn)

library(mgcv)

library(viridis)

x <- c(-12.3, 14.9, -3.8, -1.8, 6.6, 2.4, 5.8, -3.5, 8.4, -7.5, 14.1, -7, 

       -6.3, -11.2, -0.5, 8.1, 10.2, -12.3, -10.9, -11.5, 11.8, 11.6, -1.5,

       5.6, -13, 4.3, 13.8, 1.5, -8.9, 11.7, 14.3, 13.7, 10.1, 12.2, 4.8, 

       6.6, 1.9, 7.3, 10.7, -2.4, 8.7, 0, -8.5, -1.3, -10, 5.5, -1.3, -5.8,

       3.2, -7.8, 13.9, 11.1, 13.3, 6.5, -0.6, 9.1, -13.7, -6.2, 11.3, -0.5,

       9.2, 9.9, -11.5, -4.3, -2.9, -7.7, -7.4, -7.8, -11.1, -14.9, 14.4, 9.5, 

       -5.3, 0.5, -14.7, 0.8, 7.7, -11.9, 1.4, -5.2, -12.1, -4, 8.3, -5.1, 

       -13.5, -4, -12.2, 14.2, 11.1, 1.3, -9.9, -3.7, -13.2, -6.6, 8.4, 1.8,

       -8, 1.9, -1.4, 8)

y <- c(9.3, 8.4, -11.6, -1.1, -4.6, 9.7, -1, -12.4, -8.1, -0.8, 7.1, -3, -3.9,

       8.5, 6.6, -9.9, -4.4, 7.8, 7.2, 9.4, 9, 4.9, -2.5, -0.9, 9, -1, 9.7, 

       8.9, 3.4, 8.2, 7.3, 11, -6.6, 9.2, -0.1, -3.1, 7.8, -11.1, -4.9, -8.3, 

       -10.2, 8.3, 2.8, 5.3, 4.5, -2.5, 5.2, -5.2, 5.3, -1.7, 4, 6.4, 4.7, 

       -6.8, 9.8, -9.2, 3.9, -6.7, 2.9, 9.6, -10.6, -7.5, 7.7, -9.9, -9.7, 

       -0.5, -0.8, 1.9, 6.9, -5.7, 8.3, -14, -6.5, 8.1, -4.2, 10.7, -6.8, 9, 

       10.6, -5.1, 11.4, -9.5, -10.4, -5.7, 4.8, -11, 11.6, 7.2, 0, 9.7, 5.6,

       -7.9, 8.6, -6.3, -7.2, 8.1, 0.1, 9.3, 1.2, -9.1)

d <- data.frame(x = x,  y = y)

plotn(y ~ x, data = d)#図1の描画 

------------------------------------------------------

まあ、突然の出来事に私ならPCをぶん投げたくなるが（大袈裟）、世の中にはこのタイプのデータもはびこっている。明らかに周期的な挙動なので、どうせ三角関数をかませたような出力になっているのだろうが、話はシンプルではなさそうだ。よく見ると、波が上昇する時と下降する時で、波の幅が違いそうだ。上昇時は急峻、下降時は緩やか、そんな感じに見える。

　さて、非線形回帰で対処してやりたいところなのだが、果たして上記の条件を満たしてくれそうなうまい関数を設定できるだろうか？　非線形回帰を行うためには、私たちが数式を書く下してやらねばならない。思い切って「単純な三角関数で回帰してしまえ！」という気持ちにもなってくるが、どうにかうまく対処したい。

　そこで、このようなバックに潜む関数を明示的に指定しなくても、うまいこと回帰を行う方法を紹介する。それが、本稿のタイトルになっている一般化加法モデルGenelized additive model (GAM)である。

一般化加法モデルGenelized additive modelのこころ

　まずは、GAMの「こころ」について簡単に紹介する。その後、解析を行って、最後に詳細な解説をすることにしよう。まず、GAMは下記のようなモデル式を考えている。

ここでβが推定するべきパラメータである。一方、b(x)は基底関数basis functionと呼ぶ。この定義で特に基底関数の中身を定義していないが、後述するように多項式、特に3次関数をベースにした関数を基底関数とすることが多い。この基底関数をp個用意して、説明変数xを代入して変換し、その和で被説明変数yが表現できると考える。この時、「p個の基底関数をどれくらいの割合で混ぜ合わせるか」が推定されるパラメータである。このように、複数の関数の和で被説明変数を表現することから、このモデルを加法モデルadditive modelと呼ぶ。この時、誤差は正規分布に従うと考えている。さらに正規分布以外の確率分布に対してもリンク関数のさらなる変換を通じて対応できるようにしたものが一般化加法モデルGenelized additive model (GAM)である。ざっくりいえばここから先の内容は、p種類の3次関数を用意して、それを足し合わせることでぐにゃぐにゃな曲線を表現しようとする試みである。もし、基底関数が恒等関数ならそれは線形回帰やGLMに他ならない。GAMはゆえに線形回帰やGLMのさらなる拡張である。

・自然3次スプライン基底natural cubic spline basis

　今回はRのmgcvパッケージを利用することでGAMを実行するのだが、そこでは「自然3次スプライン基底の罰則項付き最小二乗法」が用いられ、この方法を指してGAMと呼ぶこともある。さて、意味不明な用語、「自然3次スプライン基底natural cubic spline basis」なるものが登場した。これがGAMを理解するうえで重要なカギになっている。

まず、「3次」の由来に関して紹介しよう。そもそも、私たちは非自明な応答を何とかして解析しようとしている。けれども、関数を指定せずして、どうやってぐにゃぐにゃな曲線にフィットさせるというのだろうか。ところで、突然だが、なにか得体のしれない関数を見たとき「より簡単な関数で近似（フィッティング）できないか」ということを考えてみたい。例えば、1次式、2次式、3次式は、非常に簡単な多項式関数だ。これらを基底関数とすることを考える。

　でも、同じ多項式関数でも、低い次数ではできないことが多い。1次式、つまり直線は曲線を表現できない。次に、2次式は曲線や極値を表現できるが、変曲点は表現できない。3次式以上なら曲線も極値も変曲点も表現が可能である。一方で、次数が多ければ多いほど、過学習の危険がある。これらのトレードオフを考慮した時に3次関数で近似することを思い至る。実際、3次関数を用いて基底関数を定義することから、これが3次スプラインにおける「3次」の由来である。

　3次関数で近似するのだとしたら、最大で極大値一つと極小値一つの関数までしか近似できないじゃないか。それはその通りである。そこでさらに次のことを考える。今推定したい非自明な応答を示すデータに関して、適当な区間で区切る。この区切りをノットknotと呼ぶ。この比較的狭い区間内だけで3次関数によるフィッティングを行うのだ。

　単純に区間で区切り、それぞれ独立に3次関数によるフィッティングを行うだけだと、区間の境界では推定された値が連続的にならない。そこで、3次関数の導関数を用意し、「3次関数が区間の境界で同じ値をとる」かつ「導関数が区間の境界で同じ値をとる」＝「共通の接線を持つよう」に制約をつけることで、全区間にわたって滑らかに曲線をつなぐことができる。この滑らかな曲線をスプラインsplineと呼び、3次関数を使ってスプラインを構築する方法を3次スプラインcubic splineと呼ぶ。

　こうやって、推定していくと、データの両端の区間はその先にデータがなく、「3次関数が区間の境界で同じ値をとる」かつ「導関数が区間の境界で同じ値をとる」制約がなくなってしまう。このとき、データのない領域を推定しよう（これを外挿とよぶ）とすると、3次関数にフィットしているわけであるから、データの予測値は非常に大きく振れることになる。このような大振れを避けるため、最後の制約として、この境界の両端において3次関数の二次導関数が0となるような制約をつけて推定する。二次導関数 = 0ということは、データの両端点は変曲点になるということであり、つまりその瞬間は上に凸でも下に凸でもなく直線的変動になる制約を設けたことになる。こうすれば、極端な挙動が避けられる。この制約を課したスプラインを自然スプラインnatural splineと呼び、以上をすべて満たしたスプラインを自然3次スプラインと呼ぶわけだ。この自然3次スプラインは平滑化スプラインsmoothing splineの1種である。

　最もわかりやすい例を示すことで、自然3次スプライン基底による回帰の有用性を示そう。今、ノットがK個あり、ノットの位置のxの値をa1 < a2 <  …… < aK-1 < aKとする。a1 = min(x)、aK = max(x)である。m（m = 1, 2, ……, K-2）に対応するノットをamとし、このam対応する自然3次スプライン基底をb2+m(x)とすると、b2+m(x)の最もわかりやすいものの一つ（唯一ではないことに注意）は、

　　x < amでは0、

　　am ≦ x < aK-1では3次関数、

　　aK-1 ≦ x < aKでは別の3次関数、

 　　aK < xでは1次関数

となる滑らかな連続関数である（具体的な数式は後述）。例えば、K = 4のとき、下記のような見た目である。

このように定義される自然3次スプライン基底が、K-2個存在することになる。なお、b1(x) = 1、b2(x) = xとしておく。これで、b1(x) ~ bK(x)のK個の基底関数を作ることができた。このK個の基底関数の混ぜ合わせでK個のノットで区切られたデータを推定する。

　連続関数の和は連続関数になることから、K個の自然3次スプライン基底の混ぜ合わせは、a1以下の領域は1次関数、ノット間で区切られた領域はそれぞれが異なる3次関数、aK以上の領域は1次関数として、それぞれが滑らかにつながるように推定することと等しくなる。例えば、K = 4のときは、下記のようなことを想定していることになる。確かに、区間を分けて3次式による回帰を行いつつも、連続で、かつデータの定義域外では1次関数となり、満たしてほしい性質を兼ね備えていることがわかる。

上記では最もわかりやすい自然3次スプライン基底を示したが、mgcvパッケージのgam関数では別の自然3次スプライン基底を使っているので注意してほしい。

　最後、ここで扱うGAMは「自然3次スプライン基底の罰則項付き最小二乗法」なわけだが、何を罰則項とするのかを紹介しよう。天下り的だが、当てはめた自然3次スプラインの二次導関数は元のスプラインの曲がり具合を表す関数である。二次導関数が常に0である関数といえば定数関数や1次関数のような直線で、これらの関数は全くまがっていない。「高次の関数」や「同じ次数でも最高次にかかってる係数が大きい」ほど二次導関数の変数にかかる係数は大きいから、この値の絶対値が大きいほど、より極端にぐにゃぐにゃしていることを示している。よりぐにゃぐにゃしているほど、柔軟に当てはまる一方、過学習の原因となる。そこで、ぐにゃぐにゃ具合を制御するために、二次導関数の2乗のデータ全体での積分を罰則項とすることで、ぐにゃぐにゃしすぎることを防ぐ。

1変数のGAM

　では、上記のデータを解析してみよう。以下のデータを改めて提示する。

------------------------------------------------------

x <- c(-12.3, 14.9, -3.8, -1.8, 6.6, 2.4, 5.8, -3.5, 8.4, -7.5, 14.1, -7, 

       -6.3, -11.2, -0.5, 8.1, 10.2, -12.3, -10.9, -11.5, 11.8, 11.6, -1.5,

       5.6, -13, 4.3, 13.8, 1.5, -8.9, 11.7, 14.3, 13.7, 10.1, 12.2, 4.8, 

       6.6, 1.9, 7.3, 10.7, -2.4, 8.7, 0, -8.5, -1.3, -10, 5.5, -1.3, -5.8,

       3.2, -7.8, 13.9, 11.1, 13.3, 6.5, -0.6, 9.1, -13.7, -6.2, 11.3, -0.5,

       9.2, 9.9, -11.5, -4.3, -2.9, -7.7, -7.4, -7.8, -11.1, -14.9, 14.4, 9.5, 

       -5.3, 0.5, -14.7, 0.8, 7.7, -11.9, 1.4, -5.2, -12.1, -4, 8.3, -5.1, 

       -13.5, -4, -12.2, 14.2, 11.1, 1.3, -9.9, -3.7, -13.2, -6.6, 8.4, 1.8,

       -8, 1.9, -1.4, 8)

y <- c(9.3, 8.4, -11.6, -1.1, -4.6, 9.7, -1, -12.4, -8.1, -0.8, 7.1, -3, -3.9,

       8.5, 6.6, -9.9, -4.4, 7.8, 7.2, 9.4, 9, 4.9, -2.5, -0.9, 9, -1, 9.7, 

       8.9, 3.4, 8.2, 7.3, 11, -6.6, 9.2, -0.1, -3.1, 7.8, -11.1, -4.9, -8.3, 

       -10.2, 8.3, 2.8, 5.3, 4.5, -2.5, 5.2, -5.2, 5.3, -1.7, 4, 6.4, 4.7, 

       -6.8, 9.8, -9.2, 3.9, -6.7, 2.9, 9.6, -10.6, -7.5, 7.7, -9.9, -9.7, 

       -0.5, -0.8, 1.9, 6.9, -5.7, 8.3, -14, -6.5, 8.1, -4.2, 10.7, -6.8, 9, 

       10.6, -5.1, 11.4, -9.5, -10.4, -5.7, 4.8, -11, 11.6, 7.2, 0, 9.7, 5.6,

       -7.9, 8.6, -6.3, -7.2, 8.1, 0.1, 9.3, 1.2, -9.1)

d <- data.frame(x = x,  y = y)

------------------------------------------------------

　mgcvパッケージのgam関数を使うことでGAMを使うことができる。引数の指定の仕方はlmやglmと非常に類似しているが、説明変数を式中でsmooth termとして指定するs関数に噛ませて、入力する。s関数内のbs引数は、平滑化基底smoothing basisを指定、crは3次スプラインである。

------------------------------------------------------

fit1 <- gam(y ~ s(x, bs = "cr"), data = d)

fit1

## 

## Family: gaussian 

## Link function: identity 

## 

## Formula:

## y ~ s(x, bs = "cr")

## 

## Estimated degrees of freedom:

## 8.84  total = 9.84 

## 

## GCV score: 7.185584

fit2 <- summary(fit1)

fit2

## 

## Family: gaussian 

## Link function: identity 

## 

## Formula:

## y ~ s(x, bs = "cr")

## 

## Parametric coefficients:

##             Estimate Std. Error t value Pr(>|t|)   

## (Intercept)   0.8590     0.2545   3.375  0.00109 **

## ---

## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 

## Approximate significance of smooth terms:

##        edf Ref.df     F p-value    

## s(x) 8.843  8.992 82.59  <2e-16 ***

## ---

## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 

## R-sq.(adj) =  0.882   Deviance explained = 89.3%

## GCV = 7.1856  Scale est. = 6.4783    n = 100

------------------------------------------------------

大体、summary関数に噛ませた後に出力されるものに共通しているので、fit2の中身を確認してゆこう。

　Parametric coefficientでは、smooth termではないパラメータの推定値が得られている。今回、明示的に表記していないが、いつものごとく、切片項が自動で推定される。切片推定値は有意に0とは異なることがp値からわかる。

　Approximate significance of smooth termsでは、smooth termが有意かどうかの解析がなされている。通常の線形回帰などのように、ここに具体的な係数が表示されるわけでもなく、何の値が有意に0から異なるのかが判然としない。推定された係数はcoef関数を使って確認できる。

------------------------------------------------------

coef(fit1)

## (Intercept)      s(x).1      s(x).2      s(x).3      s(x).4      s(x).5 

##   0.8590000   9.7540274  -0.3050497  -9.1540824   1.4268564   9.8898512 

##      s(x).6      s(x).7      s(x).8      s(x).9 

##  -5.6685070  -8.8493424   6.0360021   6.4016260

------------------------------------------------------

Intercept項は良いとして、残りのs(x).1 ~ 9が説明変数xに対する基底関数にかかる推定パラメータである。9種類（切片項も含めて考えてもよいので10種類）の自然3次スプライン基底が作られ、それにかかるパラメータが0ではない可能性が高いならば、有意であるとしてApproximate significance of smooth termsで検出される。

　edfは有効自由度effective degree of freedomで、そもそも9種類の基底関数は互いに独立ではないので推定される切片以外のパラメータも独立ではない。独立ではないことを加味した自由度がedfで8.84程度となっている。通常の線形回帰なら、推定されたパラメータは9個だから自由度が9になることを思い出そう。切片は独立なので切片の自由度は1である。ゆえに、fit1の中身のEstimated degrees of freedomではこれらの和の9.84が登場する。ref.dfは参照自由度と呼ばれるものだが、mgcvパッケージの作者はこの値にあまり意味はないと表現している。

　GCVは一般化交差検証値generalized cross validationと呼ばれる重要な値で、この値が最小になるようにβが推定されている。この値は、罰則項の時と類似の発想で、係数推定に使っていないデータに対して推定されたモデルを当てはめたときの当てはまりの悪さを示すものである。この値が小さいということは、未知のデータに対しても当てはまりが良いことを示す、といってもよいだろう。

　GAMのresult内にはほかにも重要な情報が入っている。例えば、下記にはノットをどこにしたかの情報が入っている。

------------------------------------------------------

fit1$smooth[[1]]$xp

##         0%  11.11111%  22.22222%  33.33333%  44.44444%  55.55556%  66.66667% 

## -14.900000 -11.544444  -7.722222  -5.133333  -1.444444   1.633333   6.533333 

##  77.77778%  88.88889%       100% 

##   9.266667  11.844444  14.900000

------------------------------------------------------

3次スプラインではノットの個数10になる（mgcvパッケージ、smooth.construct.cr.smooth.specのhelpより。）。ノットの位置の最小値はxの最小値、最大値はxの最大値である。K個のノットを定めるとすると、ノットの位置はxの値に関して、0/K × 100、1/K × 100……、K/K × 100パーセンタイルの点になっている。つまり、各ノット間に挟まれるxの個数が均等になるように設定されるのがデフォルトである。

　では回帰曲線を描いてみることにする。今まで通りの方法なら自然3次スプライン基底をもとに回帰曲線を描くが、やり方は後程紹介するが、煩雑なのでpredict関数と使うのが便利だろう。今回は区間推定も併せて行ってしまう。誤差構造は正規分布を仮定しているから、線形回帰の信頼区間と予測区間の推定を使うことができる。予測区間推定にはfit1$sigにデータの分散の不偏推定が記録されているのでそれを使おう。

------------------------------------------------------

fit1$sig2

## [1] 6.478292

sum(fit1$residuals^2)/df.residual(fit1)

## [1] 6.478292

xx <- seq(min(d$x), max(d$x), length = 200)

newdata <- data.frame(x = xx)

pred <- predict(fit1, newdata = newdata, se.fit = T)

c_upr <- pred$fit + qt(p = 0.975, df = df.residual(fit1)) * pred$se.fit

c_lwr <- pred$fit - qt(p = 0.975, df = df.residual(fit1)) * pred$se.fit

p_upr <- pred$fit + qt(p = 0.975, df = df.residual(fit1)) * 

  sqrt(pred$se.fit^2 + fit1$sig2)

p_lwr <- pred$fit - qt(p = 0.975, df = df.residual(fit1)) * 

  sqrt(pred$se.fit^2 + fit1$sig2)

lim <- range(c(p_upr, p_lwr, d$y))

plotn(y ~ x, data = d, ylim = lim)#図2の描画 

overdraw(polygon(c(xx, rev(xx)), c(c_upr, rev(c_lwr)), 

                 border = F, col = "#00000030"),

         polygon(c(xx, rev(xx)), c(p_upr, rev(p_lwr)), 

                 border = F, col = "#00000030"),

         points(xx, pred$fit, type = "l", col = "red"))

------------------------------------------------------

一応、ノットの位置も示しておこう。

------------------------------------------------------

plotn(y ~ x, data = d, ylim = lim)#図3の描画 

overdraw(polygon(c(xx, rev(xx)), c(c_upr, rev(c_lwr)), 

                 border = F, col = "#00000030"),

         polygon(c(xx, rev(xx)), c(p_upr, rev(p_lwr)), 

                 border = F, col = "#00000030"),

         points(xx, pred$fit, type = "l", col = "red"))

for(i in 1:length(fit1$smooth[[1]]$xp)){

  overdraw(abline(v = fit1$smooth[[1]]$xp[i]))

}

------------------------------------------------------

では、もし、xの定義域外まで予測を伸ばしたらどうなるだろうか。

------------------------------------------------------

xx <- seq(-20, 20, length = 200)

newdata <- data.frame(x = xx)

pred <- predict(fit1, newdata = newdata, se.fit = T)

c_upr <- pred$fit + qt(p = 0.975, df = df.residual(fit1)) * pred$se.fit

c_lwr <- pred$fit - qt(p = 0.975, df = df.residual(fit1)) * pred$se.fit

p_upr <- pred$fit + qt(p = 0.975, df = df.residual(fit1)) * 

  sqrt(pred$se.fit^2 + fit1$sig2)

p_lwr <- pred$fit - qt(p = 0.975, df = df.residual(fit1)) * 

  sqrt(pred$se.fit^2 + fit1$sig2)

lim <- range(c(p_upr, p_lwr, d$y))

plotn(y ~ x, data = d, xlim = range(xx), ylim = lim)#図4の描画 

overdraw(polygon(c(xx, rev(xx)), c(c_upr, rev(c_lwr)), 

                 border = F, col = "#00000030"),

         polygon(c(xx, rev(xx)), c(p_upr, rev(p_lwr)), 

                 border = F, col = "#00000030"),

         points(xx, pred$fit, type = "l", col = "red"))

------------------------------------------------------

確かに、定義域外ではxは直線的な挙動の推定となっている。上記のデータに関して、明らかに周期的であり、おそらく定義域外でも増減を繰り返すだろう。予測を見るに、本解析のGAMでは周期性を考慮しているわけではないことがわかる。

　では、最後、基底関数の「挙動」を取り出して確認し、それを使って回帰曲線を描く。model.matrix関数にgamオブジェクトを入力するとモデル内でのxに対する「重みづけ」を取り出せるのだが、これが基底関数に対応している。

------------------------------------------------------

m_mat1 <- model.matrix(fit1)

for(i in 1:ncol(m_mat1)){

  tmp <- data.frame(x = d$x, basis.f = m_mat1[,i])

  tmp <- tmp[order(tmp$x),]

  plotn(basis.f ~ x, data = tmp, type = "l")#図5～14の描画 

}

------------------------------------------------------

以上のように、初めに紹介した自然3次スプライン基底と比べて、GAM内で動いている基底関数はかなり複雑な挙動である。だがこれらの関数も同様に、足し合わせるとノット間で3次関数に、ノット外で1次関数になるように働く。

　これらの基底関数にgam関数で推定した係数を掛けて併せて足しあげることで回帰曲線が出来上がる。

------------------------------------------------------

res <- sweep(m_mat1, 2, coef(fit1),`*`)

tmp <- data.frame(x = d$x, res = rowSums(res))

tmp <- tmp[order(tmp$x), ]

plotn(y ~ x, data = d)#図15の描画 

overdraw(points(tmp$x, tmp$res, type = "l", col = "red"))

------------------------------------------------------

もし、バックで動いている基底関数の計算方法がわかれば、上記のようなカクカクの回帰曲線ではなく滑らかな回帰曲線を描けるが、私の力不足でまだわかっていない。いつか計算できたら、ここに追記しようと思う。以上のようにGAMによって、こちらが関数を定義しなくても非線形回帰を行うことができた。

　ちなみに今回は、以下のような関数を定義して、これに正規分布に従う誤差をつけて生成した。この関数は増加する時と減少する時で、全域で微分可能だが異なる三角関数から値を生成する。f1とf2で三角関数の周波数、aで振幅、thetaで位相を定義する。

------------------------------------------------------

sstw <- function(x, f1 = 1, f2 = 1, a = 1, theta = 0){

  res <- sapply(x, FUN = function(org) {

    tmp <- org - theta

    if(tmp <= 0){

      while(tmp <= 0){

        tmp <- tmp + (1/f1 + 1/f2) * pi

        }

    } else {

      while(tmp > (1/f1 + 1/f2) * pi){

        tmp <- tmp - (1/f1 + 1/f2) * pi

        }

    }

    if(0<= tmp & tmp < pi/f1){

      a*cos(f1*tmp)

    } else {

      a*cos(f2*tmp - (f2/f1 - 1)*pi)

    }

  }

  )

  return(unlist(res))

}

xx <- seq(0, 20, length = 2000)

plotn(xx, sstw(xx, f1 = 1/3, f2 = 1), type = "l")#図16の描画 

------------------------------------------------------

例えば、以下のように生成したデータも解析できる。

------------------------------------------------------

n <- 100

s <- 2

x <- round(runif(n, -15, 15), digits = 1)

y <- rnorm(n = n, mean = sstw(x, f1 = 2/3, f2 = 3/4, a = -10), sd = s) + x

x <- c(-4.2, -4.4, 11.5, -6.1, 2.3, 10.4, -0.4, -0.7, -5.8, -1.4, -2.4, 11.6, 9.9, 13.5, 2.1, -6.5, -1.2, 3.7, 12, 6, 9.6, -7.6, 0.6, -10.9, 3.4, -3.4, 13.7, -5.8, 7.1, 14.9, 0.4, 7.4, -9.5, -7.3, -7.5, -5.4, -9.3, 8.1, 3.7, -8.9, 10.5, 0, -11, -10.1, -7.4, 4.7, -0.6, 1.6, -1.3, 3.9, 14.9, 12.5, -8.1, -4.4, -8.9, -6, -9.2, -13, -5.4, 10.4, 10.1, 8.9, -4.2, 14.3, 5.9, -9.9, -6.8, -5.7, -9, -2, -1, -7.5, 14.3, 10.7, 12.9, -10.2, -0.5, -9, -1.9, 1.6, 1.8, -10.8, -6.7, -2.7, 4, -11.1, 6.7, -14.2, 1.2, 2.8, -10.3, 1.7, -0.4, 7.9, 1.3, 0.6, 4.5, -4.1, -4.6, -12.3)

y <- c(4.9, 6.9, 12.7, -4.1, 2.2, 7.7, -8.9, -9, 0, -4.9, 3.1, 16.6, 0.2, 24.5, 0.1, -6.5, -9.3, 12.6, 15.6, 10.8, -1.7, -15.3, -10.7, -12.3, 9.2, 2, 21, -2.4, 5.1, 19.5, -8, 0.8, -18.3, -14.7, -14.1, 4.8, -23.9, -0.6, 10.6, -21.3, 4.9, -6.9, -9, -13.5, -13.9, 12.8, -12.6, -2.4, -9.8, 12.6, 21.1, 20.6, -13.9, 6.3, -17.7, -0.5, -18.9, -2.6, 0.3, 5.8, 3.7, 2.9, 4.5, 23.6, 14.1, -16.6, -5.8, -1.3, -15.8, -1, -7.3, -15.2, 24, 6.4, 19.9, -15.4, -10.2, -22.5, -1.4, -5.3, 0.3, -15.2, -9.1, 4, 14.3, -10.9, 5.7, -9.9, -6.5, 7.3, -15.1, 0.7, -14.3, 0.7, -6.1, -8.9, 11.1, 6.7, 2.3, -5.6)

d <- data.frame(x = x,  y = y)

plotn(y ~ x, data = d)#図17の描画 

------------------------------------------------------

解析は以下の通りで、予測区間などの図示を行う。

------------------------------------------------------

fit1 <- gam(y ~ s(x, bs = "cr"), data = d)

fit2 <- summary(fit1)

fit2

## 

## Family: gaussian 

## Link function: identity 

## 

## Formula:

## y ~ s(x, bs = "cr")

## 

## Parametric coefficients:

##             Estimate Std. Error t value Pr(>|t|)    

## (Intercept)  -0.9960     0.2792  -3.568 0.000579 ***

## ---

## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 

## Approximate significance of smooth terms:

##        edf Ref.df   F p-value    

## s(x) 8.925  8.998 181  <2e-16 ***

## ---

## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 

## R-sq.(adj) =  0.943   Deviance explained = 94.8%

## GCV = 8.6519  Scale est. = 7.7932    n = 100

xx <- seq(min(d$x), max(d$x), length = 200)

newdata <- data.frame(x = xx)

pred <- predict(fit1, newdata = newdata, se.fit = T)

c_upr <- pred$fit + qt(p = 0.975, df = df.residual(fit1)) * pred$se.fit

c_lwr <- pred$fit - qt(p = 0.975, df = df.residual(fit1)) * pred$se.fit

p_upr <- pred$fit + qt(p = 0.975, df = df.residual(fit1)) * 

  sqrt(pred$se.fit^2 + fit1$sig2)

p_lwr <- pred$fit - qt(p = 0.975, df = df.residual(fit1)) * 

  sqrt(pred$se.fit^2 + fit1$sig2)

lim <- range(c(p_upr, p_lwr, d$y))

plotn(y ~ x, data = d, ylim = lim)#図18の描画 

overdraw(polygon(c(xx, rev(xx)), c(c_upr, rev(c_lwr)), 

                 border = F, col = "#00000030"),

         polygon(c(xx, rev(xx)), c(p_upr, rev(p_lwr)), 

                 border = F, col = "#00000030"),

         points(xx, pred$fit, type = "l", col = "red"))

------------------------------------------------------

では、こんなデータはいかがだろうか。

------------------------------------------------------

n <- 100

s <- 2

x <- round(runif(n, -15, 15), digits = 1)

y <- rnorm(n = n, mean = sstw(x, f1 = 2, f2 = 1, a = 0.5), sd = s) + x

x <- c(-4.7, -6.3, -4, 10.2, -2.7, -10.8, -0.3, 14.7, -10.6, -14, -3.3, 11.8, 9.1, -5, -2.2, -8, -11.1, -12.2, -6.6, 1, -3.4, 0.2, -9.8, 9.5, 2.5, -5, -0.4, -7.4, 5.4, -10.9, -14.9, 4.6, -11.1, -2.2, -4.7, 8.9, -13.6, 9.8, 7.2, 5, -0.5, -9.7, -3.7, -2, 2, 11.8, 12.8, 2.7, 5.5, -10.9, -2.3, 11.5, 14.4, -12.2, 12.5, 1.3, 1.1, 11.6, -4.1, -2.9, 14.4, -0.8, -8.7, -8.3, -9.2, -8.6, -8, 0.8, 0.5, 11.1, -13.3, 5.7, -4.1, -13.7, -7.4, 0.7, 13, 4.6, -3.2, 5.9, 1.1, 7.3, -10, 13.1, -7.1, 4.6, -13.8, -10.2, 10.5, -9.3, -5.2, -5.2, 8.8, 7.6, -13, -13.5, 11.6, -10.3, -6.3, -8.1)

y <- c(0.4, 0.1, 0.1, 0, -0.5, 0.1, 0.5, 0.2, 0.1, 0.6, -0.5, -0.3, 0.5, 0.5, -0.3, -0.6, -0.2, -0.5, -0.1, -0.2, -0.4, 0.4, 0.4, 0.6, -0.4, 0.4, 0.6, -0.3, 0.2, 0.1, 0.3, 0.5, -0.1, -0.2, 0.4, 0.3, 0.3, 0.5, -0.5, 0.4, 0.5, 0.6, -0.3, -0.3, -0.3, -0.4, 0.1, -0.2, 0.1, 0, -0.4, -0.3, 0.4, -0.4, -0.1, -0.7, -0.3, -0.4, 0.1, -0.5, 0.3, 0.3, 0.1, -0.3, 0.6, -0.1, -0.4, 0, 0.3, -0.3, -0.2, -0.2, 0.3, 0.2, -0.2, 0.2, 0.1, 0.5, -0.4, -0.4, -0.4, -0.3, 0.6, 0.5, -0.4, 0.3, 0.3, 0.3, -0.3, 0.5, 0.5, 0.3, 0.4, -0.2, -0.4, 0.1, -0.5, 0.3, 0, -0.5)

d <- data.frame(x = x,  y = y)

plotn(y ~ x, data = d)#図19の描画 

------------------------------------------------------

データ生成過程はcosカーブから生成されているのだが、全体的に散らばっていて、うがった目で見れば波……？かな、という感じである。デフォルトの解析結果では以下の通り。

------------------------------------------------------

fit1 <- gam(y ~ s(x, bs = "cr"), data = d)

fit2 <- summary(fit1)

fit2

## 

## Family: gaussian 

## Link function: identity 

## 

## Formula:

## y ~ s(x, bs = "cr")

## 

## Parametric coefficients:

##             Estimate Std. Error t value Pr(>|t|)

## (Intercept)  0.02100    0.03668   0.573    0.568

## 

## Approximate significance of smooth terms:

##        edf Ref.df    F p-value

## s(x) 1.039  1.077 0.13   0.797

## 

## R-sq.(adj) =  -0.00863   Deviance explained = 0.195%

## GCV = 0.13734  Scale est. = 0.13454   n = 100

------------------------------------------------------

デフォルトでは、smooth term項は有意でないので、非線形な効果は小さい可能性が高い。図示して見ると、

------------------------------------------------------

xx <- seq(min(d$x), max(d$x), length = 200)

newdata <- data.frame(x = xx)

pred <- predict(fit1, newdata = newdata, se.fit = T)

c_upr <- pred$fit + qt(p = 0.975, df = df.residual(fit1)) * pred$se.fit

c_lwr <- pred$fit - qt(p = 0.975, df = df.residual(fit1)) * pred$se.fit

p_upr <- pred$fit + qt(p = 0.975, df = df.residual(fit1)) * 

  sqrt(pred$se.fit^2 + fit1$sig2)

p_lwr <- pred$fit - qt(p = 0.975, df = df.residual(fit1)) * 

  sqrt(pred$se.fit^2 + fit1$sig2)

lim <- range(c(p_upr, p_lwr, d$y))

plotn(y ~ x, data = d, ylim = lim)#図20の描画 

overdraw(polygon(c(xx, rev(xx)), c(c_upr, rev(c_lwr)), 

                 border = F, col = "#00000030"),

         polygon(c(xx, rev(xx)), c(p_upr, rev(p_lwr)), 

                 border = F, col = "#00000030"),

         points(xx, pred$fit, type = "l", col = "red"))

------------------------------------------------------

ほとんど線形回帰と同じ結果である。

　上記のデータのように細かく振動する割にデータが疎だとsmooth termが有意にならない。なぜかと言えば、デフォルトだとノットの数は大体10で、ノット間が3次関数で近似できる形にならないからだ（例えば、ノット間に含まれるデータが山が2つあるみたいになる）。そこでノットを十分に増やしてやると、上記のデータを非線形的に捉えることも可能になる。s関数にはいくつかの引数があり、k引数がノットの数を制御する。例えば、下記のようにノットを20個にできる。

------------------------------------------------------

fit1 <- gam(y ~ s(x, bs = "cr", k = 20), data = d)

fit2 <- summary(fit1)

fit2

## 

## Family: gaussian 

## Link function: identity 

## 

## Formula:

## y ~ s(x, bs = "cr", k = 20)

## 

## Parametric coefficients:

##             Estimate Std. Error t value Pr(>|t|)

## (Intercept)  0.02100    0.01365   1.538    0.128

## 

## Approximate significance of smooth terms:

##        edf Ref.df     F p-value    

## s(x) 18.58  18.96 32.51  <2e-16 ***

## ---

## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 

## R-sq.(adj) =   0.86   Deviance explained = 88.7%

## GCV = 0.023175  Scale est. = 0.018636  n = 100

xx <- seq(min(d$x), max(d$x), length = 200)

newdata <- data.frame(x = xx)

pred <- predict(fit1, newdata = newdata, se.fit = T)

c_upr <- pred$fit + qt(p = 0.975, df = df.residual(fit1)) * pred$se.fit

c_lwr <- pred$fit - qt(p = 0.975, df = df.residual(fit1)) * pred$se.fit

p_upr <- pred$fit + qt(p = 0.975, df = df.residual(fit1)) * 

  sqrt(pred$se.fit^2 + fit1$sig2)

p_lwr <- pred$fit - qt(p = 0.975, df = df.residual(fit1)) * 

  sqrt(pred$se.fit^2 + fit1$sig2)

plotn(y ~ x, data = d, ylim = lim)#図21の描画 

overdraw(polygon(c(xx, rev(xx)), c(c_upr, rev(c_lwr)), 

                 border = F, col = "#00000030"),

         polygon(c(xx, rev(xx)), c(p_upr, rev(p_lwr)), 

                 border = F, col = "#00000030"),

         points(xx, pred$fit, type = "l", col = "red"))

------------------------------------------------------

すると、smooth termは有意になり、非線形な効果を検出することができた。ノットの数が少なすぎると未学習になるし、ノットの数が増えるほど過学習が生じやすくなる。一方で、ノットが少なければノット間のデータ数が増えるから回帰の信頼度が上がるし、多ければ信頼度が下がる。じゃあ、ノットの数はどう決めてやればいいのだろうか。そこで、下記のようにノットの数を変えながら、そのときのsmooth termの有効自由度、GCV、AIC、そしてそのときにsmooth termが有意だったかどうかを記録する。

------------------------------------------------------

edfs <- NULL

gcvs <- NULL

aics <- NULL

pchs <- NULL

alpha <- 0.05

xx <- c(3, 5, 10, 20, 30, 40, 50, 60, 70, 80)

for(i in xx){

  fit1 <- gam(y ~ s(x, bs = "cr", k = i), data = d)

  fit2 <- summary(fit1)

  edfs <- c(edfs, fit2$edf)

  gcvs <- c(gcvs, fit1$gcv.ubre)

  aics <- c(aics, fit1$aic)

  if(fit2$s.pv > alpha){

    tmp <- 4#有意でないなら×シンボル

  } else {

    tmp <- 16有意でないなら●シンボル

  }

  pchs <- c(pchs, tmp)

}

------------------------------------------------------

それぞれ図示してみる。

------------------------------------------------------

plotn(edfs ~ xx, pch = pchs)#図22の描画 

plotn(gcvs ~ xx, pch = pchs)#図23の描画 

plotn(aics ~ xx, pch = pchs)#図24の描画 

------------------------------------------------------

有効自由度、GCV、AICはノットの数に対して線形に変化するのではなく、途中で頭打ちになることがわかる。ということは、ノットの数を増やしすぎても、途中から当てはまりに大した改善が見られなくなることを意味する。どこを基準に考えるかは難しいが、いくつかの解析の経験としてGCV、AICのいずれかがサチったときの最小のノット数を選ぶだろうか。であれば、k = 30くらいを選ぶ。

------------------------------------------------------

fit1 <- gam(y ~ s(x, bs = "cr", k = 30), data = d)

fit2 <- summary(fit1)

fit2

## 

## Family: gaussian 

## Link function: identity 

## 

## Formula:

## y ~ s(x, bs = "cr", k = 30)

## 

## Parametric coefficients:

##             Estimate Std. Error t value Pr(>|t|)  

## (Intercept)  0.02100    0.01063   1.975    0.052 .

## ---

## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 

## Approximate significance of smooth terms:

##        edf Ref.df     F p-value    

## s(x) 26.27   28.2 38.14  <2e-16 ***

## ---

## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 

## R-sq.(adj) =  0.915   Deviance explained = 93.8%

## GCV = 0.015543  Scale est. = 0.011304  n = 100

xx <- seq(min(d$x), max(d$x), length = 200)

newdata <- data.frame(x = xx)

pred <- predict(fit1, newdata = newdata, se.fit = T)

c_upr <- pred$fit + qt(p = 0.975, df = df.residual(fit1)) * pred$se.fit

c_lwr <- pred$fit - qt(p = 0.975, df = df.residual(fit1)) * pred$se.fit

p_upr <- pred$fit + qt(p = 0.975, df = df.residual(fit1)) * 

  sqrt(pred$se.fit^2 + fit1$sig2)

p_lwr <- pred$fit - qt(p = 0.975, df = df.residual(fit1)) * 

  sqrt(pred$se.fit^2 + fit1$sig2)

lim <- range(c(p_upr, p_lwr, d$y))

plotn(y ~ x, data = d, ylim = lim)#図25の描画 

overdraw(polygon(c(xx, rev(xx)), c(c_upr, rev(c_lwr)), 

                 border = F, col = "#00000030"),

         polygon(c(xx, rev(xx)), c(p_upr, rev(p_lwr)), 

                 border = F, col = "#00000030"),

         points(xx, pred$fit, type = "l", col = "red"))

------------------------------------------------------

こんな風に自由にノットを変えられるなら、本当はそんな関係がないのに恣意的なノットの指定でそう見せかけられたりしないだろうか。本当にそのような効果がないデータを試しに解析してみよう。下記のように単純な正規ノイズに従うデータを回帰してみよう。

------------------------------------------------------

n <- 100

s <- 2

x <- round(runif(n, -15, 15), digits = 1)

y <- rnorm(n = n, mean = 0, sd = s)

x <- c(-6, -8.9, 4.2, -3.1, 0, 14.2, 5.5, -4.2, 7, 0.8, -6, -11, -10.3, 6.4, -13.2, 13.4, -5.3, 8.2, -12, 1.7, -11.1, -12.6, -1.4, -6.7, -6, 4.2, 14.2, 3.1, -4.8, -11.7, -2.7, 14.4, -0.7, 5, 10.5, 9, -4.3, -14.5, -4.8, -10.9, 0.3, -13.9, -8, -5.8, -4.4, -0.1, 6.7, 7.1, 13.6, 6.5, -1.9, -12.8, -12.5, -14.8, -4.5, 3.7, 1, -6.5, 9.6, 1.6, -4.6, 2.2, -8.7, 8.1, -1.2, 8.3, -3.1, 0.8, -8.5, 2.9, 3.8, -7.9, 10.6, -4.3, 3.7, 12, 0.2, 13.9, -9, -3.9, -14.1, 13.5, 6.3, 5.7, -6.9, 10.2, -5.8, -1.8, -11.6, -0.6, 8.5, -5.8, 5.9, -10.8, 2.1, 0.5, 1.7, -13.9, 4.3, 4.7)

y <- c(0.4, -2.4, 1.5, -0.4, -0.1, -0.5, 3.2, 0.8, -1.9, -2.3, -0.3, -2.7, 2.5, -0.6, 1.4, -3.2, -1.4, -0.1, -2.9, 1.4, 0.2, -0.6, 2.1, 1.1, 3.8, 1.4, 2.6, 4.4, 2.2, -0.4, -1.8, -1, -0.1, 0.8, 3, 2.9, -0.8, 2.5, 1.4, 0.9, -0.4, 1.9, -2.8, -2.1, -3, 0.3, -0.9, -1.9, 0.6, 0.4, 3.6, 2.8, -1, -0.4, 1.3, 0.6, -2.1, 1.4, -1, 3.1, 1.4, 0.4, -0.3, -0.5, 1.5, 0.7, -1, 1.3, -3.2, 1.2, 5.7, -1, -3.1, -1.9, 0, 0.5, 1, 0.8, 0.6, 1.7, -0.1, 5.2, 5.2, -0.9, 0.5, -2.3, -2.4, -0.6, -4.8, 3.4, -3.1, 1.6, 0.1, -1.3, -2.8, 1.9, 1.4, 3.7, -2.8, 1.2)

d <- data.frame(x = x,  y = y)

plotn(y ~ x, data = d)#図26の描画 

------------------------------------------------------

これをGAMで解析する。

------------------------------------------------------

fit1 <- gam(y ~ s(x, bs = "cr"), data = d)

fit2 <- summary(fit1)

fit2

## 

## Family: gaussian 

## Link function: identity 

## 

## Formula:

## y ~ s(x, bs = "cr")

## 

## Parametric coefficients:

##             Estimate Std. Error t value Pr(>|t|)

## (Intercept)   0.2630     0.2074   1.268    0.208

## 

## Approximate significance of smooth terms:

##        edf Ref.df     F p-value

## s(x) 4.202  5.147 0.955   0.463

## 

## R-sq.(adj) =  0.0303   Deviance explained = 7.15%

## GCV = 4.5394  Scale est. = 4.3032    n = 100

xx <- seq(min(d$x), max(d$x), length = 200)

newdata <- data.frame(x = xx)

pred <- predict(fit1, newdata = newdata, se.fit = T)

c_upr <- pred$fit + qt(p = 0.975, df = df.residual(fit1)) * pred$se.fit

c_lwr <- pred$fit - qt(p = 0.975, df = df.residual(fit1)) * pred$se.fit

p_upr <- pred$fit + qt(p = 0.975, df = df.residual(fit1)) * 

  sqrt(pred$se.fit^2 + fit1$sig2)

p_lwr <- pred$fit - qt(p = 0.975, df = df.residual(fit1)) * 

  sqrt(pred$se.fit^2 + fit1$sig2)

lim <- range(c(p_upr, p_lwr, d$y))

plotn(y ~ x, data = d, ylim = lim)#図27の描画 

overdraw(polygon(c(xx, rev(xx)), c(c_upr, rev(c_lwr)), 

                 border = F, col = "#00000030"),

         polygon(c(xx, rev(xx)), c(p_upr, rev(p_lwr)), 

                 border = F, col = "#00000030"),

         points(xx, pred$fit, type = "l", col = "red"))

------------------------------------------------------

回帰曲線は無理やりあてはめられているが、smooth termは有意にはならない。では、上記と同様にノットの数を変えながら、そのときのsmooth termの有効自由度、GCV、AIC、そしてそのときにsmooth termが有意だったかどうかを記録する。

------------------------------------------------------

edfs <- NULL

gcvs <- NULL

aics <- NULL

pchs <- NULL

alpha <- 0.05

xx <- c(3, 5, 10, 20, 30, 40, 50, 60, 70, 80)

for(i in xx){

  fit1 <- gam(y ~ s(x, bs = "cr", k = i), data = d)

  fit2 <- summary(fit1)

  edfs <- c(edfs, fit2$edf)

  gcvs <- c(gcvs, fit1$gcv.ubre)

  aics <- c(aics, fit1$aic)

  if(fit2$s.pv > alpha){

    tmp <- 4#有意でないなら×シンボル

  } else {

    tmp <- 16有意でないなら●シンボル

  }

  pchs <- c(pchs, tmp)

}

plotn(edfs ~ xx, pch = pchs)#図28の描画 

plotn(gcvs ~ xx, pch = pchs)#図29の描画 

plotn(aics ~ xx, pch = pchs)#図30の描画 

------------------------------------------------------

ノット数が増えるにつれて、同様にサチって行くものの、どの段階でもsmooth termは有意にならない。こういう時は、smooth termを除去したほうが良いだろう。実際、切片項だけの線形回帰にするとAICは435.8程度になり、最も予測性の良いモデルとなる。もちろん、見たい現象があるときに、不正にならない範囲でノットを調整する必要はあるだろう。

2変数のGAM―線形＆非線形―

　続いて、2変数が含まれるGAMにチャレンジしよう。

------------------------------------------------------

x1 <- c(9.7, 14.3, -5.3, 13.1, 9.5, 7.2, -7.5, -5.2, 4.9, -7.4, -9.1, -1.2,

        -6.7, -10.2, 1.6, 6.7, -6.4, -14.7, -12.1, -8.6, -6.4, -9.5, -7.4, 

        13.3, 1.1, 8.2, -1.6, 14.2, -8, -12, 7.9, 11.1, -1.9, 7.9, -9.1, 12, 

        9, 9.7, 3, -11.5, 7.4, 0.7, -5.2, 13.5, -2.1, -3.3, -10.4, -10.1, 

        -9.8, -11.7, 8.2, 0.6, 6.7, 10.4, 6.6, 8.5, -0.8, 8.3, -12.8, -0.1,

        -1, 10.8, -2, -1.2, 13.3, -14.7, 13.8, -2.7, -7.3, 13.5, 11.1, -10.5, 

        -8.4, -7.8, 14.6, 3.4, -2, 0.7, -1, 11.8, -11.1, -13.1, 7.8, -12.2, 

        13.9, 12.4, 14.7, 1.8, 8, -6.4, 10.9, 0.7, 13, -13.7, 11.8, 8.7, -0.8,

        -13, 9.7, -3.9)

x2 <- c(8.1, 6, 0.6, 12.6, 6.3, -5.1, 8.7, 6.6, -4.2, -7.9, 10.8, 8.7, -6.8, 

        0.5, 5.2, 2.4, 13.8, 9.6, -3.7, -10.3, -13.3, -7.8, 6.5, -9.5, -6.5, 

        -7.4, -5.1, 11.4, -4.3, 5.4, 6.8, -14.9, -11, 10.7, -7, 9, 6.7, 2.9, 

        -5.7, -13.6, 9.2, 7.3, 6.2, 4.6, -6.3, -13.6, -14.5, 10.2, 14.8, -5.5,

        8.1, -3, -8.5, -1.7, 6, 13.3, -3, 14.8, 1.6, -6.3, 11, -6.1, -6.1, 

        11.4, -10.9, 0.3, -4.4, -7.4, -14, 4.2, -8, -9.7, 10.9, 2.8, -8.9, 6.6,

        2.4, 2.4, 0.9, -11.6, 7.4, 4.4, -13, -3.6, 11.2, 1.2, -8.7, -9.9, -7.5,

        13.6, -9.6, 4.1, 1.1, -3.5, -1.2, -10.8, -4.5, 6.3, 5.7, -7.4)

y <- c(13.4, -19.4, -14.2, -27.8, 14.6, 6.9, 8, -20.7, -17.2, 15.1, 13.6, 7.7,

       7.2, 11.3, 8.1, -3.9, -11.9, -17, -13.3, 23.9, 4.5, 22.1, 6, -15.8, 

       19.2, 22.2, 9.6, -24.5, 21.5, -18.2, 9.2, 13, 12.7, 9.5, 23.8, -16.5,

       16, 13.9, -6.3, 0.5, 3, 14, -19.5, -22.5, 1.8, -9.5, 16.7, 6.7, 4.5, 

       -6.9, 14.2, 18.9, 1.8, 12.4, -8.3, 12.1, 16.9, 8.7, -22, 24.9, 9.6, 9.2,

       5.9, 5.4, -12.9, -8.8, -17, -2.8, 16.3, -21.9, 8.7, 12.6, 11.5, 13.6,

       -8.5, -17.2, -0.6, 15.9, 13.9, -0.5, -5, -23.9, 22.1, -14.3, -25.7, 

       -13.3, -9, 13.2, 19.9, -11.2, 9.4, 14.9, -19.8, -15.6, -7.9, 26.3, 

       17.4, -23.2, 12.9, -15)

d <- data.frame(x1 = x1,  x2 = x2, y = y)

plotn(y ~ x1, data = d)#図30の描画 

plotn(y ~ x2, data = d)#図31の描画 

------------------------------------------------------

明らかに、yとx1の関係は周期性（非線形）があるが、x2とは線形の関係がある。この時、gam関数は非線形な効果と線形な効果の両方を組み込むことが可能である。下記のように、非線形を仮定する変数はs関数を噛ませて、そうでない関数はそのままformulaに記述すればよい。

------------------------------------------------------

fit1 <- gam(y ~ s(x1, bs = "cr") + x2, data = d)

fit2 <- summary(fit1)

fit2

## 

## Family: gaussian 

## Link function: identity 

## 

## Formula:

## y ~ s(x1, bs = "cr") + x2

## 

## Parametric coefficients:

##             Estimate Std. Error t value Pr(>|t|)    

## (Intercept)  1.54193    0.50577   3.049  0.00303 ** 

## x2          -0.51068    0.06187  -8.254 1.31e-12 ***

## ---

## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 

## Approximate significance of smooth terms:

##         edf Ref.df    F p-value    

## s(x1) 8.972      9 77.9  <2e-16 ***

## ---

## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 

## R-sq.(adj) =  0.884   Deviance explained = 89.6%

## GCV = 28.728  Scale est. = 25.576    n = 100

------------------------------------------------------

上記では非線形な効果（x1）と線形の効果（x2）について両方が有意であると判定された。ではこれらを踏まえて回帰曲線、信頼区間と予測区間を図示しよう。重回帰で紹介したように、多変数関数の回帰曲線（曲面）の図示はやや面倒で、解釈もそれなりに困難を伴う。いくつかの図示の方法があるが、まず片方の変数を代表値に固定して、もう一方の変数の挙動を確認する図示を示す。非線形な効果を固定するときは、代表値の選び方に気を付けたほうが良い。今回は図の見た目から、大体データ全体の中央付近を通るようなx1の代表値を選んだ。

------------------------------------------------------

xx1 <- seq(min(d$x1), max(d$x1), length = 200)

xx2 <- mean(d$x2)

newdata <- data.frame(x1 = xx1, x2 = xx2)

pred <- predict(fit1, newdata = newdata, se.fit = T)

c_upr <- pred$fit + qt(p = 0.975, df = df.residual(fit1)) * pred$se.fit

c_lwr <- pred$fit - qt(p = 0.975, df = df.residual(fit1)) * pred$se.fit

p_upr <- pred$fit + qt(p = 0.975, df = df.residual(fit1)) * 

  sqrt(pred$se.fit^2 + fit1$sig2)

p_lwr <- pred$fit - qt(p = 0.975, df = df.residual(fit1)) * 

  sqrt(pred$se.fit^2 + fit1$sig2)

lim <- range(c(p_upr, p_lwr, d$y))

plotn(y ~ x1, data = d, ylim = lim)#図32の描画 

overdraw(polygon(c(xx1, rev(xx1)), c(c_upr, rev(c_lwr)), 

                 border = F, col = "#00000030"),

         polygon(c(xx1, rev(xx1)), c(p_upr, rev(p_lwr)), 

                 border = F, col = "#00000030"),

         points(xx1, pred$fit, type = "l", col = "red"))

xx1 <- -2.5

xx2 <- seq(min(d$x2), max(d$x2), length = 200)

newdata <- data.frame(x1 = xx1, x2 = xx2)

pred <- predict(fit1, newdata = newdata, se.fit = T)

c_upr <- pred$fit + qt(p = 0.975, df = df.residual(fit1)) * pred$se.fit

c_lwr <- pred$fit - qt(p = 0.975, df = df.residual(fit1)) * pred$se.fit

p_upr <- pred$fit + qt(p = 0.975, df = df.residual(fit1)) * 

  sqrt(pred$se.fit^2 + fit1$sig2)

p_lwr <- pred$fit - qt(p = 0.975, df = df.residual(fit1)) * 

  sqrt(pred$se.fit^2 + fit1$sig2)

lim <- range(c(p_upr, p_lwr, d$y))

plotn(y ~ x2, data = d, ylim = lim)#図33の描画 

overdraw(polygon(c(xx2, rev(xx2)), c(c_upr, rev(c_lwr)), 

                 border = F, col = "#00000030"),

         polygon(c(xx2, rev(xx2)), c(p_upr, rev(p_lwr)), 

                 border = F, col = "#00000030"),

         points(xx2, pred$fit, type = "l", col = "red"))

------------------------------------------------------

区間予測の表示をあきらめて、yとx1の関係についてx2が変化した時の回帰曲線の挙動、あるいはとx2の関係についてx1が変化した時の回帰曲線の挙動を図示する方法も示す。

------------------------------------------------------

l1 <- 200

l2 <- 10

xx1 <- seq(min(d$x1), max(d$x1), length = l1)

xx2 <- quantile(d$x2, probs = seq(0, 1, length = l2))

newdata <- data.frame(x1 = rep(xx1, times = l2), x2 = rep(xx2, each = l1))

pred <- predict(fit1, newdata = newdata, se.fit = T)

xxt <- quantile(d$x2, probs = seq(0, 1, length = 2*l2-1))

xxt <- xxt[!xxt%in%xx2]

d$g <- sapply(d$x2, function(tmp){

  LETTERS[sum(xxt < tmp) + 1]

})

lim <- range(pred$fit)

plotn(y ~ x1, data = d, ylim = lim, mode = "m", group = "g", legend = T, 

      leg.lab = xx2, leg.title = "x2", leg.sp = 4)#図34の描画 

for(i in 1:l2){

  tmp <- pred$fit[newdata$x2 == xx2[i]]

  overdraw(points(xx1, tmp, type = "l", col = .default_col[i]))

}

l1 <- 10

l2 <- 200

xx1 <- quantile(d$x1, probs = seq(0, 1, length = l1))

xx2 <- seq(min(d$x2), max(d$x2), length = l2)

newdata <- data.frame(x1 = rep(xx1, each = l2), x2 = rep(xx2, times = l1))

pred <- predict(fit1, newdata = newdata, se.fit = T)

xxt <- quantile(d$x1, probs = seq(0, 1, length = 2*l1-1))

xxt <- xxt[!xxt%in%xx1]

d$g <- sapply(d$x1, function(tmp){

  LETTERS[sum(xxt < tmp) + 1]

})

lim <- range(pred$fit)

plotn(y ~ x2, data = d, ylim = lim, mode = "m", group = "g", legend = T, 

      leg.lab = xx1, leg.title = "x1", leg.sp = 4)#図35の描画 

for(i in 1:l1){

  tmp <- pred$fit[newdata$x1 == xx1[i]]

  overdraw(points(xx2, tmp, type = "l", col = .default_col[i]))

}

------------------------------------------------------

おおむね、回帰曲線はその周りのデータの良い推定になっていそうだ。

　あるいは下記のように、回帰曲面を3次元グラフやヒートマップを使って表示することも可能である。

------------------------------------------------------

l <- 100

xx1 <- seq(min(d$x1), max(d$x1), length = l)

xx2 <- seq(min(d$x2), max(d$x2), length = l)

newdata <- data.frame(x1 = rep(xx1, times = l), 

                      x2 = rep(xx2, each = l))

pred <- predict(fit1, newdata = newdata)

mat_pred <- matrix(pred, nrow = l)

persp(xx1, xx2, mat_pred, 

      xlab = "x1", ylab = "x2", zlab = "predicted y",

      theta = 140, phi = 20, expand = 0.5, ticktype="detailed")#図36の描画 

image(xx1, xx2, mat_pred, col = viridis(256), xlab = "x1", ylab = "x2", las = 1)#図37の描画 

contour(xx1, xx2, mat_pred, add = T)

------------------------------------------------------

なお、上記のような手計算以外にも、vis.gam関数にgamオブジェクトを入力することで簡単に表示することもできる。

------------------------------------------------------

vis.gam(fit1, color="cm", theta = 140, phi = 20)#図38の描画 

vis.gam(fit1, color="cm", plot.type="contour")#図39の描画 

------------------------------------------------------

2変数のGAM―非線形＆非線形―

　続いて、非線形の2変数が含まれるGAMにチャレンジしよう。

------------------------------------------------------

x1 <- c(-4.9, -5.3, 15, 8.3, -7.8, -5.2, 5.8, 13.8, 9.8, -7.6, 1.4, -10.7, 14.5, 2.7, -8.9, -9.3, 13.3, 5.4, -11.5, 9.2, 0.5, 4, 5.9, -1, -4.5, 1.4, 5.1, -1.6, -3.9, 3.8, 11.8, 14.1, 7.4, 12.5, 10.1, 7.2, 9.5, 9, -1.3, -8.5, 5.1, 1.7, -6.9, 6.3, 12.7, -6.4, 9.5, 0.2, -2.7, 6.1, -2.3, -12.4, 0.2, 3.4, 9.4, -9.6, 10.9, 11.8, -6.5, 4.9, -0.2, -0.9, 7.9, 3.6, 14.1, 7.8, 0.2, -7.4, -13.6, 3.5, 1.5, 4.1, -9.8, 5.3, -7.3, -10.5, 0.7, 7.2, 9.5, 5, -5.4, 13.8, -9.3, -9.7, -0.5, 7, 13, -5.2, -9.9, 1, -12.7, 2.6, 2.6, -6.8, -2.3, -5.6, -14.4, -3.9, 14.5, 1.7)

x2 <- c(-10.9, 6.6, -12, 4.3, 11.1, -7.5, 10.1, -5.2, 13.9, 7.9, -5.3, 14.6, -6.6, -7.7, -11.9, -13.9, 3.8, -1, -11.2, 5.6, -10.7, -4.2, -11.8, -2.9, 14.3, 13.6, -12.3, -7.9, 8.6, -0.8, -6.2, 13.9, 10.4, -14.9, -11.8, 8.3, -11.7, 0.1, -14.1, -14.8, 12.9, 2, -14.4, -5.8, 13.7, -1.1, 14.6, 5.2, 2.1, 15, 14.7, -4.4, -10, 3.7, -7.1, -0.6, -1, 3.2, -4.4, -2.5, -14.5, -12.4, -5.7, 2.4, 2.1, -4.1, 9.1, -1.7, 1.9, 2.7, -4.6, -9.7, 1.4, -4.4, -14.9, 1.5, -11.4, -5.1, 14.3, -11.9, -8.3, 9.9, -7.6, 13.6, 4.4, 4.9, -10.3, -1.8, -7.8, -2.5, 1.8, -13.7, -14.2, 1.6, 10.1, 13.2, 10.5, 11.8, 13.3, 8.6)

y <- c(-8.7, -15, -3.2, 9.6, 1.3, -16.6, -17.1, -19.6, 15.8, -2.6, -4, 9.6, -17.8, -12.1, 15.1, 18.8, -6.2, 1, -3.3, 6.7, 9.5, -15.1, -4.3, 8.3, -3.4, 10.9, -6.1, -6.1, -18.1, 0.4, -14.9, -2.6, -5, 3.4, 12, -6.4, 11.9, 21.1, 12.8, 17.7, -8.7, 11.4, 11.8, -12.6, -2.6, 7.6, 17.1, 8.8, 3.1, 4.2, 8.3, -14.5, 5.5, -3.7, -0.4, 17.9, 11.9, 2.5, -5.7, -6.4, 19.9, 12.6, -2, -1.5, -1.3, 2.3, 0.3, 11.8, -2.5, -1.4, -0.3, -14.2, 15.5, -16, 15.4, 11.8, 10.5, -6, 15, -6.2, -12.7, -20, -0.2, 12.8, 12.2, 0.8, -13.3, -2.8, -2.2, 10.9, -1.8, 7.1, 8.4, 10.9, -10, -5.3, -14.2, -12, -4.8, -5)

d <- data.frame(x1 = x1,  x2 = x2, y = y)

plotn(y ~ x1, data = d)#図40の描画 

plotn(y ~ x2, data = d)#図41の描画 

------------------------------------------------------

明らかに、yはx1、x2の両方について非線形の関係がある。両方とも（交互作用のない）非線形を仮定するなら、下記のように解析できる。

------------------------------------------------------

fit1 <- gam(y ~ s(x1, bs = "cr") + s(x2, bs = "cr"), data = d)

fit2 <- summary(fit1)

fit2

## 

## Family: gaussian 

## Link function: identity 

## 

## Formula:

## y ~ s(x1, bs = "cr") + s(x2, bs = "cr")

## 

## Parametric coefficients:

##             Estimate Std. Error t value Pr(>|t|)  

## (Intercept)   0.5570     0.2236   2.491   0.0148 *

## ---

## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 

## Approximate significance of smooth terms:

##         edf Ref.df     F p-value    

## s(x1) 8.964  8.999 113.5  <2e-16 ***

## s(x2) 8.300  8.838 108.3  <2e-16 ***

## ---

## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 

## R-sq.(adj) =  0.956   Deviance explained = 96.4%

## GCV = 6.1194  Scale est. = 5.0017    n = 100

------------------------------------------------------

yとx1の関係についてx2が変化した時の回帰曲線の挙動、あるいはとx2の関係についてx1が変化した時の回帰曲線の挙動を図示する方法を示す。

------------------------------------------------------

l1 <- 200

l2 <- 10

xx1 <- seq(min(d$x1), max(d$x1), length = l1)

xx2 <- quantile(d$x2, probs = seq(0, 1, length = l2))

newdata <- data.frame(x1 = rep(xx1, times = l2), x2 = rep(xx2, each = l1))

pred <- predict(fit1, newdata = newdata, se.fit = T)

xxt <- quantile(d$x2, probs = seq(0, 1, length = 2*l2-1))

xxt <- xxt[!xxt%in%xx2]

d$g <- sapply(d$x2, function(tmp){

  LETTERS[sum(xxt < tmp) + 1]

})

lim <- range(pred$fit)

plotn(y ~ x1, data = d, ylim = lim, mode = "m", group = "g", legend = T, 

      leg.lab = xx2, leg.title = "x2", leg.sp = 4)#図42の描画 

for(i in 1:l2){

  tmp <- pred$fit[newdata$x2 == xx2[i]]

  overdraw(points(xx1, tmp, type = "l", col = .default_col[i]))

}

l1 <- 10

l2 <- 200

xx1 <- quantile(d$x1, probs = seq(0, 1, length = l1))

xx2 <- seq(min(d$x2), max(d$x2), length = l2)

newdata <- data.frame(x1 = rep(xx1, each = l2), x2 = rep(xx2, times = l1))

pred <- predict(fit1, newdata = newdata, se.fit = T)

xxt <- quantile(d$x1, probs = seq(0, 1, length = 2*l1-1))

xxt <- xxt[!xxt%in%xx1]

d$g <- sapply(d$x1, function(tmp){

  LETTERS[sum(xxt < tmp) + 1]

})

lim <- range(pred$fit)

plotn(y ~ x2, data = d, ylim = lim, mode = "m", group = "g", legend = T, 

      leg.lab = xx1, leg.title = "x1", leg.sp = 4)#図43の描画 

for(i in 1:l1){

  tmp <- pred$fit[newdata$x1 == xx1[i]]

  overdraw(points(xx2, tmp, type = "l", col = .default_col[i]))

}

------------------------------------------------------

回帰曲面を3次元グラフやヒートマップを使って表示する。

------------------------------------------------------

l <- 100

xx1 <- seq(min(d$x1), max(d$x1), length = l)

xx2 <- seq(min(d$x2), max(d$x2), length = l)

newdata <- data.frame(x1 = rep(xx1, times = l), 

                      x2 = rep(xx2, each = l))

pred <- predict(fit1, newdata = newdata)

mat_pred <- matrix(pred, nrow = l)

persp(xx1, xx2, mat_pred, 

      xlab = "x1", ylab = "x2", zlab = "predicted y",

      theta = 140, phi = 20, expand = 0.5, ticktype="detailed")#図44の描画 

image(xx1, xx2, mat_pred, col = viridis(256), xlab = "x1", ylab = "x2", las = 1)#図45の描画 

contour(xx1, xx2, mat_pred, add = T)

------------------------------------------------------

2変数のGAM―非線形＆非線形かつ交互作用あり―

　最後に、非線形の2変数が含まれ、かつ交互作用を考慮したGAMにチャレンジしよう。

------------------------------------------------------

x1 <- c(4.7, -7.2, -12.5, 12.6, -1.5, -4.8, 12.1, -1.6, -5.1, -5.2, -0.3, 5.4, -6.7, -3.6, 4.9, 2.1, 0.7, 4, -10.1, 13.8, 8.6, 5.8, 11.2, 1, 5, -6, -2.7, -14.4, 11.9, 4.7, 11.2, -7.2, 9, -5.5, 6.7, 3.4, 7.7, -4.8, -3.1, -12.1, -5, -4.8, 7.1, -6.9, -7.6, -13.4, 3.1, -2.4, -2.7, 8.4, 13.5, 0.2, -1.8, 7.3, 13.3, -7.4, -13.3, 7.9, -1.5, 10.7, 12.6, -8, -13.3, 3.7, 10.8, -11.4, -3.4, 8.4, -13.9, -13.6, 7.5, 9.2, -12.2, -7.1, -9.9, 14.6, -2.7, 14, -6.6, 12.3, 8.4, -10.9, -4.2, -7.1, 4.4, -1, -11.4, -4.3, -7.6, 5, 6, -9, 5.4, 7, 10.9, -10.9, 14.1, 12.8, -11.3, 6)

x2 <- c(-12.6, 1, -9.3, 7.4, 0.1, -3.4, 10.2, 9, 5, -5, 8.5, 7.8, 4.5, -11.3, 9, 7, 14.5, 9.1, 1.9, 13.3, -5.4, 4.1, -7.2, 3.9, -4.3, -5.5, -0.3, 9, 6.8, 1.2, 4.6, -9.3, -3.9, 1.1, 14.1, 1.7, 14.5, -12, -5.6, 10.7, 0.7, -1.5, -4.8, 14, -9.4, -11.2, 0.6, -9.4, -10.7, -10.1, 12.6, 7.8, -13, -7.4, 14.6, 3.4, -2.5, 0.6, -5.7, -5.7, -4.4, -7.2, -6.5, 10.9, 11.5, 13.3, -1.2, -10, 5.9, -14.5, -6.8, -12.8, 14.3, -6.1, 5.1, 11.1, 10.9, -12.3, 5, -14.2, 0.4, -11.4, 11.5, -14.9, 2, -3.1, -2.7, -11.5, -0.4, -8.3, -0.6, -2.9, -9, 10.9, -11, -4.8, 5.2, 2.9, -6.3, 12.2)