2016/07/27

2020/05/04

R言語で線形モデルによる回帰分析の自作関数例

Rで線形モデルによる回帰分析を行う自作関数を作成ました。詳細な結果をcsv形式で出力するまで行なってくれます。

当ページは自作関数の紹介のみになっています。Rによる回帰分析の詳細が知りたい方は、入門コーナーの記事であるR言語で線形モデルによる回帰分析を参照していただければと思います。Pythonで回帰分析をやりたい方はこちら→Python3で線形モデルによる回帰分析とプロット

↓これが今回、解析に使うデータです。

sample-data

今回紹介する関数は、「単（重）回帰分析を行い、必要なデータを全て取り出し、まとめて、その結果をエクセルファイルに出力するというような関数」

となっております。

※R言語入門のトップページはこちら

1 自作関数”reg_fun()”
2 ”reg_fun()”で単回帰分析をしてみる
3 for文によって一気に全ての列に関して単回帰分析
4 決定係数も出力する回帰分析関数
5 今回使ったコマンドや関数

自作関数”reg_fun()”

以下がその関数になります。引数には

解析用データ、出力ファイルの名前、目的変数の名前

を入れてあげてください。

プログラム

reg_fun <- function(dat,filename,response)
{
ans <- lm(dat$Y~.,data=dat)
s.ans=summary(ans)
coe <- s.ans$coef
CI.low <- coe[,1]-1.96*coe[,2]
CI.up <- coe[,1]+1.96*coe[,2]
aic <- AIC(ans)
N <- nrow(dat)
result <- cbind(coe,CI.low,CI.up,aic,N)
result[2:nrow(result),c(7,8)]=””
res <- paste(“Y = “,response)

write.table(res,filename,append=T,quote=F,sep=”,”,row.names=F,col.names=F)
write.table(matrix(c(“”,colnames(result)),nrow=1),filename,append=T,quote=F,sep=”,”,row.names=F,col.names=F)
write.table(result,filename,append=T,quote=F,sep=”,”,row.names=T,col.names=F)
write.table(“”,filename,append=T,quote=F,sep=”,”,row.names=F,col.names=F)
}

ここで目的変数の列名を”Y”と指定していることに気をつけましょう。関数を使う前に、解析用データの目的変数の列名を”Y”に書き変える必要があります。

”reg_fun()”で単回帰分析をしてみる

今作ったreg_fun()をRに読み込んで使えるようにします。そして、実際にこれを使って回帰分析を行ってみましょう。

解析からcsvファイル出力まで関数にして解析が楽になったので、せっかくだから、病気を目的変数にとり、それ以外の全ての列に関して単回帰分析を行ってみましょう。

まずは出力ファイル名を決めます。

プログラム

filename <- “全列単回帰分析_目的変数は病気.csv”

今回はこのようにします。

for文によって一気に全ての列に関して単回帰分析

以下ようにfor分を使う事によって、全ての列に関しての単回帰分析を一気に行うことが出来ます。

プログラム

for(i in 1:ncol(df))
{
if(i != 5)
{
dat <- df[,c(5,i)]
colnames(dat)[1] = “Y”
reg_fun(dat,filename,”病気”)
}
}

ncol(df)はdfの列数を表していて、for分をdfの列数分繰り返すという意味です。if分は説明変数と目的変数が両方病気になることを防いでいます。そして、目的変数の列名を”Y”に変更するのも忘れずに行いましょう。

結果はこのようになります。

これだけの単回帰分析をこれだけのプログラムで出来るのはなかなか便利ですよね。今回の結果から、有意水準0.05ならば目的変数が体重の場合のみ有意な結果になるということが分かります。

決定係数も出力する回帰分析関数

決定係数と調整済み決定係数もcsvファイルに出力するバージョンです。是非お役立て下さい。

プログラム

reg_fun <- function(dat,filename,response)

{

ans <- lm(dat$Y~.,data=dat)

s.ans=summary(ans)

coe <- s.ans$coef

CI.low <- coe[,1]-1.96*coe[,2]

CI.up <- coe[,1]+1.96*coe[,2]

aic <- AIC(ans)

R.squared = s.ans$r.squared

Adj.R.squared = s.ans$adj.r.squared

N <- nrow(dat)

result <- cbind(coe,CI.low,CI.up,aic,R.squared,Adj.R.squared,N)

result[2:nrow(result),c(7:10)]=””

res <- paste(“Y = “,response)

write.table(res,filename,append=T,quote=F,sep=”,”,row.names=F,col.names=F,fileEncoding=”CP932″)

write.table(matrix(c(“”,colnames(result)),nrow=1),filename,append=T,quote=F,sep=”,”,row.names=F,col.names=F,fileEncoding=”CP932″)

write.table(result,filename,append=T,quote=F,sep=”,”,row.names=T,col.names=F,fileEncoding=”CP932″)

write.table(“”,filename,append=T,quote=F,sep=”,”,row.names=F,col.names=F,fileEncoding=”CP932″)

}

今回使ったコマンドや関数

今回使ったRのコマンドや関数を一覧にしてまとめました。

コマンド・関数	使い方・機能
lm()	線形モデルを用いた回帰分析。lm(目的変数のベクトル~.,解析データ)で、解析を実行。単回帰、重回帰共に実行できます。
summary()	回帰分析の結果を引数にすることで、より詳細な結果を表示する。 ※行列やベクトルにも使用可能。
$coefficient	今回の回帰分析の場合、summary()によって詳細を抽出した回帰分析の結果に対して、回帰係数,標準誤差,t値,p値を抽出する。
exp(coe[,1]-1.96coe[,2]) exp(coe[,1]+1.96coe[,2])	回帰分析における、95％信頼区間を求めている。
AIC	回帰分析のAICを抽出。
plot()	plot(x,y)でベクトルx,yについて2次元のグラフ上に散布図をプロット。
reg_fun()	当サイトによる自作関数。引数に解析データ、出力ファイル名、目的変数名を取る。線形モデルによる単（重）回帰分析を行い、回帰係数、標準誤差、t値、p値、RR、95%信頼区間、AIC、データの要素数をcsvファイルに出力。