里の桜、都会の桜

追分市民の森の菜の花畑と山桜









ドウダンツツジ


森の中で、ヒトリシズカが咲きだしていました。




目黒川の桜




花筏

---

花筏とか見てたら、一般化線形モデル（ Generalized linear model ）を、かっこよく表現したいなぁと思って
こんなの作りました。

Ｒのcarsデータ・セット

こんなのを、
こんな感じに

評判の悪い３D円グラフと同じで、見た目のかっこ良さを第一の目的にしましたが、このモデルの予測データの確率分布を感じていただければよいかと

rglというパッケージを使うと、アニメーションGIFやWebGLにもできます。





WebGLの出力例はこちら（804KB) 、スマホのブラウザーでは動かないかも


一般化線形モデルって、せっかくカッコイイ名前がついているのに、自分の出番を、こんな感じで

簡単な直線回帰モデルでごまかされちゃったり...
散布図に回帰直線（上図中の赤線）引かれると、それだけでなんとなく納得しちゃいそうですけど、よく見るとXの値が小さいときの予測値の95%信頼区間（上図中の破線）があやしかったり、回帰直線のまわりの残差の分散も不均一です。



「マンガ 確率・統計が驚異的によくわかる」という入門書にも残差の分散が不均一のときは、モデルを訂正してくださいと書いてあります。（この本では分散不均一を ヘテロセダスティック : Heteroscedasticと書いています。）

ということで、こんなときはアスピリン２錠ではなく、かっこいい一般化線形モデルです。

library(MASS)
require(graphics)

#---
# linear model

par(mfrow=c(1,1))

plot(cars, xlab="speed (mph)", ylab="distance (ft)")
mdl.lm=lm(dist~speed,data=cars)

lm.predict <- predict(mdl.lm, interval="prediction")
lines(cars$speed, lm.predict[,2], lwd=1, col="blue", lty = 2)
lines(cars$speed, lm.predict[,3], lwd=1, col="blue", lty = 2)
lines(cars$speed, lm.predict[,1], lwd=2, col="red")

par(mfrow=c(2,2))
plot(mdl.lm)

#---
# Generalized Linear Model:Gamma Regression

par(mfrow=c(1,1))

mdl.glm <- glm(dist ~ I(speed + speed^2), family=Gamma(link="identity"),
    data = cars)
mdl.glm2 <- glm(dist ~ I(speed^2), family=Gamma(link="identity"), data = cars)
mdl.glm3 <- glm(dist ~ speed, family=Gamma(link="log"), data = cars)

a <- gamma.shape(mdl.glm)$alpha
x <- seq(1, max(cars$speed), by = 2)
xd <- data.frame(speed = x)
mdl.pred <- predict(mdl.glm, xd, type="response")
scale <- mdl.pred / a
mdl.lower <- qgamma(0.025, shape = a, scale=scale)
mdl.upper <- qgamma(0.975, shape = a, scale=scale)

plot(cars$speed, cars$dist, xlim=c(0, max(cars$speed) * 1.1), 
    ylim=c(0, max(cars$dist) * 1.1), 
    xlab="speed (mph)", ylab="distance (ft)")
lines(x, mdl.pred, col="red")
lines(x, mdl.lower, col="blue", lty = 2)
lines(x, mdl.upper, col="blue", lty = 2)

par(mfrow=c(2,2))
plot(mdl.glm)

####
# Generalized Linear Model:Gamma Regression (plot3d)
#
library(MASS)
library(rgl)

par(mfrow=c(1,1))

mdl.glm <- glm(dist ~ I(speed + speed^2), family=Gamma(link="identity"), 
    data = cars)

a <- gamma.shape(mdl.glm)$alpha

x <- seq(1, max(cars$speed))
xd <- data.frame(speed = x)
mdl.pred <- predict(mdl.glm, xd, type="response")
scale <- mdl.pred / a
mdl.lower <- qgamma(0.025, shape = a, scale=scale)
mdl.upper <- qgamma(0.975, shape = a, scale=scale)

y <- seq(1, max(cars$dist))

pd <- matrix(0, length(x), length(y))

for(i in x) {
  pd[i,] <- dgamma(y, shape=a, scale=scale[i])
}

pd.z <- dgamma(y, shape=a, scale=scale)

open3d(scale=c(5,1,300))
  
# clear scene:
clear3d("all")

spheres3d(cars$speed, cars$dist,rep(0,length(cars$speed)),
    radius=1,color="black")
surface3d(x,y, pd,alpha=0.1,front="lines")
for(i in 1:length(x) - 1) {
  segments3d(c(x[i],x[i+1]), c(mdl.pred[i], mdl.pred[i+1]), z=rep(0, 2), col=2, lwd=2)
  segments3d(c(x[i],x[i+1]), c(mdl.lower[i], mdl.lower[i+1]), z=rep(0, 2), col=3, lwd=1)
  if(mdl.upper[i+1] <= max(y)) {
    segments3d(c(x[i],x[i+1]), c(mdl.upper[i], mdl.upper[i+1]), z=rep(0, 2), col=3, lwd=1)
  }
}

axis3d("x")
axis3d("y")
title3d(xlab="speed (mph)", ylab="distance (ft)")
surface3d(x,y, pd,col=c("green", "blue")[1+(pd > 0.001)],alpha=0.3,front="lines")

#writeWebGL()

mdl.glm2 <- glm(dist ~ I(speed^2), family=Gamma(link="identity"), data = cars)
mdl.glm3 <- glm(dist ~ speed, family=Gamma(link="log"), data = cars)
#EOF

#いくつかのモデルでAICを比べてみる。
> library(MASS)
> library(rgl)
> 
> mdl.glm <- glm(dist ~ I(speed + speed^2), family=Gamma(link="identity"), data = cars)
> 
> mdl.glm2 <- glm(dist ~ I(speed^2), family=Gamma(link="identity"), data = cars)
> mdl.glm3 <- glm(dist ~ speed, family=Gamma(link="log"), data = cars)
> AIC(mdl.glm, mdl.glm2, mdl.glm3)
         df      AIC
mdl.glm   3 410.7229
mdl.glm2  3 411.0588
mdl.glm3  3 415.5611

マンガ 確率・統計が驚異的によくわかる

• 作者: ラリー ゴニック
• 出版社/メーカー: 白揚社
• 発売日: 1995/12
• メディア: 単行本