#####################################################################
#####################################################################
#南風原朝和（2002）心理統計学の基礎 有斐閣   
#第3章 相関関係の把握と回帰分析 
#####################################################################
#####################################################################

#-------------------------------------------------------------------#
#Step1: データハンドリング
#-------------------------------------------------------------------#

##(a) データと変数名リストの読み込み 
setwd("Z:/workshop2006/haebara")
dat      <- read.table("tab1_1.csv", header=T, sep=",")
labeldat <- read.table("tab1_1label.csv", header=T, sep=",", na="")


##(b) 質的変数の代入
attach(labeldat)
  y <- as.character(label[!is.na(code1)])     #選択肢がある変数名をyとおく 
  nalt <- paste("code", 1:2, sep="")          #選択肢の変数名をnaltとおく
  for(i in 1:length(y)){
    alternative <- labeldat[label == y[i], nalt]  #選択肢がある変数名の選択肢をalternativeとおく
    alternative <- as.character(t(alternative)) 
    dat[, y[i]] <- factor(dat[, y[i]], labels=alternative)  #質的変数を代入する 
  }
detach(labeldat)


##(c) idが一致するものを展開する
datwide <- reshape(dat, idvar="id", timevar="time", direction="wide", v.names="point")


##(d) 変化量の得点を算出する 
datwide$diff <- datwide$point.中2 -  datwide$point.小6
datwide         #表1.1 



#-------------------------------------------------------------------#
#Step2：共分散と相関係数 (pp.43--51)
#-------------------------------------------------------------------#

#【宿題 1】_____________________________________________
#・plot()を利用して図3-1を再現すること。
#  ただし， 図中の斜線と数式を書くことは，要件に含めない。 
#・この際，グラフィックパラメータは，
#  pch, xlim, ylim, xlab, ylab, xaxt = "n", yaxt = "n", tcl, を変更すること。
#・さらに，axis(), seq(), abline() 関数を利用すること 
#・(a)(b)(d) 欄を読み，(c) の部分にplot()などを利用して足りない情報を埋めること。 
#_______________________________________________________

##(a) 宿題1-4で利用するデータ 
attach(datwide)
  x <- point.小6[sex == "男"]
  y <- point.中2[sex == "男"]  
  xbar <- mean(x)
  ybar <- mean(y)
detach(datwide)


##(b) dir.create(path="フォルダ名", showWarnings=F 警告を出さない)
##    でフォルダを作成する 
dir.create(path="fig", showWarnings=F)


##(c) windows(width = インチ, height = インチ)
##    で図を描画する大きさを決定する。デフォルトは，幅，高さ共に6インチ。 
windows(width=12, height=7)
plot()
axis()
axis()
abline()


##(d) savePlot(filename=ファイル名, type="wmf" 保存形式)
##    で図を保存する。
##    図の保存形式はMS Wordを利用している場合，wmfかemfを利用する場合が多い。
##    LaTeXを利用している場合は，psよりも，
##    mediabb.styを利用してpdfで保存する方が経験的に良い (コンパイルが早い)。 
##   さらに，相対パスと絶対パスの指定の仕方を覚えること。 
##   たとえば，http://www.shoshinsha.com/hp/1hour/know/pass.html を参照すること。 
##   基本的には，相対パスでの指定を常用する方が不都合が少ない。 
  getwd()           #現在の作業ディレクトリを確認する
  savePlot(filename= "Z:/workshop2006/haebara/fig/fig3_1a", type="wmf")    #絶対パスで保存する 
  savePlot(filename= "./fig/fig3_1b", type="wmf")                          #相対パスで保存する
graphics.off()                                                             #図を閉じる 



#【宿題 2】_____________________________________________
#以下の3つの問いに答えるプログラムを書くこと。 
#(1) xとyともに平均値を超える点は何個あるのか。
#(2) xとyともに平均未満の点は何個あるのか。
#(3) 一方の変数で平均を超え，他方の変数で平均を下回る点は何個あるのか。 
#・関数としてsum()，
#  比較演算子として >, < 
#  論理演算氏として &, | 
# を利用すること。 
#・これらの使い方は，Chap02.Rの【宿題 21】と似ている。 
#_______________________________________________________



#【宿題 3】_____________________________________________
#・共分散を求める関数 scov(3.1式)を作成すること 
#  引数は，x, yとすること 
#_______________________________________________________
scov <- function(x, y){ 
      return()
      }
scov(x, y)



#【宿題 4】_____________________________________________
#・不偏共分散を求める3.2式を書き，
#  関数cov()と同一の解が得られることを確認すること 
#_______________________________________________________
cov(x, y)



#【宿題 5】_____________________________________________
#・plot()を利用して図3-2を再現すること。
#・この際，グラフィックパラメータは，
#  pch, xlim, ylim, xlab, ylab, xaxt = "n", yaxt = "n", tcl, を変更すること。
#・さらに，axis(), seq(), abline(), lm() 関数を利用すること 
#・また，作図した結果を相対パスを利用して，
#  「fig3_2」というファイル名で，
#  何らかのファイル形式で保存すること。
#  保存先は，「fig」フォルダ内に保存すること。 
#_______________________________________________________

##(a) 宿題5で利用するデータ 
x <- c(4, 6, 7, 9, 10, 12, 13, 14, 15, 16, 18, 20)
y <- c(9.4,10.6,11.2,12.4,13.0,14.2,14.8,15.4,16.0,16.6,17.8,19.0)


##(b) 図示 
plot()
axis()
axis()
abline()

graphics.off()



#【宿題 6】_____________________________________________
#・xとyの共分散の絶対値がとりうる最大の値を求めること(3.3式)
#・さらに，実際に得られた共分散の値は，
#  その最大値の約6割の大きさであることを確認するプログラムを書くこと
#・関数は，ssd()，scov()を利用すること。 
#_______________________________________________________

##(a) 宿題6-7で利用するデータ 
attach(datwide)
  x <- point.小6[sex == "男"]
  y <- point.中2[sex == "男"]  
detach(datwide)


##(b) 宿題6で利用する関数 (2.13 式)
ssd <- function(x){
  N <- length(x)
  return(sqrt(1/N * sum((x - mean(x))^2)))
}


##(c) 3.3式など 



#【宿題 7】_____________________________________________
#・相関係数を求める3.4式を書き，
#  関数cor()と同一の解が得られることを確認すること 
#・ 関数は，ssd()，scov()を利用すること。 
#_______________________________________________________
cor(x, y)



#【宿題 8】_____________________________________________
#・plot()を利用して図3-3を再現すること。
#・(a)の部分は理解する必要はない。(b)，(c)の部分を読み，
#  (d)の部分にpar(), for()などを利用して図示すること 
#・グラフィックパラメータは，
#  sub, pch, xaxt, yaxt, mfcolを変更すること。
#・また，作図した結果を相対パスを利用して，
#  「fig3_3」というファイル名で，
#  何らかのファイル形式で保存すること。
#  保存先は，「fig」フォルダ内に保存すること。 
#_______________________________________________________

##(a)gendat2（標本サイズ，相関係数）で，好みの相関係数を持つ2変量生データ作成できる 
##  ただし，相関が-1と1の場合は，何故か使えないので，近似値 (例えば，-.99と.99)を使う
##  gendat2関数の中身（計算の意味）を理解する必要はない
gendat2 <- function(nc, r){
    z <- matrix(rnorm(2*nc), ncol=2)
    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))
}


##(b)gendat2()の効果を確認する 
gendat2(nc=300, r=.5)   #任意の相関係数(r=.5)を持つ標本サイズ300の2変量データを作成
x <- gendat2(300,.5)    #データをxとして保存 
cor(x)                  #2変量の相関係数を確認 
plot(x, sub="r = .5", xlab="", ylab="", xaxt="n", yaxt="n")
graphics.off()


##(c) 宿題8で利用するデータ
r <- c(.999, .9, .8, .7, .6, .5, .4, .3, .2, .1, 0, -.8)  #相関係数 

#formatC(数値, digits=有効桁数, format="f")で任意の有効桁数の文字列を作成できる 
sub <- formatC(r, digits=1, format="f")     

#paste()で 先に作成した文字列に，"r ="を付ける
sub <- paste("r = ", sub, sep="")


##(d) par(mfcol=c()), for()を利用して描画する
##    for() で，ループ変数をi，リストを「1からrの長さ」とする。
##    <繰り返す式>の中には，gendat2を利用して，
##    任意の相関係数を持つ2変量データを作成する。
##    ループ変数は，rとsubに利用する 
windows(height=8, width=7)
par() 
  for(i in ){
    plot()
  }

graphics.off()



#【宿題 9】____________________________________________
#・3.8式と3.9式が成立することを確認すること
#・(b) 欄に記載してあるx, yを, 
#  (a) 欄のlinear.transで適当な値に線形変換した上で
#   x1, y1というオブジェクト名で保存して，
#   3.8式と3.9式の左辺と右辺を確認すること。
#・なお，3.8式は(c)欄に，3.9式は(d)欄にプログラムを記載すること 
#・linear.trans()は，Chap02.Rの【宿題 19】で作成したものである 
#_______________________________________________________

##(a) 宿題9で利用する関数 (2.20式)
linear.trans <- function(x, c, d){
    return(c * x + d)
}


##(b) 宿題9で利用するデータ 
attach(datwide)
  x <- point.小6[sex == "男"]
  y <- point.中2[sex == "男"]  
detach(datwide)

x1 <- linear.trans()
y1 <- linear.trans()


##(c) 3.8式 
                    #左辺 
                    #右辺 


##(d) 3.9式 
                    #左辺 
                    #右辺 



#【宿題 10】____________________________________________
#・3.10式が成立することを確認すること
#_______________________________________________________

##(a) 宿題10で利用するデータ 
sim1 <- gendat2(100, .4) #相関 r= .4，標本サイズ100のデータ
sim2 <- gendat2(100, .8) #相関 r= .8，標本サイズ100のデータ
x <- sim1[,1]
y <- sim1[,2]
z <- sim2[,1]
w <- sim2[,2]
xy <- x + y
zw <- z + w


##(b) 3.10式 
                                                    #左辺 
                                                    #右辺 




#-------------------------------------------------------------------#
#Step3：回帰直線のあてはめ (pp.52--57)
#-------------------------------------------------------------------#

#【宿題 11】____________________________________________
#・(a) にあるデータを基に，
#  3.13式と3.14式を参考に，回帰係数と切片を求め，
#  オブジェクトb, aとして保存すること。 
#・また，(c) 欄に，lm()とcoef()を利用して，
#  同様の結果が得られることを確認すること  
#  lm()の結果は，lm1というオブジェクト名で保存すること 
#・なお，丸め誤差により，一部，教科書で算出される値と相違が出る。 
#_______________________________________________________
#-------------------------------------------------------
#lm()の利用法 
#・lm1 <- lm(従属変数 ~ 独立変数, data=データセット名)とした結果を保存する。
#・従属変数と独立変数は，データセットの中の変数名を直接利用できる。
#　つまり，
#　lm(point.中2~point.小6, data=datwide)
#　と書くと，全員のデータの回帰分析を行うことができる。
#・その後，
#  summary(lm1), coef(lm1), residuals(lm1), predict(lm1)などの関数で結果を呼び出す。
#  詳しくは，マニュアルとR-Tipsを参照すること。
#・条件にあった被験者だけを取り出して分析したい場合は，
#  subset=(条件式)を指定する。
#  たとえば，男性だけを取り出したい場合は， 
#  lm(point.中2~point.小6, data=datwide, subset=(sex=="男"))
#  lm(point.中2~point.小6, data=datwide, subset=(sex!="女"))
#  と指定する。
#・条件が複数ある場合も，もちろん対応できる。
#  たとえば，変化量の絶対値が4以上の男性だけを取り出して分析したい場合は，
#  lm(point.中2~point.小6, data=datwide, subset=(sex=="男" & abs(diff) >= 4))
#  と指定する 
#-------------------------------------------------------

##(a) 宿題11-12で利用するデータ 
attach(datwide)
  x <- point.小6[sex == "男"]
  y <- point.中2[sex == "男"]  
detach(datwide)

##(b) 回帰係数と切片 
b                       #回帰係数 
a                       #切片 


##(c) lm()とcoef()



#【宿題 12】_____________________________________________
#・plot()を利用して図3-4を再現すること。
#  ただし， 図中の点線と数式を書くことは，要件に含めない。 
#・この際，グラフィックパラメータは，
#  pch, xlim, ylim, xlab, ylab, xaxt = "n", yaxt = "n", tcl, を変更すること。
#・さらに，axis(), seq(), abline()を利用すること
#・回帰直線は，点(xの平均, yの平均) を通ることを確認するために，
#  abline()を利用して，yの平均を水平線，xの平均を垂線として追加すること 
#・また，作図した結果を相対パスを利用して，
#  「fig3_4」というファイル名で，
#  何らかのファイル形式で保存すること。
#  保存先は，「fig」フォルダ内に保存すること。 
#_______________________________________________________
plot()
axis()
axis()
abline()

graphics.off()

 

#【宿題 13】____________________________________________
#・lm(formula, data=datwide, subset=sex=="男") と
#  いう関数を利用して，
#  逸脱行動データの女子の場合の回帰直線を求めること
#・女子のデータは，lm2というオブジェクト名で保存すること  
#_______________________________________________________
lm1 <- lm(point.中2~point.小6, data=datwide, subset=sex=="男")
lm2



#【宿題 14】_____________________________________________
#・plot()を利用して図3-5を再現すること。
#・この際，グラフィックパラメータは，
#  pch, xlim, ylim, xlab, ylab, xaxt = "n", yaxt = "n", tcl, を変更すること。
#・さらに，axis(), seq(), abline() 関数を利用すること
#・また，作図した結果を相対パスを利用して，
#  「fig3_5」というファイル名で，
#  何らかのファイル形式で保存すること。
#  保存先は，「fig」フォルダ内に保存すること。 
#_______________________________________________________
plot(, data=datwide)
axis()
axis()

graphics.off()




#-------------------------------------------------------------------#
#Step4：回帰分析における予測値と残差の性質 (pp.57--67)
#-------------------------------------------------------------------#

#【宿題 15】____________________________________________
#・(a) のデータを基に，(b)-(j)の問いに答えること
#_______________________________________________________

##(a) 宿題15,17で利用するデータ 
attach(datwide)
  x <- point.小6[sex=="男"]
  y <- point.中2[sex=="男"]
  b <- cor(x, y) * ssd(y) / ssd(x)  #3.13式 
  a <- mean(y) - b * mean(x)        #3.14式 
  lm1 <- lm(y ~ x)
detach(datwide)


##(b) 予測値 yhatを求めること (3.11式) 
yhat


##(c) 予測値yhatの平均値は従属変数yの平均値と等しいことを
##    確認すること (3.16式)
                        #左辺 
                        #右辺 


##(d) predict(lm.object) で予測値yhatを算出できること
##   を確認すること 
predict()
yhat                    #等しい


##(e) 従属変数yと予測値yhatのズレ(残差e)を求めること(3.17式)
e


##(f) residuals(lm.object)で残差eを算出できることを
##    確認すること 
residuals()
e

##(g) 残差eの平均は0となることを確認すること(3.18式) 
#     注)e-16と標記されるのは，10^(-16) = ほとんど0を意味する 


##(h) 残差eと独立変数xの相関は0となることを確認すること (3.19式) 


##(i) 図3-6を再現すること。
##    なお，作成した図をfig3_6というファイル名で保存する。 
plot()
axis()
axis()

graphics.off()


##(j) 予測値yhatと残差eの間の相関は0となることを確認すること (3.20式) 



#【宿題 16】____________________________________________
#・library(UsingR)にあるdata(fat)の
#  weight（体重）, height（身長）, body.fat（体脂肪率）
#  の3変数を使用して(b)-(e) の問いに答えること 
#_______________________________________________________

##(a) 宿題16で利用するデータ 
library(UsingR)
data(fat)
?fat

##(b) 身長と体重の散布図を描くこと 
plot(, data=fat)


##(c) 身長から体重を予測する回帰分析を行い，
##    その結果を，lm3として保存すること 
lm3


##(d) 体重と体脂肪率の相関を求めること


##(e) 身長から体重を予測する回帰分析の残差と
##    体脂肪率の相関を求めること 



#【宿題 17】____________________________________________
#・(b)-(h)の問いに答えること
#_______________________________________________________

##(a) 宿題17で利用する関数 (2.12式)
svar <- function(x){
  N <- length(x)
  return(1/N * sum((x - mean(x))^2) )
}


##(b) 従属変数yの分散は，
##    予測値yhatの分散と残差eの分散の和に等しいことを
##    確認すること (3.22式) 


##(c) 予測値yhatの分散は，3.23式になることを確認すること


##(d) 残差eの分散は3.24式になることを確認すること 


##(e) 残差の標準偏差は3.25式になることを確認すること 


##(f) 図3-7を再現すること 
##    なお，作成した図をfig3_7というファイル名で保存すること。 
plot()
axis()
axis()

graphics.off()


##(g) 予測の標準誤差3.26式を書くこと 


##(h) 予測の標準誤差は
##    summary(lm.object)
##    predict(lm.object, se.fit=T)で算出できることを確認すること 
summary()
predict(, se.fit=T)




#-------------------------------------------------------------------#
#Step5：相関係数と回帰係数の性質の違い (pp.67--72)
#-------------------------------------------------------------------#

#【宿題 18】____________________________________________
#・(a) のデータを基に，(b)-(l)の問いに答えること 
#_______________________________________________________

##(a) 宿題18で利用するデータ 
fig3.8 <- read.table("fig3_8.csv", sep=",", header=T)
head(fig3.8)
tail(fig3.8)  #父親と，息子の知能指数（N=200）のデータ 


##(b) 知能指数の平均値，標準偏差を求めること 


##(c) 父親と息子の知能指数の相関係数を求めること 


##(d) 父親の知能指数から息子の知能指数を予測する回帰分析を行い
##    その結果を，lm4というオブジェクトに保存すること 
lm4 <- lm(, data=fig3.8)


##(e) 図3.8を再現すること。
##    回帰直線は，変数yの変数xへの回帰直線だけ描けばよい
##    なお，作成した図をfig3_8というファイル名で保存すること。 
plot()
abline()
graphics.off()


##(f) predict(lm.object, newdata)
##    で任意の値の，予測値を求められることを確認すること 
predat1 <- data.frame(father=c(115, 85, 100))
predict(, newdata=predat1)


##(g) 息子の知能指数から父親の知能指数を予測する回帰分析を行い
##    その結果を，lm5というオブジェクトに保存すること 
lm5 <- lm(, data=fig3.8)


##(h) 息子の知能指数が，115, 85, 100の場合の，
##    父親の知能指数の予測値を求めること 
predat2 <- data.frame()
predict(, newdata=predat2)


##(i) subset()を利用して，父親の知能指数が100以上のデータを取り出し，
##    fig3.9というオブジェクト名で保存すること 
fig3.9 <- subset(fig3.8, )


##(j) fig3.9のデータを基に
##    父親の知能指数から息子の知能指数を予測する回帰分析を行い
##    その結果を，lm6というオブジェクトに保存すること 
##    息子の知能指数から父親の知能指数を予測する回帰分析を行い
##    その結果を，lm7というオブジェクトに保存すること 
lm6    <- lm(, data=fig3.9)
lm7    <- lm(, data=fig3.9)


##(k) 図3.9を再現すること 
##    なお，作成した図をfig3_9というファイル名で保存すること。 
plot()
abline()
graphics.off()


##(l) 表3.1 を再現すること 
##    なお，結果は，tab3.1というオブジェクトに代入すること 
tab3.1 <- data.frame(matrix(ncol=2, nrow=7))
rownames(tab3.1) <- c("xの平均", "yの平均", "xの標準偏差", "yの標準偏差", "相関係数", "回帰係数 b","回帰係数 b'")
colnames(tab3.1) <- c("選抜前", "選抜後")



#【宿題 19】____________________________________________
#・library(DAAG)にあるdata(socsupport)の
#  BDI（抑うつ）, emotional（情緒的サポート利用可能性）
#  の2変数を使用して(b)-(c) の問いに答えること 
#_______________________________________________________

##(a) 宿題19で利用するデータ 
library(DAAG)
data(socsupport)
head(socsupport)
tail(socsupport)
?socsupport


##(b) 抑うつと情緒的サポートの
##    ・相関係数，
##    ・散布図，
##    を求めること
##　　さらに，抑うつから情緒的サポートを予測する回帰分析を行うこと 


##(c) 抑うつの得点が10以上の者を選抜した際の，
##    抑うつと情緒的サポートの
##    ・相関係数，
##    ・散布図，
##    を求めること
##　　さらに，抑うつから情緒的サポートを予測する回帰分析を行うこと
socsupport2 <- subset(socsupport, )




#-------------------------------------------------------------------#
#Step6：測定の妥当性と信頼性 (pp.76--84)
#-------------------------------------------------------------------#

#【宿題 20】____________________________________________
#・(a) のデータを基に，(c)-(h)の問いに答えること 
#_______________________________________________________

##(a) 宿題20で利用するデータ 
t <- rnorm(500)           #真値をtとする 
x <- jitter(t, amount=1)  #jitter(x, amount=1)でノイズを発生させ，
                          #測定誤差を含むxを作成する。
                          #amountの大きさを変えることにより，
                          #ノイズの量が大きくなる。 
e <- x - t                #測定値xと真値tの差は誤差eとなる 


##(b) 宿題で利用する関数 (2.12式)
svar <- function(x){
  N <- length(x)
  return(1/N * sum((x - mean(x))^2) )
}


##(c) 測定値の分散は，
##    真値の分散と測定値の分散の和に等しくなることを確認すること(3.33式)
##    ただし，測定誤差と真値の相関が0とみなせない程度に応じて
##    3.33式は厳密には成り立たない 


##(d) 測定値xの信頼性を求めること (3.34式) 
##    その結果は，オブジェクトrxとして保存すること 
rx


##(e) jitter()を利用して，
##    3.32式であらわされる測定値xと同じ真値を持つ，
##    別の測定値x.primeを作成すること。
##    また，そのx.primeの測定誤差をe.primeとして保存すること 
x.prime
e.prime


##(f) xとx.primeの共分散の展開を確認すること(3.36式の1つ上の式)


##(g) 3.36式の関係が成り立つことを確認すること
##    ただし，測定誤差と真値の相関が0とみなせない程度に応じて
##    3.36式は厳密には成り立たない 


##(h) 3.37式の関係が成り立つことを確認すること
##    ただし，測定誤差の分散が等しくない程度に応じて
##    3.37式は厳密には成り立たない 



#【宿題 21】____________________________________________
#・(a) のデータを基に，(b)-(f)の問いに答えること 
#_______________________________________________________


##(a) 宿題21で利用するデータ 
sim1 <- gendat2(500, .4)  #相関が.4，N=500のデータを作成 
t <- sim1[,1]             ##真値 
x <- jitter(t, amount=1)  ##測定誤差を含むx 
e <- x - t                ##誤差 

t.prime <- sim1[,2]             ##真値 
y <- jitter(t.prime, amount=1)  ##測定誤差を含むy
e.prime <- y - t.prime          ##誤差 


##(b) x, yの共分散の展開を確認すること(3.38式の1つ上の式)


##(c) 測定値間の共分散が，
##    その真値間の共分散に等しくなることを確認すること (3.38式)  
##    ただし，3.38式の1つ上の式の第2項以下の測定誤差を含む共分散を
##    ゼロとみなせない程度に応じて
##    3.38式は厳密には成り立たない 


##(d) 測定値x, yの信頼性を求めること (3.34式) 
##    その結果は，オブジェクトrx, ryとして保存すること 
rx
ry


##(e) 測定値x,yの標準偏差と，信頼性rx, ryと真値の標準偏差の関係
##    を確認すること (3.39式) 


##(f) 測定値間の相関係数は，3.40式となることを確認すること 



#【宿題 22】____________________________________________
#・以下のプログラムは，2変数間の信頼性が等しい場合に，
#  信頼性と真値間の相関係数の値に応じて，
#  測定値の相関係数が変化することを示したものです (3.40式の一部を改変)
#・何をしているか，プログラムを読み解くこと 
#・描かれた図は何を意味するか理解すること 
#_______________________________________________________
rx <-  seq(0,  1, length=100)  #rxを信頼性，rtを真値間の相関係数とする
rt <-  seq(-1, 1, length=100)

eq3.40 <- function(rt, rx){
      return(rt * sqrt(rx * rx))  #3.40式 
      }
z <- outer(rt, rx, eq3.40)        #z方向の大きさを求める 
                                  #The R Tipsのpersp()の説明部分を参照 
#図示 
persp(rt, rx, z, 
      ticktype="detailed", 
      phi=30, theta=240, expand=.8,
      col="lightblue", 
      xlab="真値間の相関係数", ylab="信頼性",zlab= "測定値の相関係数")
graphics.off()