はじめに
Rによる重回帰分析を方法と、プロットを出力するまでの流れをまとめます。
自分が調べたかぎり、2020/11時点では重回帰分析の結果をプロットしてくれる機能はなさそうだったので、力ずくで行いました。
1 重回帰分析とは
ある「目的変数」を複数の「説明変数」で表すために行います。
身近な例として、目的変数を「体育の評定」とします。
仮に、50m走が速い学生は体育の成績が高い傾向があるとします。
しかし50m走が遅くても握力が大きい学生、1500mが速い学生も存在しており、そうした学生も高い評定を得る可能性があります。
このような場合に"50m走"、"1500m走"、"握力"、"反復横跳び"、"立ち幅跳び"、"前屈"、"ハンドボール投げ"といった複数の項目から評価することで、より精度の高い予測を行うことができます。
今回は"遺伝子発現"を目的変数、"外部環境要因"を説明変数とした重回帰分析を行いました。
2 重回帰分析
setwd("/Volumes/~/~") #作業ディレクトリの指定
data <- read.csv("Gene_Environment.csv", sep = ",", header = TRUE) #用いるデータの読み込み
z <- scale(data) #scaleで標準化（平均0、標準偏差1に）
z <- data.frame(z) #データフレーム形式に変換
m <- lm(MEpurple~., z) #lm関数で重回帰
summary(m) #結果を表示
m2<-step(m) #AICをもとに説明変数の選択を自動で行ってくれる
summary(m2)#変数選択後の結果を表示
sink("multi_regression.txt") #summaryをテキスト形式で作業中のディレクトリに保存 #" "の文字を変えることでつけたい名前をつけられる
print(summary(lm(m2))) #出力内容をR上で確認
sink() #出力を解除
3 プロット表示の考え方
得られた説明変数の係数に対して各説明変数を乗じたものの和が、各サンプルにおける予測値である。
全てのサンプルにおいて予測値を算出し、予測値-目的変数の散布図を描くことが目的。
a<- c() #空のベクトルを用意
for (i in 2:11){
x <- 説明変数係数[i,1]*サンプル固有値[1,i]
a<- c(a, c(x))
}
a
sum(a)
4 重回帰予測値プロット作成
write.csv(z, file = "standarized_gene_env.csv") #標準化されたデータセットを用意、書き出しの必要は特にない
co <- data.frame(m2$coefficients) #重回帰分析した結果の中から係数のデータセットを抽出し,"co"に格納
ad <- c() #空のベクトルを用意
for (j in 1:100){ #今回用いたものは、サンプル1~100まで存在するため1:100
for (i in 2:11){ #説明変数が10種類存在した #1行目には説明変数、2~11行目に目的変数があるため2:11
x <- co[i,1]*z[j,i] #係数*サンプル固有値を全サンプルでループして行う。(sample1の目的変数10項目を算出 -> sample2の目的変数10 -> sample3.... -> sample100)
ad<- c(ad, c(x)) #計算過程をベクトルに格納
}
}
ad #結果表示、この時点で計算ミスないか検算したがいい
ad2 <- matrix(ad, nrow=100, ncol=10, byrow=T) #現在は1列のベクトルデータになっているため100行×10列の行列に変換、byrowにより左から右へ格納していく（今回は1行目10列分を埋めたら、２行目に移る…の繰り返し）
ad3 <- rowSums(ad2) #行の和を算出 (10個の値の総和)
ad3 <- data.frame(ad3) #データフレームにする
ad4 <- ad3 + (co[1,1]) #切片のデータを全てのデータに足す（これにより各サンプルの予測値の準備完了）
ad5 <- cbind(z[,1], ad4 ) #ModuleEigengeneと予測値データをまとめた行列の作成
colnames(ad5) <- c("eigengene", "Prediction") #各列のデータラベルを指定
plot(Module_eigengene ~ Prediction, data = ad5,pch = 21,
bg = "pink",main= "Multiple regression analysis of eigengenes with factors",
xlab = "Prediction by multiple regression", ylab = "eigengenes",bty="L" ,cex = 1.5)
result <- lm(Module_eigengene ~ Prediction, data = ad5) # 回帰分析を行う
abline(result,lty = "dotted", lwd = 2) # 推定回帰直線を描く
summary(result)
もし重回帰分析後のデータからすぐにプロットを出力できる関数、機能等ありましたら教えて頂けるとありがたいです。