相関係数の区間推定
サンプル:日本統計学会 (編集); 日本統計学会公式認定 統計検定2級対応 統計学基礎, 東京図書, 2012,p125
身長のデータ
大学生
x<-c(172,167,184,175,176,175,170,180,170,179,167,175,174,162,165,163,170,169,165,175)
父親
y<-c(165,165,178,176,150,171,172,175,170,156,163,170,165,160,163,170,163,165,160,172)
cor.test関数を使用します
同じ数の要素をもつベクトル2つの相関関係、また相関係数の 95% 信頼区間,p値を求めます.Pearsonの相関係数をもとめるのであれば次の式を実行するのみでOKです.
cor.test(x,y)
Pearson's product-moment correlation
data: x and y
t = 1.4587, df = 18, p-value = 0.1619
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.1371017 0.6711054
sample estimates:
cor
0.3251458
cor.test(x,y)の引数
x, y:同じ大きさの数字ベクトル
alternative(hypothesis): "two.sided", "greater", "less".
method : "pearson(相関係数cor)", "kendall"(相関係数tau), "spearman"(相関係数rho).
(r-de様のコメントを参考に修正しました)
次に上記関数結果の95%信頼区間[-0.1371017 0.6711054]の求め方を忘れないように残しておきます
大学生の標準偏差
Vx<-(sum( (x-mean(x) )^2) )/20
Sx<-sqrt(Vx)
Sx
父親の標準偏差
Vy<-(sum( (y-mean(y) )^2) )/20
Sy<-sqrt(Vy)
Sy
大学生と父親の共分散
Sxy<-(sum((x-mean(x) )*(y-mean(y) ) ) )/20
標本の相関係数r
Sxy/(Sx*Sy)
r=0.325
標本の相関係数rの分布を考えてみます.
r=0の場合
t=r*(sqrt(n-2)/sqrt(1-r^2) )は自由度n-2のt分布に従うr≠0の場合(FisherのZ変換)
rの分布はρが0から離れると非対称性が強くなります.
Z={log((1+r)/(1-r) ) }/2と変換した場合、Zの分布は対称になります.
またZは近似的にN({log( (1+ρ)/(1-ρ) ) }/2,1/(n-3) )に従うことが知られています.
r=0.325をZ={log((1+r)/(1-r) ) }/2に代入してZ=0.337を求めます.
{log((1+ρ)/(1-ρ) )=Rとおくと
(0.337-R)=±1.96×標準偏差
0.337-1.96/√(17)<R<0.337+1.96/√(17)
2×(0.337-1.96/√(17))<log((1+ρ)/(1-ρ)<2×(0.337+1.96/√(17))
‐0.138×2< log((1+ρ)/(1-ρ) )<0.671×2
したがってρの95%信頼区間は
(exp(-0.276)-1)/(exp(-0.276)+1)
-0.137
(exp(1.626)-1)/(exp(1.626)+1)
0.671
また勉強して更新していきますのでコメントお願いします.
正規分布のグラフ
グラフ初歩の初歩
手探り状態でやってます
赤文字のみRで実行
正規分布
#1正規分布の乱数を100個生成
x<-rnorm(100)
#ヒストグラム,freq=Fで確率密度,ylimでy軸範囲
hist(x,freq=F,ylim=c(0,0.6))
#枠の処理、lty線種は実践、btyで左と下のみ
box(lty=1,bty = "l")
#確率密度のグラフ,add=Tで重ねる、x軸範囲-3~3
curve(dnorm(x),add=T,-3,3)
これらをオリジナルの関数として記憶させます
mynorm<-function(n,g)
{
x<-rnorm(n)
hist(x,freq=F,ylim=c(0,0.6),xlab=g)
box(lty=1,bty = "l")
curve(dnorm(x),add=T,-3,3)
}
mynorm(100,"n=100")
#四つのグラフを2行2列のまとめてみます
par(mfrow = c(2,2))
mynorm(30,"n=30")
mynorm(100,"n=100")
mynorm(300,"n=300")
mynorm(1000,"n=1000")
グラフに色をつける
library(RColorBrewer)
#RColorBrewerパッケージのサンプル
display.brewer.all()
ヒストグラムを塗ってみる
どの色セットを使用するかを指定する
cols <- brewer.pal(8,"Pastel1") # brewer.pal(何色、パレット名)
y<-c(1,2,3,4,5,6,7)
p<-c(2,3,4,5,4,3,2)
q<-c(2,3,4,5,4,3,2)
par(mfrow = c(3,3),mar = c(5, 4, 1, 4)) #余白 底辺、左、上、右の順
pos.x <- barplot(q,ylim=c(0,6),col=cols[1])
pos.x <- barplot(q,ylim=c(0,6),col=cols[2])
pos.x <- barplot(q,ylim=c(0,6),col=cols[3])
pos.x <- barplot(q,ylim=c(0,6),col=cols[4])
pos.x <- barplot(q,ylim=c(0,6),col=cols[5])
pos.x <- barplot(q,ylim=c(0,6),col=cols[6])
pos.x <- barplot(q,ylim=c(0,6),col=cols[7])
pos.x <- barplot(q,ylim=c(0,6),col=cols[8])
pos.x <- barplot(q,ylim=c(0,6),col=cols[9])
棒グラフの中央に散布図をプロットする方法 (pos.x)
y<-c(1,2,3,4,5,6,7)
p<-c(2,3,4,5,4,3,2)
q<-c(2,3,4,5,4,3,2)
par(mar = c(5, 4, 1, 4)) #余白 底辺、左、上、右の順
pos.x <- barplot(q,ylim=c(0,6))
points(pos.x, p)