今度はバブルチャートと一般化線形モデルをやってみようと思います。
データは相変わらずirisを使います。
狙いは、一般化線形モデルを2項分布でやるというもので、具体的にはirisのPetal.WidthからSpceciesを推定します。
データの準備
x軸にPetal.Width、y軸にSpeciesを表示させます。
箱ひげ図
まずは基本として箱ひげ図。
g<-ggplot(iris,aes(x=Petal.Width, y=Species))+
geom_boxplot()
plot(g)
散布図
続いて箱ひげ図に散布図も入れてみる。ついでに箱ひげ図、散布図ともに種ごとの色を設定。
g<-ggplot(iris,aes(x=Petal.Width, y=Species))+
geom_boxplot(aes(color=Species))+
geom_point(aes(color=Species))
plot(g)
ただ、単純に散布図だとプロットが重なってしまって、分布具合がわかりにくくなってしまいます。
バブルチャート
重なってる数に合わせて点の大きさを変えるという形で散布具合を可視化するチヤートをバブルチャートと言います。
それは以下のようにできます。
g<-ggplot(iris,aes(x=Petal.Width, y=Species))+
geom_boxplot(aes(color=Species))+
geom_count(aes(color=Species),show.legend=FALSE)
plot(g)
使っているgeom_countはhelpを見ると、「geom_pointのヴァリエーションの1つで、各ロケーションでの観測点数を数え上げて表示する」とあります。
ちなみに、今回のデータは完全な生データなので、geom_countを使って描画の際に数え上げを行っていますが、最初から度数分布が与えられている場合には、geom_point()でsizeに度数が入った変数を指定することでバブルチャートを描くことができます。
一般化線形モデルをやる
一般化線形モデルとは
普通の線形モデリング（＝回帰分析）は、観測された目的変数は正規分布をしてるということが前提となっています。
一方、一般化線形モデリングでは、観測された目的変数の分布を以下の任意の分布に設定することができます。
| 分布名 | family(デフォルトのリンク関数) | 設定可能な他のリンク関数 |
|---|---|---|
| 二項分布 | binomial(link = "logit") | probit, cauchit, log, cloglog |
| 正規分布（ガウス分布） | gaussian(link = "identity") | log, inverse |
| ガンマ分布 | Gamma(link = "inverse") | identity and log |
| 逆ガウス分布（ワルド分布） | inverse.gaussian(link = "1/mu^2") | inverse, identity, log. |
| ポアソン分布 | poisson(link = "log") | identity, sqrt |
| 疑似尤度 | quasi(link = "identity", variance = "constant") | logit, probit, cloglog, inverse, log, 1/mu^2, sqrt, power() |
| 疑似二項分布 | quasibinomial(link = "logit") | 同上 |
| 疑似ポアソン分布 | quasipoisson(link = "log") | 同上 |
（正直、疑似尤度、疑似に項分布、疑似ポアソン分布はまだ勉強不足で良くわからないです。（苦笑））
また、線形モデリングでは期待値$E(y)=\alpha + \beta x$という単純な直線の式でフィッティングをします。
それに対して、GLMでは、目的変数の期待値をリンク関数Gで変換したものに対して、直線の式をフィッティングします。
つまり、$G(E(y)) =\alpha+\beta x$ということです。
分布をガウス分布にし、リンク関数をidentity（同一）に設定すると通常の回帰分析になります。
r.lm <- lm(Sepal.Width~ Sepal.Length,iris)
r.lm$coefficients
# ----------------------------------
# (Intercept) Sepal.Length
# 3.4189468 -0.0618848
# ----------------------------------
r.glm <-glm(Sepal.Width~Sepal.Length,family = gaussian(link = "identity"),iris)
r.lm$coefficients
# ----------------------------------
# (Intercept) Sepal.Length
# 3.4189468 -0.0618848
# ----------------------------------
AIC(r.lm,r.glm)
# ----------------------------------
# df AIC
# r.lm 3 179.4644
# r.glm 3 179.4644
# ----------------------------------
注意：変数変換とリンク関数による変換の違い
Rのlmに与えるformulaには変換を加えることができます。例えば、以下の例は観測変数に対数変換を掛けたものです。
r.lm <- lm(log(Sepal.Width)~ Sepal.Length,iris)
r.lm$coefficients
# (Intercept) Sepal.Length
# 1.20866016 -0.01732247
一方、glmで分布をガウス分布、リンク関数をlogに設定すると以下のようになります。
r.glm <-glm(Sepal.Width~Sepal.Length,family = gaussian(link = "log"),iris)
r.lm$coefficients
# -----------------------------
# (Intercept) Sepal.Length
# 1.23902617 -0.02081587
ということで、似てはいるものの、結果が微妙に異なっています。
上記したように、
回帰分析では$E(y)=\alpha + \beta x$なので、変数変換するということは
$log(E(y))=\alpha + \beta x$。
一方、リンク関数をlogとしたGLMでは
$log(E(y)) =\alpha+\beta x$。
両者は見た目同じになります。
ただ、両者は根本的に違っているのは、変数変換はあくまで観測された値に対して掛けるということです。
つまり、変数変換は、誤差も含めて変換されるのに対して、
GLMではあくまで期待値部分だけを変換し、誤差は変換の外に置かれているということです。
すなわち、数式的には、
回帰分析での対数変換は正しくは
$log(y)=log(E(y)+\epsilon)=\alpha+\beta x+\epsilon$
$\therefore y=e^{\alpha+\beta x+\epsilon}=e^{\alpha+\beta x}\cdot e^{\epsilon}$
というフィッティングを掛けることになります。
一方、リンク関数はあくまで期待値に対して掛けるものになります。
したがって、この場合の観測変数は
$y=G^{-1}(E(y))+\epsilon=e^{\alpha+\beta x}+\epsilon$
となります。（$G^{-1}$は逆関数という意味です）
ということで、両者で誤差項$\epsilon$の扱いが異なっています。
これは、誤差が従っている分布が異なっているということです。
つまり、前者の場合、$e^{\epsilon}$が正規分布に従う→$\epsilon$は対数正規分布に従うわけですが
一方、後者の場合、あくまで$\epsilon$はfamilyで設定したガウス分布（正規分布）に従うとうことになります。
どちらを使うべきかは、誤差としてどういうものを想定するのかによります。
参考：https://mrkm-a.hatenablog.com/entry/20140513/p1
花びらの幅から種を予想するGLM
ということで、もともとの狙いである、irisのPetal.WidthからSpeciesを予想するGLMを行います。
とりあえずデモということで、データはversicolorとvirginicaの2つに絞ることにします。ということで分布は二項分布になります。
# setosaだけ除いたデータにしたうえで、新たにSという変数を作って、
# versicolorを0、virginicaを1というダミー数字を割り当てる
iris.glm <- subset(iris,Species!="setosa")
iris.glm$S <- c(rep(0, 50),rep(1,50))
# 描画
g<-ggplot(iris.glm,aes(x=Petal.Width,y=S))+
geom_count(aes(color=Species))
plot(g)
# glm
r.glm <- glm(S~Petal.Width, family = binomial, iris.glm)
結果は
> summary(r.glm)
Call:
glm(formula = Species ~ Petal.Width, family = binomial, data = iris.glm)
Deviance Residuals:
Min 1Q Median 3Q Max
-2.13868 -0.16468 -0.00929 0.13000 2.46891
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -21.126 4.597 -4.596 4.31e-06 ***
Petal.Width 12.947 2.843 4.554 5.25e-06 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 138.629 on 99 degrees of freedom
Residual deviance: 33.421 on 98 degrees of freedom
AIC: 37.421
Number of Fisher Scoring iterations: 7
となります。
なお、今回、versicolorに0、virginicaに1を割り当てましたが、これはggplotでの描画のためのものであって、glmを実行だけなら明示的に割り当てしなくても、glmの方の内部処理として勝手に0と1が割り当てられます。
結果の信頼区間を求める。
confint()という関数で、パラメータの信頼区間を算出できます。
第1引数はglmの結果、第2引数はlevel=という形で信頼区間の範囲を0から1までの数値で与えます。
coef(r.glm)
# (Intercept) Petal.Width
# -21.12559 12.94748
confint(r.glm)
# Waiting for profiling to be done...
# 2.5 % 97.5 %
# (Intercept) -32.144431 -13.68017
# Petal.Width 8.355608 19.79404
参考：http://triadsou.hatenablog.com/entry/20101022/1287715244
GLMの結果の描画 その1
ggplotのgeom_smoothを使うと簡単に描画できます。
2つ目の引数がゴチャゴチャしてますが、glmに渡す分布を指示しています。
g<-ggplot(iris.glm,aes(x=Petal.Width,y=S))+
geom_count(aes(color=Species))+
geom_smooth(method = "glm",method.args = list(family = binomial(link = "logit")))
plot(g)
GLMの結果の描画 その2
任意の関数を描画できるstat_functionを使って、明示的に関数を与えて描画します。
以下ではgeom_smoothと結果が一致していることが分かるように、geom_smoothに重ねて描画しています。
r.glm.coef <- coef(r.glm)
g<-ggplot(iris.glm,aes(x=Petal.Width,y=S))+
geom_count(aes(color=Species))+
geom_smooth(method = "glm",method.args = list(family = binomial(link = "logit")))+
stat_function(fun=function(x) 1/(1+exp(-(r.glm.coef[1]+r.glm.coef[2]*x))), color="red")
plot(g)
