Rの例としてよく使用されるirisデータを用いて、残差を簡単に求める。
特に今回は例として、Sepal.Length(がくの長さ)とSepal.Width(がくの幅)について、残差を求める
head(iris) ## Sepal.Length Sepal.Width Petal.Length Petal.Width Species ## 1 5.1 3.5 1.4 0.2 setosa ## 2 4.9 3.0 1.4 0.2 setosa ## 3 4.7 3.2 1.3 0.2 setosa ## 4 4.6 3.1 1.5 0.2 setosa ## 5 5.0 3.6 1.4 0.2 setosa ## 6 5.4 3.9 1.7 0.4 setosa
resultに線形回帰の結果を格納する
result <- lm(Sepal.Length ~ Sepal.Width, data=iris)
例えば、タンパク質の摂取量と総エネルギーの一次回帰式を求めたい場合はlm(prot ~ energy, data=***)とすれば良い
残差はresiduals()で求める。(今回はirisの後ろに新しい「resid1」という列を作成する)
iris$resid1 <- residuals(result) head(iris[,c(1,2,6)]) ## Sepal.Length Sepal.Width resid1 ## 1 5.1 3.5 -0.6444588 ## 2 4.9 3.0 -0.9561394 ## 3 4.7 3.2 -1.1114672 ## 4 4.6 3.1 -1.2338033 ## 5 5.0 3.6 -0.7221227 ## 6 5.4 3.9 -0.2551144
ここでがくの幅が平均値をとるときのがくの長さの期待値を求める
iris$mean1 <- coefficients(result)[1] + coefficients(result)[2]*mean(iris$Sepal.Width)
最後に、resid1とmean1を足し合わせると、残差補正されたSepal.Lengthになる
iris$new1 <- iris$mean1 + iris$resid1 head(iris[,c(1,2,6,7,8)]) ## Sepal.Length Sepal.Width resid1 mean1 new1 ## 1 5.1 3.5 -0.6444588 5.843333 5.198874 ## 2 4.9 3.0 -0.9561394 5.843333 4.887194 ## 3 4.7 3.2 -1.1114672 5.843333 4.731866 ## 4 4.6 3.1 -1.2338033 5.843333 4.609530 ## 5 5.0 3.6 -0.7221227 5.843333 5.121211 ## 6 5.4 3.9 -0.2551144 5.843333 5.588219
最後に作成したnew1が残差補正されたがくの長さということになる。
これは栄養疫学でよく用いられるWilletらによって提唱された残差法(総エネルギー摂取量を補正する方法)の一例である。
最後に複数の項目について、残差補正するような場合を考えたサンプルを紹介します
これはirisの一列目を基準にして、2~4列目までの項目で残差補正した値を求めるコードです。
for(i in 2:4) { y <- iris[,i] x <- iris[,1] result <- lm(y ~ x) resid <- residuals(result) mean <- coefficients(result)[1] + coefficients(result)[2]*mean(x) new <- resid + mean iris <- cbind(iris, new) }
以上
20160523
RF
20210330 改定済(一部表記ミスがありました)