統計学と疫学と時々、助教生活

疫学を専門とする助教の研究に関する備忘録的ページ。

Rで残差補正された値を求める

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 改定済(一部表記ミスがありました)