ずいぶんと間が空きました。学会の季節になると統計分析のご依頼やご相談や査読者への対応など立て込んで参ります。
当研究所は2000年にインターネット上で開設しましたが、当時、気軽に統計分析を利用できるWebサイトは少なかったと思います。お陰様で沢山の経験をさせて頂きました。分析費用が1桁違うのでは・・、と言われたりします。しかし、古い言葉かも知れませんが"苦学生"などの懐具合を思えば商売下手でいいと思ってやっています。
当研究所は2000年にインターネット上で開設しましたが、当時、気軽に統計分析を利用できるWebサイトは少なかったと思います。お陰様で沢山の経験をさせて頂きました。分析費用が1桁違うのでは・・、と言われたりします。しかし、古い言葉かも知れませんが"苦学生"などの懐具合を思えば商売下手でいいと思ってやっています。
今回は、重回帰分析についてですが、通常(最小二乗法)の回帰分析ではなく「ロバスト回帰分析」について、その手法のほんの一例をご紹介したいと思います。
「すぐに役立つ統計のコツ」(オーム社)の94ページでは、普通(最小二乗法)の回帰分析しかご紹介していません。
そこで、
ガンマグロブリン(x)と反応する蛋白化学反応(y)の関係を例題にして、データ解析環境「R」での手法をご紹介します。
そこで、
ガンマグロブリン(x)と反応する蛋白化学反応(y)の関係を例題にして、データ解析環境「R」での手法をご紹介します。
「R」の環境において、
新しいスクリプトを立ち上げ、次のプログラムを「コピー&ペースト」して下さい。
新しいスクリプトを立ち上げ、次のプログラムを「コピー&ペースト」して下さい。
***
Y<- c(4.3,6.2,16.2,11,7,17.1,8.5,16,29.5,12.1,10.6,8.9,46.7,28.1,3.9,5.4,20.1,20.8,21.7,45.0)
X<- c(0.77,1.04,2.46,1.48,0.92,2.03,0.97,2.02,2.57,1.53,1.48,1.27,6.4,2.67,0.83,0.6,2.7,2.98,2.6,1.56)
Y<- c(4.3,6.2,16.2,11,7,17.1,8.5,16,29.5,12.1,10.6,8.9,46.7,28.1,3.9,5.4,20.1,20.8,21.7,45.0)
X<- c(0.77,1.04,2.46,1.48,0.92,2.03,0.97,2.02,2.57,1.53,1.48,1.27,6.4,2.67,0.83,0.6,2.7,2.98,2.6,1.56)
dat<- data.frame(y=Y, x=X)
dat # 原データの確認
***
dat # 原データの確認
***
外れ値(x=1.56, y=45.0)を除き"dat1"とします。
***
dat1<- dat[-20, ] # 原データから20行のデータを削除
dat1 # 外れ値を除いたデータの確認
***
***
dat1<- dat[-20, ] # 原データから20行のデータを削除
dat1 # 外れ値を除いたデータの確認
***
原データ(dat)と外れ値削除後データ(dat1)の回帰直線を見てみましょう。
***
windows(width=4, height=6)
plot(dat$x, dat$y, pch=20, col="red", xlab="dat$x", ylab="dat$y") # 外れ値(赤色)
par(new=T)
plot(dat1$x, dat1$y, pch=20, xlab="", ylab="")
abline(lm(dat$y~ dat$x), lty=1, col="red")
abline(lm(dat1$y~ dat1$x), lty=1) # 外れ値の削除
***
***
windows(width=4, height=6)
plot(dat$x, dat$y, pch=20, col="red", xlab="dat$x", ylab="dat$y") # 外れ値(赤色)
par(new=T)
plot(dat1$x, dat1$y, pch=20, xlab="", ylab="")
abline(lm(dat$y~ dat$x), lty=1, col="red")
abline(lm(dat1$y~ dat1$x), lty=1) # 外れ値の削除
***
図1 原データ(dat)と削除後データ(dat1)の回帰直線
図1の赤線は原データの、黒線は削除後の回帰直線です。
表1 削除後データ(dat1)によるlm法(普通の方法)の結果
----------------------------------------------------------
fit.lm<- lm(dat1$y~ dat1$x)
summary(fit.lm)
fit.lm<- lm(dat1$y~ dat1$x)
summary(fit.lm)
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.2049 1.3947 0.147 0.885
dat1$x 7.7762 0.5941 13.089 2.63e-10
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.2049 1.3947 0.147 0.885
dat1$x 7.7762 0.5941 13.089 2.63e-10
----------------------------------------------------------