草薙の研究ログ

英語の先生をやってます。

GAMLSSとBAMLSSで分布をデータにフィット

GAMLSSとBAMLSS

新幹線に乗ってコンビニにいくようなものだけど,GAMLSSとそのベイズ的モデルであるBAMLSSを使って,観測のデータに対して分布をフィットさせる。
GAMLSSというのは,「位置,尺度,形状,それぞれの母数のための一般化加法モデル(Generalized additive model for location, scale and shape)」のこと。詳しくはwikiや,開発者らのページを。

Generalized additive model for location, scale and shape - Wikipedia

gamlss | for statistical modelling

大雑把にいえば,分布ベースの回帰モデルのひとつ。
基本的には,(ひとつの)応答変数になんらかの分布を仮定する。正規分布とかガンマとか。
たとえばガンマ分布だったら,α(形状)とβ(比率)母数のふたつがあると。
で,応答変数のそれぞれの値は,これらの母数によってランダムに発生していると。
GAMLSSでは,これら,それぞれの母数の値自体を予測する。
予測変数の値という条件における,応答変数が従うであろう分布の母数を予測して,そこから応答変数ができてるみたいな。

この母数に対する予測は,線形な関係であってもいいし,加法モデルというくらいだから,スプラインとかでもいい。
基本,4母数までの分布ならできるらしく(μ,σ,ν,τ),Rのgamlssパッケージには莫大な数の分布がサポートされている。
変量効果も入れられる。とっても自由。

BAMLSSというのは,これとほぼ同じようなことをベイズ的にやるんだそうだ(細かい事前分布の設定とかアルゴリズムはよくわからない)。

これが面白いのは,平均値を予測するだけでなくて,散布度や形状を予測することができるようになるところ。
たとえば,実験条件Aは,応答変数の平均値には影響を及ぼさないけど,分散を大きくする効果をもつ,とか。
歪んだ分布を仮定して,実験条件Bは,応答変数の形状に影響を及ぼす,とか。そういうモデリングができるようになる。
「分位点回帰のパラメトリック版を非線形でもできるようにしてる」,みたいな感じ。

…それは置いておいて,このパッケージは特に便利なので,まずは単変量の分布の母数をこれらのパッケージで推定する。

ある分布に従う単変量の母数の推定

MASSパッケージのfitdistrやfitdistrplusパッケージでもoptim関数でも,なんなら場合によってはglm関数でも,同じことをやることができる場合は多いのだけど,gamlssパッケージやbamlssパッケージは,普通に単変量のデータの分布を推定できる。
まずは,データにワイブル分布を当てはめて母数を推定する例をやってみる。

#読み込み
library(gamlss)
#データの作成
set.seed(92)
d1<-rweibull(1000,2,1)

#gamlssの基本関数を使って推定
#最初の回帰式がμ(最初の母数),次の回帰式をσと当てている使用
WEIは,ワイブル分布
m1<-gamlss(d1~1,sigma.formula=~1,family="WEI")

#これでダン
#サマリーしてみる

summary(m1)

切片の推定値はそれぞれ,μ(ここでは尺度母数) = -0.01,σ(ここでは形状母数) = .677くらい。
これにはリンク関数としてlogが噛んでいるので,係数をexpする。
すると,まあ,1と2くらい。よく推定できている。もちろん,これはただの最尤推定と同じなので,他のパッケージとかとも一致する。

で,もっと楽な関数があって,gamlssMLという最尤推定の関数。

m2<-gamlssML(d1,"WEI")
summary(m2)

これもおなじことになる。わざわざ切片からの回帰式とか書かなくてもよい。

もっと便利なのがあって,histDistという関数。
これは,ヒストグラムとフィットした確率密度曲線を描いてくれるだけでなくて,上のサマリーも返してくれる。

histDist(d1,"WEI")

f:id:kusanagik:20170919175751p:plain

図はこんな感じ。オプションもあるので,色々いじれるけど,自分はやらないだろうな。

BAMLSSでも基本は同じ。少なくともこの状態で,それぞれの事前分布は基本的に無情報っぽい。
MCMCサンプラーJAGSとか色々指定できるみたいだけど,デフォルトではGMCMCというのを使っているそう(これはよくわからない)。

library(bamlss)

#いろいろなオプションを省略して…
m3<-bamlss(list(mu=d1~1,sigma=~1),family="weibull")
summary(m3)

それぞれの母数の95%信用区間(パーセンタイル法)も出してくれる。
以下のようにしたら,MCMCサンプルも可視化できる。

par(mfrow=c(1,2))
plot(density(exp(unlist(m3$s[,1]))),main="Scale",xlab="Value")
plot(density(exp(unlist(m3$s[,5]))),main="Shape",xlab="Value")

f:id:kusanagik:20170919194203p:plain

簡単なモデル

えっと,せっかくだから,ただ生データをフィットさせるだけでなくて,GAMLSSやBAMLSSの中でも,もっとも簡単な回帰モデル的なものをやってみよう。
アイリスデータの,がく片の長さ(sepal length)を応答変数,あやめの種類(3水準)を予測変数としよう。がく片の長さは正規分布に従うとして,種類は,がく片の長さの平均と標準偏差にそれぞれどのような影響をあたえるだろうか。ここであやめの種類はfactor型。
データとしてはこんな感じ。

f:id:kusanagik:20170919195037p:plain

見た感じ,平均はsetosa < versicolor < virginicaという大きさのようになっている。標準偏差も同じ。
普通の線形モデルなら,

summary(lm(iris[,1]~iris[,5]))
anova(lm(iris[,1]~iris[,5]))

みたいにする。

一方,gamlssでは,こんな感じにする。

m4<-gamlss(iris[,1]~iris[,5],sigma.formula=~iris[,5],family="NO")

summary(m4)

******************************************************************
Family:  c("NO", "Normal") 

Call:  gamlss(formula = iris[, 1] ~ iris[, 5], sigma.formula = ~iris[,  
    5], family = "NO") 

Fitting method: RS() 

------------------------------------------------------------------
Mu link function:  identity
Mu Coefficients:
                    Estimate Std. Error t value Pr(>|t|)    
(Intercept)          5.00600    0.04935  101.44   <2e-16 ***
iris[, 5]versicolor  0.93000    0.08751   10.63   <2e-16 ***
iris[, 5]virginica   1.58200    0.10179   15.54   <2e-16 ***
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

------------------------------------------------------------------
Sigma link function:  log
Sigma Coefficients:
                    Estimate Std. Error t value Pr(>|t|)    
(Intercept)          -1.0528     0.1000 -10.528  < 2e-16 ***
iris[, 5]versicolor   0.3814     0.1414   2.697  0.00783 ** 
iris[, 5]virginica    0.5900     0.1414   4.172  5.2e-05 ***
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

------------------------------------------------------------------
No. of observations in the fit:  150 
Degrees of Freedom for the fit:  6
      Residual Deg. of Freedom:  144 
                      at cycle:  2 
 
Global Deviance:     206.9715 
            AIC:     218.9715 
            SBC:     237.0353 
******************************************************************

普通の線形モデルと違って,標準偏差の大きさの違いも同時にモデリングしている。
まあ,見た感じの結果になっている。これは,こんな感じに可視化するとわかりやすい。
term.plotという関数もある。

par(mfrow=c(1,2))
term.plot(m4,"mu",main="Mu")
term.plot(m4,"sigma",main="Sigma")

f:id:kusanagik:20170919195536p:plain

もちろん,診断プロットもある。

plot(m4)

それぞれの予測された母数から,確率密度曲線を描くとこんな感じ。

f:id:kusanagik:20170919200250p:plain

もちろん,標準偏差には影響を及ぼさないとか,標準偏差には影響を及ぼすけど,平均値には及ぼさない,といったモデルを作り,適合度を比較したり尤度比検定をしたりすることもできる。
まあでも,これ予測変数が因子型だったら,教師ありの混合分布モデルみたいなもん。

blmssパッケージでは,

m5<-bamlss(list(mu=iris[,1]~iris[,5],sigma=~iris[,5]),family="gaussian")

で,(当然だけど)同じような結果を得ることができる。MCMCサンプルを可視化したり,信用区間を構築してもよい。
とっても便利。

平均は変わらないけど,標準偏差だけ違う2群のモデル

回帰モデルは自由に作れるので,ある因子が標準偏差だけに影響を及ぼしているモデルとかも作れる。
たとえば,こんなデータだとしよう。

f:id:kusanagik:20170919205647p:plain

group<-as.factor(c(rep("A",100),rep("B",100)))
score<-c(rnorm(100,50,10),rnorm(100,50,30))
d2<-data.frame(group,class)

m6<-gamlss(score~1,sigma.formula=~group,data=d2)
summary(m6)

term.plot(m6,"sigma",main="Sigma")

結果はこんな感じ。

f:id:kusanagik:20170919210104p:plain

ちょっと複雑なモデル

えっと,実質的に意味があるかはわからないけど(ないだろう),がく片の長さの分布を正規分布だと仮定して,がく片の幅と花びらの長さから,その予測値という条件下における分布の母数(平均,標準偏差)を推定するとしよう。で,その予測値と予測される母数の関係は非線形だと。

#ここでは罰則付きスプラインの一種
m7<-gamlss(Sepal.Length~pb(Sepal.Width)+pb(Petal.Length),sigma.formula=~pb(Sepal.Width)+pb(Petal.Length),data=iris,family="NO")
summary(m6)


par(mfrow=c(2,2))
term.plot(m7,"mu",main="Mu")
term.plot(m7,"sigma",main="Sigma")

可視化するとこんな感じ。

f:id:kusanagik:20170919203712p:plain


がく片の幅と花びらの長さが長くなれば,がく片の長さの分布における平均はそれぞれ単調に増加しそう(花びらの長さはちょっと指数的に見える)。
がく片の幅は,がく片の長さの分布における標準偏差とはあまり関係なさそうだ。しかし,花びらの長さは,がく片の長さの分布における標準偏差非線形の影響を及ぼしているかも(意味ないと思うけど)。…みたいな。

で,これに種類を変量効果と考えて…というようにやっていくと,大分自由なモデリングができそう。
最尤推定による分布のフィットから書き始めて,やっとここまで来たって感じ)