久しぶりの青空 [相鉄沿線]
土曜日、久しぶりに青空が広がりました。
青い空に、真っ白な入道雲、このあとすごい夕立になりました。
今週もカブトムシに会えました。
ジメジメした日が続いていたからでしょうか、大きなコガネニカワタケ（黄金膠茸）が生えていました。
今週の夏休みの自由研究
BIOCONNECTOR ( http://bioconnector.org/workshops/r-stats.html#our_data:_nhanes ) のサイトにあるNHANES (the National Health and Nutrition Examination Survey )データ・セットを使って、身長と体重の関係を線形モデルと一般化線形モデルで表してました。
最初に、統計ソフトのサンプルプログラムにありがちな、身長と体重の線形モデル
散布図を表示して、線形モデルによる予測値の中央値を赤の実線、95%区間を青の破線で表します。
残差の評価
線形モデルを使って身長から予測した体重の値が大きくなるほど、残差も大きくなっています。身長が高い人たちほど体重の差は大きくなると考えられるので、この結果は容易に想像がつきます。
「マンガ 確率・統計が驚異的によくわかる」にも、「気づかないところからの手強い不意打ちの典型（身長/体重のデータの中にある）」と書かれています。
身長と体重の関係は、線形モデルでは、うまく表せないことは、わかりましたが、一般化線形モデルを使うとして、線形予測子や確率分布関数には、何を指定したら良いのか考えてみます。
身長と体重の関係を表すものとしては、BMI (body mass index) を算出する際に使われる「標準体重」の計算式 が広く知られているようです。
標準体重(kg) = 身長(m)2 × 22
このことから、体重は身長の自乗に比例すると考えるのが良さそうなので、線形予測子を
"Weight ~ I(Height^2) "
に、体重は0以上の連続値なので確率分布にはGamma分布を指定してモデルを作ります。
散布図を表示して、一般化線形モデルによる予測値の中央値を赤の実線、95%区間を青の破線で表します。
線形モデルと一般化線形モデルのAICを比較します。
体重が身長の2乗比例するモデルの方が、予測値の実測値に対するフィッティングも良く、AICも小さくなっています。
でも、まだ残差が少し不均一です。
BMIは、身長の低い人は値が小さくなり過ぎるため、英オックスフォード大学の数学者Nick Trefethen教授は、体重が身長の2.5乗に比例するNEW BMIという指標を提唱してるそうです。（参照 : NEW BMI https://people.maths.ox.ac.uk/trefethen/bmi.html ）
（身長が低い学童期ではローレル指数という、体重が身長の３乗に比例する計算式が使われているそうです。）
" Weight ~ I(Height^2.5) "としたときのAICも計算してみました。
AICの値は、2.5乗を使った方が、すこしだけ小さくなっています。
散布図でモデルのフィッティングを確認しましたが、２乗を使ったモデルとあまり変わらないようです。
もうすこし工夫しないといけないかな...
青い空に、真っ白な入道雲、このあとすごい夕立になりました。
今週もカブトムシに会えました。
ジメジメした日が続いていたからでしょうか、大きなコガネニカワタケ（黄金膠茸）が生えていました。
今週の夏休みの自由研究
BIOCONNECTOR ( http://bioconnector.org/workshops/r-stats.html#our_data:_nhanes ) のサイトにあるNHANES (the National Health and Nutrition Examination Survey )データ・セットを使って、身長と体重の関係を線形モデルと一般化線形モデルで表してました。
（NHANES データ・セットは、 https://ww2.amstat.org/publications/jse/v4n1/datasets.johnson.html からダウンロードできます。）
library(MASS)
require(graphics)
nh <- read.csv("datasets/nhanes.csv")
hw = nh[, c("Height", "Weight")]
# 欠損値(NA)を削除します。
hw = subset(hw, !is.na(hw$Height))
hw = subset(hw, !is.na(hw$Weight))
最初に、統計ソフトのサンプルプログラムにありがちな、身長と体重の線形モデル
散布図を表示して、線形モデルによる予測値の中央値を赤の実線、95%区間を青の破線で表します。
#--- # linear model par(mfrow=c(1,1)) plot(hw, xlab="height (cm)", ylab="weight (kg)") mdl.lm=lm(Weight~Height,data=hw) lm.predict <- predict(mdl.lm, interval="prediction") lines(hw$Height, lm.predict[,2], lwd=1, col="blue", lty = 2) lines(hw$Height, lm.predict[,3], lwd=1, col="blue", lty = 2) lines(hw$Height, lm.predict[,1], lwd=2, col="red")
残差の評価
par(mfrow=c(2,2)) plot(mdl.lm)
線形モデルを使って身長から予測した体重の値が大きくなるほど、残差も大きくなっています。身長が高い人たちほど体重の差は大きくなると考えられるので、この結果は容易に想像がつきます。
「マンガ 確率・統計が驚異的によくわかる」にも、「気づかないところからの手強い不意打ちの典型（身長/体重のデータの中にある）」と書かれています。
身長と体重の関係は、線形モデルでは、うまく表せないことは、わかりましたが、一般化線形モデルを使うとして、線形予測子や確率分布関数には、何を指定したら良いのか考えてみます。
身長と体重の関係を表すものとしては、BMI (body mass index) を算出する際に使われる「標準体重」の計算式 が広く知られているようです。
標準体重(kg) = 身長(m)2 × 22
このことから、体重は身長の自乗に比例すると考えるのが良さそうなので、線形予測子を
"Weight ~ I(Height^2) "
に、体重は0以上の連続値なので確率分布にはGamma分布を指定してモデルを作ります。
散布図を表示して、一般化線形モデルによる予測値の中央値を赤の実線、95%区間を青の破線で表します。
#---
# Generalized Linear Model:Gamma Regression
par(mfrow=c(1,1))
mdl.glm <- glm(Weight ~ I(Height^2), family=Gamma(link="identity"),
data = hw)
a <- gamma.shape(mdl.glm)$alpha
x <- seq(1, max(hw$Height), by = 2)
xd <- data.frame(Height = 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(hw$Height, hw$Weight, xlim=c(50, max(hw$Height) * 1.1),
ylim=c(0, max(hw$Weight) * 1.1), xlab="height (cm)", ylab="weight (kg)")
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)
線形モデルと一般化線形モデルのAICを比較します。
AIC(mdl.lm, mdl.glm) df AIC mdl.lm 3 41591.53 mdl.glm 3 40270.39
体重が身長の2乗比例するモデルの方が、予測値の実測値に対するフィッティングも良く、AICも小さくなっています。
でも、まだ残差が少し不均一です。
BMIは、身長の低い人は値が小さくなり過ぎるため、英オックスフォード大学の数学者Nick Trefethen教授は、体重が身長の2.5乗に比例するNEW BMIという指標を提唱してるそうです。（参照 : NEW BMI https://people.maths.ox.ac.uk/trefethen/bmi.html ）
（身長が低い学童期ではローレル指数という、体重が身長の３乗に比例する計算式が使われているそうです。）
" Weight ~ I(Height^2.5) "としたときのAICも計算してみました。
mdl.glm_NewBMI <- glm(Weight ~ I(Height^2.5),
family=Gamma(link="identity"), data = hw)
AIC(mdl.lm, mdl.glm, mdl.glm_NewBMI)
df AIC
mdl.lm 3 41591.53
mdl.glm 3 40270.39
mdl.glm_NewBMI 3 40111.38
AICの値は、2.5乗を使った方が、すこしだけ小さくなっています。
散布図でモデルのフィッティングを確認しましたが、２乗を使ったモデルとあまり変わらないようです。
もうすこし工夫しないといけないかな...
a_NewBMI <- gamma.shape(mdl.glm_NewBMI)$alpha
x_NewBMI <- seq(1, max(hw$Height), by = 2)
xd_NewBMI <- data.frame(Height = x_NewBMI)
mdl.pred_NewBMI <- predict(mdl.glm_NewBMI, xd_NewBMI, type="response")
scale_NewBMI <- mdl.pred_NewBMI / a_NewBMI
mdl.lower_NewBMI <- qgamma(0.025, shape = a_NewBMI, scale=scale_NewBMI)
mdl.upper_NewBMI <- qgamma(0.975, shape = a_NewBMI, scale=scale_NewBMI)
plot(hw$Height, hw$Weight, xlim=c(50, max(hw$Height) * 1.1),
ylim=c(0, max(hw$Weight) * 1.1), xlab="height (cm)", ylab="weight (kg)")
lines(x, mdl.pred_NewBMI, col="red")
lines(x, mdl.lower_NewBMI, col="blue", lty = 2)
lines(x, mdl.upper_NewBMI, col="blue", lty = 2)
コメント 0