搁による濒尘関数を使わない回帰と梃子比?肠辞辞办の距离
?Download as DOCX, PDF?
0 likes?2,092 views
ちょっとマニアックですが、行列計算で回帰係数を求めてt検定してみました。 その過程で梃子比を求め、最後にcookの距離とどう違うのかを図示しています。
![2011 年 5 月作成
lm 関数を使わない回帰分析と梃子比?Cook の距離
○ テストデータ作成
誤差項に自由度 10 の t 分布乱数を用いた回帰モデルに従うデータを作成。n=1000。
rm(list=ls(all=TRUE)) # 作業領域のクリア
setwd("d:/TSR2/sim") # データファイルなどを置く作業ディレクトリ
の指定
library(sn) # for rmst （t 分布に従う乱数を発生させる
関数）
set.seed(24) # 乱数を固定するためのシード設定
a1 <- 2 # 回帰パラメータの傾きを指定
b1 <- 5 # 回帰パラメータ切片
x <- runif(1000, min=1, max=10) # 説明変数は一様乱数
y <- a1 * x + b1 + rt(1000, df= 10)
# 被説明変数は回帰モデルに従い自由度 10 の t 分布乱数残差を持つ
○ lm 関数を使わない回帰分析と梃子比
x の先頭に要素が全て 1 の列を加えて行列計算を行う。 逆行列が計算できれば回帰分析
は可能。hat 行列（ここでは hat1）の対角成分が梃子比となるが、hat 行列は説明変数の
みから算出されており、梃子比は単に説明変数の平均からの近さで決まる。
x1 <- cbind(rep(1, 1000), x) # x の先頭に全て 1 の列を加える
hat1 <- x1 %*% solve((t(x1) %*% x1)) %*% t(x1) # hat 行列の計算
teko<- diag(hat1)
(bt1 <- solve(t(x1) %*% x1) %*% t(x1) %*% y) # 回帰係数算出
res1 <- y - hat1 %*% y # 残差
summary(res1) # 残差の四分位値を含む基本統計量算出
s2 <- sum(res1^2) / (nrow(x1) - ncol(x1)) # 誤差項の不偏分散推定量
sqrt(s2) # 誤差項の標準偏差
(se1 <- (diag(s2*solve(t(x1)%*%x1)))^0.5) # 回帰??????の標準誤差
##### t 検定
(t1 <- bt1 / se1) # t値
2 * pt(-abs(t1), nrow(x1) - ncol(x1)) # 両側 t 検定の p 値
(R2 <- 1- sum(res1^2)/(sum((y-mean(y))^2))) # 決定係数
(R2r<- 1- sum(res1^2) / (sum((y-mean(y))^2))*(nrow(x1)-1)
/ (nrow(x1)-ncol(x1))) # 自由度調整済み決定係数
##### F 検定 (切片以外がゼロであるを検定する F 統計値の計算)
(f1 <-R2/(1-R2)*(nrow(x1)-ncol(x1))/(ncol(x1)-1))
1-pf(f1,ncol(x1)-1,nrow(x1)-ncol(x1)) # F 検定の p 値の計算
summary(lm(y~x)) # lm 関数の結果と比較してみる
flg<- rep(1, 1000)
flg[which(teko>boxplot(teko)$stats[3])] <- 2 # 中位数
flg[which(teko>boxplot(teko)$stats[4])] <- 3 # 四分位値
flg[which(teko== boxplot(teko)$stats[5])] <- 4 # 最大値のデータ
plot(x, y, col=c("black", "yellow", "orange", "red")[flg], pch=20)
abline(lm1, col="blue", lwd=2) # 回帰線入り](https://image.slidesharecdn.com/lmcookslideshare-121214093406-phpapp02/85/R-lm-cook-1-320.jpg)
![○ Cook の距離
Cook の距離は全てのデータを用いた際と、 １つのデータを除いた場合の回帰式による予
測値の差に関する距離尺度。 Cook の距離が大きいデータは、回帰式に与える影響が大きく、
異常値である可能性がある。
一般に、Cook の距離が 0.5 以上であれば大きいとされる。
R では関数 cooks.distance を用いて算出することができる。
cd1 <- cooks.distance(lm1) # Cook の距離算出
flg2<- rep(1, 1000) # 色分け用のフラグ
flg2[which(cd1 >= 0.0003)] <- 2 # 適当に大きな数字に色をつける
flg2[which(cd1 >= 0.001)] <- 3
plot(x, y, col=c("black","orange","red")[flg2], pch=20)
abline(lm1, col="blue", lwd=2) # 回帰線入り
25
25
20
20
y
y
15
15
10
10
5
5
2 4 6 8 10 2 4 6 8 10
x x
梃子比 Cook の距離](https://image.slidesharecdn.com/lmcookslideshare-121214093406-phpapp02/85/R-lm-cook-2-320.jpg)
搁による濒尘関数を使わない回帰と梃子比?肠辞辞办の距离
- 1. 2011 年 5 月作成 lm 関数を使わない回帰分析と梃子比?Cook の距離 ○ テストデータ作成 誤差項に自由度 10 の t 分布乱数を用いた回帰モデルに従うデータを作成。n=1000。 rm(list=ls(all=TRUE)) # 作業領域のクリア setwd("d:/TSR2/sim") # データファイルなどを置く作業ディレクトリ の指定 library(sn) # for rmst （t 分布に従う乱数を発生させる 関数） set.seed(24) # 乱数を固定するためのシード設定 a1 <- 2 # 回帰パラメータの傾きを指定 b1 <- 5 # 回帰パラメータ切片 x <- runif(1000, min=1, max=10) # 説明変数は一様乱数 y <- a1 * x + b1 + rt(1000, df= 10) # 被説明変数は回帰モデルに従い自由度 10 の t 分布乱数残差を持つ ○ lm 関数を使わない回帰分析と梃子比 x の先頭に要素が全て 1 の列を加えて行列計算を行う。 逆行列が計算できれば回帰分析 は可能。hat 行列（ここでは hat1）の対角成分が梃子比となるが、hat 行列は説明変数の みから算出されており、梃子比は単に説明変数の平均からの近さで決まる。 x1 <- cbind(rep(1, 1000), x) # x の先頭に全て 1 の列を加える hat1 <- x1 %*% solve((t(x1) %*% x1)) %*% t(x1) # hat 行列の計算 teko<- diag(hat1) (bt1 <- solve(t(x1) %*% x1) %*% t(x1) %*% y) # 回帰係数算出 res1 <- y - hat1 %*% y # 残差 summary(res1) # 残差の四分位値を含む基本統計量算出 s2 <- sum(res1^2) / (nrow(x1) - ncol(x1)) # 誤差項の不偏分散推定量 sqrt(s2) # 誤差項の標準偏差 (se1 <- (diag(s2*solve(t(x1)%*%x1)))^0.5) # 回帰??????の標準誤差 ##### t 検定 (t1 <- bt1 / se1) # t値 2 * pt(-abs(t1), nrow(x1) - ncol(x1)) # 両側 t 検定の p 値 (R2 <- 1- sum(res1^2)/(sum((y-mean(y))^2))) # 決定係数 (R2r<- 1- sum(res1^2) / (sum((y-mean(y))^2))*(nrow(x1)-1) / (nrow(x1)-ncol(x1))) # 自由度調整済み決定係数 ##### F 検定 (切片以外がゼロであるを検定する F 統計値の計算) (f1 <-R2/(1-R2)*(nrow(x1)-ncol(x1))/(ncol(x1)-1)) 1-pf(f1,ncol(x1)-1,nrow(x1)-ncol(x1)) # F 検定の p 値の計算 summary(lm(y~x)) # lm 関数の結果と比較してみる flg<- rep(1, 1000) flg[which(teko>boxplot(teko)$stats[3])] <- 2 # 中位数 flg[which(teko>boxplot(teko)$stats[4])] <- 3 # 四分位値 flg[which(teko== boxplot(teko)$stats[5])] <- 4 # 最大値のデータ plot(x, y, col=c("black", "yellow", "orange", "red")[flg], pch=20) abline(lm1, col="blue", lwd=2) # 回帰線入り
- 2. ○ Cook の距離 Cook の距離は全てのデータを用いた際と、 １つのデータを除いた場合の回帰式による予 測値の差に関する距離尺度。 Cook の距離が大きいデータは、回帰式に与える影響が大きく、 異常値である可能性がある。 一般に、Cook の距離が 0.5 以上であれば大きいとされる。 R では関数 cooks.distance を用いて算出することができる。 cd1 <- cooks.distance(lm1) # Cook の距離算出 flg2<- rep(1, 1000) # 色分け用のフラグ flg2[which(cd1 >= 0.0003)] <- 2 # 適当に大きな数字に色をつける flg2[which(cd1 >= 0.001)] <- 3 plot(x, y, col=c("black","orange","red")[flg2], pch=20) abline(lm1, col="blue", lwd=2) # 回帰線入り 25 25 20 20 y y 15 15 10 10 5 5 2 4 6 8 10 2 4 6 8 10 x x 梃子比 Cook の距離