################################################ #南風原（2002）第3章 ################################################ #Step1：データセットの作成 ################################################ ##(a)作業フォルダを指定 setwd("C:/data/haebara/chap02") ##(b)データをcsvから読み込み，data1という名前をつける ##ラベルはcsvに書き込まれているので，header=Tと指定 data1 <- read.csv("http://blue.zero.jp/yokumura/R/haebara/chap1.csv",header=T) ##(c)読み込んだデータを確認 data1 ##(d)男性だけのデータ，女性だけのデータを作成 逸脱男 <- subset(data1,性別=="男",select=c(小6,中2,変化量)) 逸脱女 <- subset(data1,性別=="女",select=c(小6,中2,変化量)) 逸脱全体 <- subset(data1,select=c(小6,中2,変化量)) ##(e)追加パッケージをなんでもかんでも追加する(20-30分時間がかかる) # 1. CRANから一覧を入手 # packs<-CRAN.packages(contriburl=contrib.url("http://cran.md.tsukuba.ac.jp/")) # 2. 一覧をもとにダウンロードとインストール # install.packages(packs[,1],contriburl=contrib.url("http://cran.md.tsukuba.ac.jp/")) ################################################ #Step2：散布図と共分散 図3-1 ################################################ ##(a)attachで，「逸脱男」の中の変数を指定できるようにする attach(逸脱男) ##(b)plot(y軸の変数~x軸の変数)でプロット ## xlim=c(x軸の最小値,x軸の最大値) ## ylim=c(y軸の最小値,y軸の最大値) ## abline(h=水平線を描く位置) ## abline(v=垂線を描く位置） plot(中2~小6, ylim=c(0,22),xlim=c(0,22) ) abline(h=mean(中2)) abline(v=mean(小6)) ##(c)共分散の定義 scov <- function(x,y){ 1/length(x)*sum((x-mean(x))*(y-mean(y))) } ##(d)不偏共分散の定義 cov2 <- function(x,y){ (1/(length(x)-1))*sum((x-mean(x))*(y-mean(y))) } ##(e)定義したscovとcov2と，元々ある関数covの結果の比較 ##元々ある関数covは不偏共分散を算出することが確認できる scov(小6,中2) #標本共分散 cov2(小6,中2) #不偏共分散 cov(小6,中2) #不偏共分散 ##(f)作業が終了したら，attachしたデータに対しdetachする detach(逸脱男) ################################################ #Step3：共分散から相関係数へ p.47-49 ################################################ ##(a)相関係数の定義 scor <- function(x,y){ 1/length(x)*sum((x-mean(x))*(y-mean(y)))/ ##分子は共分散 (sqrt(1/length(x)*sum((x-mean(x))^2))*sqrt(1/length(y)*sum((y-mean(y))^2))) } ##(b)定義したscorと元々ある関数corの結果の比較 attach(逸脱男) scor(中2,小6) cor(中2,小6) #scorとcorは同一 detach(逸脱男) ##(c)gendat2（標本サイズ，相関係数）で，好みの相関係数を持つ2変量生データ作成する ##ただし，相関が-1と1の場合は，何故か使えないので，近似値を使う ##gendat2関数の中身（計算の意味）を理解する必要はない gendat2 <- function(nc, r) { # 仮のデータ行列を作る。この時点では変数間の相関は近似的に0である。 z <- matrix(rnorm(2*nc), ncol=2) # 主成分分析を行い，主成分得点を求める。この時点で変数間の相関は完全に0である。 res <- eigen(r2 <- cor(z)) coeff <- solve(r2) %*% (sqrt(matrix(res$values, 2, 2, byrow=TRUE))*res$vectors) z <- t((t(z)-apply(z, 2, mean))/sqrt(apply(z, 2, var)*(nc-1)/nc)) %*% coeff # コレスキー分解の結果をもとに，指定された相関係数行列 を持つように主成分得点を変換する。 z %*% chol(matrix(c(1, r, r, 1), ncol=2)) } ##(d)gendat2()の効果を確認 gendat2(300,.5) x <- gendat2(300,.5) cor(x) plot(x) ##(e)par(mafrow=c(行の数,列の数))で，複数図を描くことを指定する par(mfrow=c(3,4)) plot(gendat2(300,.999999999999999),ylab="") plot(gendat2(300,.9),ylab="") plot(gendat2(300,.8),ylab="") plot(gendat2(300,.7),ylab="") plot(gendat2(300,.6),ylab="") plot(gendat2(300,.5),ylab="") plot(gendat2(300,.4),ylab="") plot(gendat2(300,.3),ylab="") plot(gendat2(300,.2),ylab="") plot(gendat2(300,.1),ylab="") plot(gendat2(300,0),ylab="") plot(gendat2(300,-.8),ylab="") ################################################ #Step4：逸脱行動データへの適応 p.55-57 ################################################ ##(a)標準偏差関数の定義(ssd)と命名 ssd <- function(x){ #BeginFnc sqrt(1/length(x)*sum((x-mean(x))^2)) } ##(b)逸脱男の回帰係数と切片の確認 attach(逸脱男) b <- cor(中2,小6)*ssd(中2)/ssd(小6) ##回帰係数 b a <- mean(中2)-b*mean(小6) ##切片 a ##注意 切片の値が微妙に合わないが，丸め誤差の影響である mean(中2)-0.59*mean(小6) ##教科書で出てくる切片 detach(逸脱男) ##(c)逸脱女の回帰係数と切片の確認 attach(逸脱女) b <- cor(中2,小6)*ssd(中2)/ssd(小6) ##回帰係数 b a <- mean(中2)-b*mean(小6) ##切片 a detach(逸脱女) ##(d)lm(従属変数~独立変数,data=データセット名)関数を使うと一発 reg1 <- lm(中2~小6,data=逸脱男) reg2 <- lm(中2~小6,data=逸脱女) ##(e)plot関数を用いて，散布図を描く plot(逸脱男$中2~逸脱男$小6, xlim=c(0,20),ylim=c(0,20), xlab="小6", ylab="中2", pch="男",col="blue", main="逸脱行動データにおける男女の回帰直線の比較" ) lines(逸脱女$中2~逸脱女$小6,type="p",pch=24,col="red") ##(f)abline(切片,回帰係数)で回帰直線を描く ##また，coef()で回帰分析の結果を取り出す abline(coef(reg1),col="blue") abline(coef(reg2),col="red") ##(g)「car」パッケージをインストールすると， ##scatterplot(従属変数~独立変数 | 質的変数, smooth=F, data=データセット名) ##で同じことが1行のプログラムで実行可能である library(car) scatterplot(中2~小6|性別,smooth=F,data=data1) ################################################ #Step6：予測値と残差の平均と相関 p.57-58 ################################################ ##(a)予測値（3.16式）を確認 attach(逸脱男) y.predict <- mean(中2)+0.5848126 *(小6-mean(小6)) #予測値の計算3.15式 mean(y.predict) #予測値の平均を確認 mean(中2) #従属変数の平均と一致する ##(b)残差（3.17式）を確認 y.res <- 中2 - y.predict ##残差の計算3.17式 mean(y.res) ##残差の平均は0になる cor(小6,y.res) ##独立変数と残差の相関は0 cor(y.predict,y.res) ##予測値と残差の相関は0 plot(y.res~小6, xlab="独立変数",ylab="残差",main="独立変数と残差の相関" ) ##図3-6を描画 abline(lm(y.res~小6)) detach(逸脱男) ##(c)実は，lm()関数の結果で，予測値と残差は算出されている reg1$fitted.values ##予測値 reg1$residuals ##残差 predict(reg1) ##予測値 residuals(reg1) ##残差 ################################################ #Step7：予測の標準誤差 p.62-63 ################################################ ##(a)標準誤差の確認 (3.25)式 attach(逸脱男) se.pred <- ssd(中2)*sqrt(1-cor(中2,小6)^2) ##(b)標準誤差の確認 (3.26)式 se.pred*sqrt(length(中2)/(length(中2)-2)) ##(c)summary()を使った標準誤差の確認 summary(reg1) ##3.26式を使ったものが「Residual standard error」として ##算出される ##(d)相関係数と予測の標準誤差の関係　図3-7 r <- seq(0,1,by=.01) ##0から1まで，0.01の間隔のベクトルを作成 se.sy <- sqrt(1-r^2) ##y軸 plot(se.sy~r, type="l",main="相関係数と予測の標準誤差の関係") detach(逸脱男) ################################################ #Step8：回帰とは p.63-65 ################################################ ##(a)ゴルトンのデータをシュミレーション golton <- gendat2(200,0.6) golton <- as.data.frame(golton) names(golton) <- c("父親","息子") golton$父親 <- 15*(golton$父親-mean(golton$父親))/ssd(golton$父親)+100 golton$息子 <- 15*(golton$息子-mean(golton$息子))/ssd(golton$息子)+100 apply(golton,2,mean) #平均 apply(golton,2,ssd) #標準偏差 cor(golton) #相関係数の確認 ##(b)回帰分析 reg3 <- lm(息子~父親,data=golton) plot(golton) abline(reg3) ##(c)平均への回帰の例 reg.ex <- function(x){ 100+0.6*(x-100) } reg.ex(115) reg.ex(85) ################################################ #Step9：回帰とは p.62-63 ################################################ ##(a)父親の知能指数が100以上の者を抽出 golton2 <- subset(golton,父親>=100) ##(b)図3-9を描画 plot(golton2,xlim=c(55,145),ylim=c(55,145)) abline(reg3) reg4 <- lm(息子~父親,data=golton2) abline(reg4) graphics.off() ##(c)表3-1を計算（数値は教科書と同一にならないが， ##元々のデータがシュミレーションなので仕方が無い） apply(golton,2,mean) apply(golton2,2,mean) apply(golton,2,ssd) apply(golton2,2,ssd) cor(golton) cor(golton2) lm(息子~父親,data=golton) lm(息子~父親,data=golton2) lm(父親~息子,data=golton) lm(父親~息子,data=golton2)