################################################ #南風原(2002)第8章 偏相関と重回帰分析 ################################################ #Step1:データと基本統計量(表8-1,図8-2,表8-2)p.226-27 ################################################ ##(a)50組の母子の協調性データを読み込み setwd("C:/data/haebara/chap08") data1 <- read.csv("http://blue.zero.jp/yokumura/R/haebara/chap08.csv",header=T) data1 #確認 ##(b)「番号」という変数は必要ないので, #使用する変数だけsubset()で抽出する data1 <- subset(data1,select=c(母親価値,通園年数,協調性)) data1 #確認(表8-2) ##(c)3変数の散布図と回帰平面(図8-2)を描く #3次元散布図を描くには,scatterplot3dライブラリを使う library(scatterplot3d) #ライブラリを呼び出し attach(data1) #使用するデータを変数名で指定できるようにする #--------------------------------------# #lm(従属変数~独立変数+独立変数...+独立変数)で重回帰分析を実行する #そして,重回帰分析の結果を, #reg0(名前は何でもいい)という新しいオブジェクトに保存する #--------------------------------------# reg0 <- lm(協調性~ 母親価値 + 通園年数) #--------------------------------------# #scatterplot3d(x,y,z,angle,main)で3次元散布図を描く #angleの数字を適当に変えることにより,角度が変化する #最初は,angleを指定せずに実行して, #次に,適当な数字を代入して行き,気に入った図が描ければOK #なお,scatterplot3d()の結果を, #s2d(名前は何でもいい)という新しいオブジェクトに保存している #--------------------------------------# s2d <- scatterplot3d(x=母親価値,y=通園年数,z=協調性, angle=20,main="3変数の散布図と回帰平面") #--------------------------------------# #「scatterplot3d()の結果を保存したオブジェクト名$plane3d()」で回帰平面を描く #ここでは,s2d$plane3d(重回帰分析の結果, col="red") #と指定することで,3次元散布図に,回帰平面を上書きしている #--------------------------------------# s2d$plane3d(reg0,col="red") graphics.off() #図を閉じる ##(d)協調性データの基本統計量(表8-2)を求める apply(data1,2,mean) #平均 apply(data1,2,sd) #不偏分散を基にした標準偏差 #標準偏差関数を定義(ssd)と命名 ssd <- function(x){ #BeginFnc sqrt(1/length(x)*sum((x-mean(x))^2)) } #EndFnc apply(data1,2,ssd) #標本分散を基にした標準偏差 cor(data1) #相関係数 ################################################ #Step2:残差への注目 p.227-29 ################################################ ##(a)残差x2|x1と子供の協調性得点との散布図(図8-3) #--------------------------------------# #母親価値(x1)によって,通園年数(x2)を予測する回帰分析を行う #回帰分析は,lm(従属変数~独立変数)で実行できる #ここでは,回帰分析の結果を,reg1というオブジェクトに保存している #--------------------------------------# reg1 <- lm(通園年数~母親価値) #--------------------------------------# #残差(通園年数から母親価値の影響を除いた成分)は, #(1) 回帰分析の結果$residuals #(2) residuals(回帰分析の結果) #で取り出せる(第3章 参照) #--------------------------------------# reg1$residuals #残差1 residuals(reg1) #残差2 cor(母親価値,reg1$residuals) #残差は母親価値(x1)とは無相関. #-3.289373e-17と表示されるが, #-3.289373 * 10^(-17)を意味する #要するに,実質「0」である. #--------------------------------------# #plot(y軸~x軸)で図8-3を再現 #plotした後,回帰直線を描くために #lm()で回帰分析を行い,その結果を保存する(reg2). #そして,abline(回帰分析の結果)で回帰直線が描ける #--------------------------------------# plot(協調性~reg1$residuals) reg2 <- lm(協調性~reg1$residuals) abline(reg2) cor(協調性,reg1$residuals) ##部分相関 cor(y, x2|x1) graphics.off() #図を閉じる ################################################ #Step3:部分相関係数 p.229-30 ################################################ ##(a)部分相関係数の公式(8.1式)を定義 #x2からx1の影響を除いたときのyとの部分相関係数 r_y(2|1) part.cor <- function(x1,x2,y){ (cor(y,x2)-cor(y,x1)*cor(x1,x2))/ #分子 sqrt(1-cor(x1,x2)^2) #分母 } part.cor(x1=母親価値,x2=通園年数,y=協調性) #確認 cor(協調性,reg1$residuals) #回帰分析を利用した結果と同一 ################################################ #Step4:偏相関係数 p.230-31 ################################################ ##(a)残差x2|x1と残差y|x1との散布図 reg3 <- lm(協調性~母親価値) reg3$residuals #協調性から母親価値の影響を除いた成分 reg1$residuals #通園年数から母親価値の影響を除いた成分 plot(reg3$residuals~reg1$residuals) reg4 <- lm(reg3$residuals~reg1$residuals) abline(reg4) #図8-4 graphics.off() #図を閉じる ##(b)図8-3と図8-4を比べる par(mfrow=c(1,2)) #1行2列の形で図を描く plot(協調性~reg1$residuals) abline(reg2) #図8-3 plot(reg3$residuals~reg1$residuals) abline(reg4) #図8-4 graphics.off() #図を閉じる ##(c)偏相関の公式(8.2式) partial.cor <- function(x1,x2,y){ (cor(y,x2)- cor(y,x1) * cor(x1, x2))/ #分子 (sqrt(1-cor(y,x1)^2) * sqrt(1-cor(x1,x2)^2)) #分母 } partial.cor(x1=母親価値,x2=通園年数,y=協調性) #確認 cor(reg3$residuals,reg1$residuals) #偏相関(回帰分析の結果を利用) ##(d) 偏相関と部分相関の関係(8-3式) part.cor(x1=母親価値,x2=通園年数,y=協調性)/sqrt(1-cor(協調性,母親価値)^2) ##(e)ggmライブラリを利用して偏相関を求める library(ggm) #ライブラリを呼び出し #--------------------------------------# #correlations(data) #で下側三角に相関係数, #上側三角に「残りの全変数の影響を除いた場合の偏相関係数」が表示される #> correlations(data1) # 母親価値 通園年数 協調性 #母親価値 1.0000000 0.3913046 0.2882904 #通園年数 0.5376770 1.0000000 [0.3290749]------>> これが,教科書で求めている偏相関係数 #協調性 0.4799263 0.5014633 1.0000000 #--------------------------------------# correlations(data1) #--------------------------------------# #pcor(u=c(使用する変数名), S=cov(data)) #「使用する変数名」で指定した, #1番目と2番目の偏相関係数が求まる. #3番目以降は,影響を除去する変数名を指定する #ここでは,母親価値の影響を除いた場合の,協調性と通園年数の偏相関を求める #--------------------------------------# pcor(u=c("協調性","通園年数","母親価値"), S=cov(data1)) ################################################ #Step5:偏回帰係数 p.232-234 ################################################ ##(a)図8-3と図8-4の散布図(再掲) 傾きが同じことが確認できる par(mfrow=c(1,2)) #1行2列の形で図を描く plot(協調性~reg1$residuals) abline(reg2) #図8-3 plot(reg3$residuals~reg1$residuals) abline(reg4) #図8-4 graphics.off() #図を閉じる ##(b)回帰係数が等しいことを数字で確認する reg2 reg4 #回帰係数は等しい ##(c) e+?の意味を確認する 1.126e+2 #1.126*(10^2)=112.6 1.126e+1 #1.126*10^1=11.26 1.126e+0 #1.126*10^0=1.126 1.624e-16 #1.624*10^(-16)=0.0000000000000001624 ##(d) 8.5式 part.cor(母親価値,通園年数,協調性) * ssd(協調性) / (ssd(通園年数)*sqrt(1-cor(母親価値,通園年数)^2)) partial.cor(x1=母親価値,x2=通園年数,y=協調性) * (ssd(協調性)*sqrt(1-cor(協調性,母親価値)^2))/(ssd(通園年数)*sqrt(1-cor(母親価値,通園年数)^2)) ################################################ #Step6:標準偏回帰係数 p.234-35 ################################################ ##(a) 8.7式 (cor(協調性,通園年数) - cor(協調性,母親価値) * cor(母親価値,通園年数) )/ (1-cor(母親価値,通園年数)^2) ##(b) 標準偏回帰係数を求めるための関数を定義(覚えなくても大丈夫) #--------------------------------------# #model <- lm(従属変数~独立変数+...) #std.est(model) #で標準偏回帰係数が算出される #ここでは,一番最初に定義した #reg0 <- lm(協調性~ 母親価値 + 通園年数) #を利用して, #std.est(reg0) #とすることにより,標準偏回帰係数を求めている #--------------------------------------# std.est <- function(model){ sdd <- c(0,sd(model$model[-1])) std.est <- coef(model)*sdd/sd(model$model[1])[1] print(std.est) } std.est(reg0) ################################################ #Step7:最小二乗法による母数の推定 p.238-40 ################################################ ##(a) 8.15式 b1 <- ((cor(協調性,母親価値) - cor(協調性,通園年数)*cor(母親価値,通園年数))/(1-cor(母親価値,通園年数)^2))* ssd(協調性)/ssd(母親価値) b1 #確認 ##(b) 8.16式 b2 <- ((cor(協調性,通園年数)- cor(協調性,母親価値)*cor(母親価値,通園年数))/(1-cor(母親価値,通園年数)^2))* ssd(協調性)/ssd(通園年数) b2 #確認 ##(c) 8.17式 a <- mean(協調性)- b1*mean(母親価値) - b2*mean(通園年数) a #確認 reg0 #関数を使った場合と値が一致することを確認 ##(d) 8.18式 b1.std <- ((cor(協調性,母親価値) - cor(協調性,通園年数)*cor(母親価値,通園年数))/(1-cor(母親価値,通園年数)^2)) b1.std #確認 ##(e) 8.19式 b2.std <- ((cor(協調性,通園年数)- cor(協調性,母親価値)*cor(母親価値,通園年数))/(1-cor(母親価値,通園年数)^2)) b2.std std.est(reg0) #関数を使った場合と値が一致することを確認 ################################################ #Step8:重相関係数 p.240-41 ################################################ ##(a) 独立変数が2個の場合の重相関係数 8-20式 MultipleR <- sqrt( (cor(協調性,母親価値)^2 + cor(協調性,通園年数)^2 - 2* cor(協調性,母親価値) * cor(協調性,通園年数)*cor(母親価値,通園年数))/ (1-cor(母親価値,通園年数)^2) ) MultipleR #確認 ##(b) summary(reg0) #lm(従属変数~独立変数1 + 独立変数2) #表示されるMultiple R-Squaredは,重相関係数の2乗 MultipleR^2 #確認 ################################################ #Step9:従属変数との相関と重相関係数 p.243-245 ################################################ ##(a) 独立変数と従属変数の相関が.7,独立変数間の相関は.54の場合の重相関係数 sqrt(((0.7^2 + 0.7^2 -2*0.7*0.7*0.54)/(1-0.54^2))) ##(b) 重相関係数の計算式(8.20式)を定義 multi.cor <- function(yx1,yx2,x1x2){ sqrt( (yx1^2 + yx2^2 - 2 * yx1 * yx2 * x1x2 )/ ( 1-x1x2^2 ) ) } multi.cor(yx1=.48,yx2=.5,x1x2=.54) #表8-1のデータ multi.cor(yx1=.7,yx2=.7,x1x2=.54) #Step9 (a) の場合 multi.cor(yx1=.48,yx2=.5,x1x2=.1) #p.244 9行目 multi.cor(yx1=.48,yx2=.5,x1x2=.0) #p.244 16行目 multi.cor(yx1=.48,yx2=.5,x1x2=-.5) #p.244 17行目 ##(c) 図8-7 独立変数間の相関係数の変化に伴う重相関係数の変化を描画 x <- seq(-1,1,by=0.01) #x軸用の値を準備 multi.cor(yx1=.6,yx2=.6,x1x2=x) #ry1 = ry2 = .6の場合 plot(multi.cor(yx1=.6,yx2=.6,x1x2=x)~x, type="l", ylim=c(0,1)) lines(multi.cor(yx1=.4,yx2=.4,x1x2=x)~x) #ry1 = ry2 = .4の場合 lines(multi.cor(yx1=.2,yx2=.2,x1x2=x)~x) #ry1 = ry2 = .2の場合 graphics.off() #図を閉じる ################################################ #Step10:平方和の分割と重相関係数の検定 p.250-55 ################################################ ##(a) 分散説明率 (決定係数)8.29式 summary(reg0) #決定係数 Multiple R-Squared summary(reg0)$r.squared #summary(lm())$r.squaredで決定係数だけを取り出す ##(b) 予測の標準誤差 8.33式 ssd(協調性)*sqrt(1-summary(reg0)$r.squared) ##(c) 予測の標準誤差 8.34式 ssd(協調性)*sqrt(1-summary(reg0)$r.squared)*sqrt(50/(50-2-1)) summary(reg0) #標準誤差はResidual standard errorに表示される summary(reg0)$sigma #summary(lm())$sigmaで標準誤差だけを取り出す ##(d) 自由度調整済み決定係数 8.35式 1-((50-1)/(50-2-1))* (1-summary(reg0)$r.squared) #自由度調整済み決定係数 summary(reg0) #Adjusted R-squaredに表示される summary(reg0)$adj.r.squared #summary(lm())$adj.r.squaredで取り出す ##(e) 8.37式 (summary(reg0)$r.squared/2)/((1-summary(reg0)$r.squared)/(50-2-1)) summary(reg0) #F-statistic 10.74 on 2 and 47 DF p-value: 0.0001440 summary(reg0)$fstatistic #summary(lm())$sigmaでF値と自由度を取り出す summary(reg0)$fstatistic[1] #f値 summary(reg0)$fstatistic[2] #分母の自由度 summary(reg0)$fstatistic[3] #分子の自由度 pf(summary(reg0)$fstatistic[1],df1=summary(reg0)$fstatistic[2],df2=summary(reg0)$fstatistic[3],lower.tail=F) ################################################ #Step11:個々の独立変数の寄与の評価 p.255- ################################################ ##(a) 独立変数の寄与の検定 8.44式 (summary(reg0)$r.squared - summary(reg3)$r.squared)/ ((1-summary(reg0)$r.squared)/(50-2-1)) anova(reg3,reg0) #anova(回帰分析の結果1, 回帰分析の結果2)で表示される ##(b) 標準誤差を用いた検定 8.47式 summary(reg0) #表示される偏回帰係数(Estimate)を標準誤差(Std Error)で割ると #t値が求まる 5.1716/1.6647 #例1 0.3027/0.1467 #例2 1.1259/0.4713 #例3