はじめに
- Rで回帰分析をした際の出力に似た出力が得られるスクリプト・ファイルを作成してみました(^^)。変なこと書いていたらごめんなさい。(2021/08/20)
- Rのスクリプトに日本語の説明を加えました。(2021/08/29)
- 似たコンセプトのサイトを見つけました。(リンク)(2021/08/31)
Rの回帰分析結果とスクリプトの実行結果
Rの回帰分析結果
> y=c(2,4,3,0,0,1,2); x=c(1,3,4,2,1,0,4)
> m=lm(y~x); summary(m)
Call:
lm(formula = y ~ x)
Residuals:
1 2 3 4 5 6 7
0.9231 1.8077 0.2500 -1.6346 -1.0769 0.4808 -0.7500
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.5192 0.8922 0.582 0.586
x 0.5577 0.3443 1.620 0.166
Residual standard error: 1.327 on 5 degrees of freedom
Multiple R-squared: 0.3441, Adjusted R-squared: 0.2129
F-statistic: 2.623 on 1 and 5 DF, p-value: 0.1662
スクリプトの実行結果
> source("OLS_estimation.R")
Residuals:
1 2 3 4 5 6 7
0.9231 1.8077 0.2500 -1.6346 -1.0769 0.4808 -0.7500
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.5192 0.8922 0.5819 0.5859
x 0.5577 0.3443 1.6196 0.1662
Residual standard error: 1.3272 on 5 degrees of freedom
Multiple R-squared: 0.3441, Adjusted R-squared: 0.2129
F-statistic: 2.6232 on 1 and 5 DF, p-value: 0.1662
スクリプト
ファイル
- OLS_estimation.R(68)
スクリプトの内容
#
# 回帰分析
#
# データの読み込み
y=matrix(c(2,4,3,0,0,1,2)) # y: 7x1 matrix
x=matrix(c(1,3,4,2,1,0,4)) # x: 7x1 matrix
# データ数の取得
n=nrow(y)
# 説明変数行列の作成
# iota=matrix(1,nrow=n,ncol=1)
iota=matrix(rep(1,times=n))
X=cbind(iota,x);
# 説明変数の数の取得
k=ncol(X)
# 回帰係数ベクトルの推定
betahat=solve(crossprod(X),crossprod(X,y))
# 残差ベクトルの計算
uhat=y-X%*%betahat
# 被説明変数の平均からの偏差の二乗和の計算
TSS=as.vector(crossprod(y-mean(y)*iota))
# 残差二乗和の計算
RSS=as.vector(crossprod(uhat))
# 誤差分散の不偏推定値の計算
s2=RSS/(n-k)
# t値からなるベクトルと対応するp値の計算(両側検定)
vcov=s2*solve(crossprod(X))
sevec=matrix(sqrt(diag(vcov)))
tvec=betahat/sevec
pvec=pt(-abs(tvec),df=n-k)+(1-pt(abs(tvec),df=n-k))
# 決定係数と自由度調整済み決定係数の計算
R2=1-RSS/TSS
adjR2=1-(RSS/(n-k))/(TSS/(n-1))
# F値と対応するp値の計算
F=(R2/(k-1))/((1-R2)/(n-k))
pval=1-pf(F,df1=k-1,df2=n-k)
#
# 分析結果出力の準備
#
Residuals=as.vector(uhat)
names(Residuals)=as.character(1:n)
coefficients=cbind(betahat,sevec,tvec,pvec)
dimnames(coefficients)=list(c("(Intercept)","x"),c("Estimate", "Std. Error", "t value", "Pr(>|t|)"))
sigma=sqrt(s2)
r.squared=R2
adj.r.squared=adjR2
fstat=c(F,k-1,n-k,pval)
names(fstat)=c("value","numdf","dendf","p-value")
#
# 分析結果の出力
#
cat("Residuals:\n"); print(round(Residuals,4)); cat("\n")
cat("Coefficients:\n"); print(round(coefficients,4)); cat("\n")
cat(sprintf("Residual standard error: %.4f on %d degrees of freedom\n", sigma, n-k))
cat(sprintf("Multiple R-squared: %.4f,\t", r.squared))
cat(sprintf("Adjusted R-squared: %.4f\n", adj.r.squared))
cat(sprintf("F-statistic: %.4f on %d and %d DF, p-value: %.4f\n", fstat[[1]],k-1,n-k,fstat[[4]]))
せっかくなので
- せっかくなので，かなり以前に作ったものですが，EXCELの回帰分析出力に似た結果を表示するGNU Octaveスクリプトも載せておきます。(2021/08/20)
EXCEL回帰分析結果
スクリプトの実行結果
----------------------------------
回帰統計
----------------------------------
重相関R 0.586607
重決定R2 0.344108
補正R2 0.212930
標準偏差 1.327230
観測数 7
----------------------------------
分散分析表
-------------------------------------------------
自由度 変動 分散 観測された分散比 有意F
-------------------------------------------------
回帰 1 4.620879 4.620879 2.623206 0.166237
残差 5 8.807692 1.761538
合計 6 13.42857
-------------------------------------------------
-----------------------------------------------------------------
係数 標準誤差 t P-値 下限95％ 上限95％
-----------------------------------------------------------------
切片 0.519231 0.892233 0.581945 0.585854 -1.774327 2.812789
X値1 0.557692 0.344333 1.619632 0.166237 -0.327443 1.442828
-----------------------------------------------------------------
スクリプト・ファイル
スクリプト・ファイルの内容
pkg load statistics % 被説明変数ベクトルの作成 y = [2 4 3 0 0 1 2]'; % 標本数の取得 T = length(y); % 説明変数行列の作成 iota = ones(T, 1); x = [1 3 4 2 1 0 4]'; X = [iota x]; % 説明変数の数の取得 k = size(X, 2); % 回帰係数ベクトルのOLS推定値ベクトルの計算 b = (X'*X) \ (X'*y); b1 = b(1); b2 = b(2); % 残差ベクトルの計算 e = y - X*b; % TSS, RSS, ESSの計算 TSS = (y - mean(y)*iota)'*(y - mean(y)*iota); RSS = e'*e; ESS = TSS - RSS; % 分散の不偏推定値の計算 s2 = RSS/(T - k); % 標準誤差の計算 s = sqrt(s2); % OLS推定量の分散共分散行列の計算 cov = s2*inv(X'*X); % 標準誤差の計算 s_b1 = sqrt(cov(1, 1)); s_b2 = sqrt(cov(2, 2)); % t値の計算 t_b1 = b(1) / s_b1; t_b2 = b(2) / s_b2; % P値の計算 p_b1 = tcdf(-abs(t_b1), T - k) + 1 - tcdf(abs(t_b1), T - k); p_b2 = tcdf(-abs(t_b2), T - k) + 1 - tcdf(abs(t_b2), T - k); % 信頼係数0.95の信頼区間の計算 tast = tinv(0.975, T-k); cilb_b1 = b1 - tast*s_b1; ciub_b1 = b1 + tast*s_b1; cilb_b2 = b2 - tast*s_b2; ciub_b2 = b2 + tast*s_b2; % 決定係数，重相関係数，自由度調整済決定係数の計算 R2 = 1 - RSS / TSS; R = sqrt(R2); adjR2 = 1 - (RSS/(T-k)) / (TSS/(T-1)); % F値の計算 F = ( (TSS - RSS)/(k-1) ) / ( RSS/(T-k) ); p_F = 1 - fcdf(F, k-1, T-k); % 結果出力の準備 out_reg = [k-1, ESS, ESS/(k-1), F, p_F]; out_rss = [T-k, RSS, RSS/(T-k)]; out_tol = [T-1, TSS]; out_b1 = [b1, s_b1, t_b1, p_b1, cilb_b1, ciub_b1]; out_b2 = [b2, s_b2, t_b2, p_b2, cilb_b2, ciub_b2]; % -------------------- EXCEL 分析ツール風の出力 ---------------------- % filename = 'octave_excel_reg.txt'; fid = fopen (filename, 'w'); fprintf(fid,'\n') fprintf(fid,'----------------------------------\n') fprintf(fid,' 回帰統計\n') fprintf(fid,'----------------------------------\n') fprintf(fid,'重相関R\t\t %8.6f \n', R) fprintf(fid,'重決定R2\t %8.6f \n', R2) fprintf(fid,'補正R2\t\t %8.6f \n', adjR2) fprintf(fid,'標準偏差\t %8.6f \n', s) fprintf(fid,'観測数\t\t %d \n', T) fprintf(fid,'----------------------------------\n') fprintf(fid,'\n') fprintf(fid,'分散分析表\n') fprintf(fid,'-------------------------------------------------\n') fprintf(fid,' 自由度 変動 分散 観測された分散比 有意F\n') fprintf(fid,'-------------------------------------------------\n') fprintf(fid,'回帰\t %d %9.6f %9.6f %9.6f %9.6f\n', out_reg) fprintf(fid,'残差\t %d %9.6f %9.6f\n', out_rss) fprintf(fid,'合計\t %d %9.5f \n', out_tol) fprintf(fid,'-------------------------------------------------\n') fprintf(fid,'\n') fprintf(fid, '-----------------------------------------------------------------\n') fprintf(fid,' 係数 標準誤差 t P-値 下限95％ 上限95％\n') fprintf(fid, '-----------------------------------------------------------------\n') fprintf(fid,'切片\t %9.6f %9.6f %9.6f %9.6f %9.6f %9.6f \n', out_b1) fprintf(fid,'X値1\t %9.6f %9.6f %9.6f %9.6f %9.6f %9.6f \n', out_b2) fprintf(fid, '-----------------------------------------------------------------\n') fprintf(fid,'\n') fclose (fid);
OLS_estimation.R(68)