理学療法士がまとめたノート

統計学備忘録 since2016

Rを使って統計学を勉強するブログです

相関係数の区間推定

サンプル:日本統計学会 (編集); 日本統計学会公式認定 統計検定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

 

また勉強して更新していきますのでコメントお願いします.

なぜt検定を繰り返してはいけないのか?

sample は「石村卓夫;分散分析のはなし,東京図書,1992,p137」

A<-c(12,14,16)
B<-c(13,15,17)
C<-c(15,17,19)
D<-c(17,19,21)

x<-c(rep("A",3),rep("B",3),rep("C",3),rep("D",3))
y<-c(A,B,C,D)
stripchart(y~x,vertical=T)

f:id:yoshida931:20170224165734p:plain

視覚的にはAD間、BD間には差があるように見えます

 

重要な前提です
Aの母平均をμa
Bの母平均をμb
Cの母平均をμc
Dの母平均をμd とします
「μa=μb=μc=μd」を検定する場合にt検定を繰り返してはいけないのです.
ある特別なペアの違いなどを確認したい場合にはt検定を使用します.研究内容により異なるので間違いないように.

 

参考 http://www.gen-info.osaka-u.ac.jp/MEPHAS/ave.html
2種類の既存薬AとBを組み合わせた配合薬Cの配合効果を評価する場合.この場合、既存薬A,Bのそれぞれの効果と配合薬Cの効果を比較します.ここで いいたいのはCが既存薬A、Bの両方よりも効果があるということです.「CがAよ りも優れている、かつCがBよりも優れている」ということを 示します.つまり帰無仮説はC=AかつC=Bとなり、A=B=Cとは異なります.このような場合には2標本t検定を繰り返して用いるほうが適切だと考えられます.

 

一元配置分散分析
oneway.test(y~x,var.equal = T)
#F = 3.6875, num df = 3, denom df = 8, p-value = 0.06214
#有意水準5%で差が「ない」という結果になりました

 

#全ての群の母分散は等しいと仮定してt検定を繰り返してみます
AB<-t.test(A,B,paired=FALSE, var.equal=T, conf.level=0.95)
AC<-t.test(A,C,paired=FALSE, var.equal=T, conf.level=0.95)
AD<-t.test(A,D,paired=FALSE, var.equal=T, conf.level=0.95)
BC<-t.test(B,C,paired=FALSE, var.equal=T, conf.level=0.95)
BD<-t.test(B,D,paired=FALSE, var.equal=T, conf.level=0.95)
CD<-t.test(C,D,paired=FALSE, var.equal=T, conf.level=0.95)

 

#それぞれの検定結果のP値のみを取り出してみます
str(AB3);str(AC3);str(AD3);str(BC3);str(BD3);str(CD3)
AB 0.573
AC 0.14
AD 0.0376
BC 0.288
BD 0.0705
CD 0.288

 やはりAD間には有意水準5%で「差がある」という結果になりました.しかし分散分析では差がないという結果でした.

 

分散分析の帰無仮説
有意水準0.5において「μa=μb=μc=μd」と仮定する

t検定の帰無仮説
有意水準0.5において「μa=μb」、
かつ有意水準0.5において「μa =μc」、
かつ有意水準0.5において「μa =μd」、
かつ有意水準0.5において「μb=μc」、
かつ有意水準0.5において「μb=μd」、
かつ有意水準0.5において「μc¬=μd」と仮定した検定になります.

したがってt検定を使った「μa=μb=μc=μd」の確率は(1-0.05)^6となります.少なくとも一組には差がある確率は1-(1-0.05)^6=0.2649081となり、有意水準26.5%の検定になってしまいます。つまりαエラーの危険性が高くなるのです.このような矛盾を解消するための手法が多重比較です(しかし、本当にt検定は繰り返して使ってはいけないのだろうか? 本当に有意水準の高い検定になっているのだろうか? 2群間比較にはt検定しかないのだが…疑問が残ります)

多重比較(Tukey)を行ってみます

TukeyHSD(aov(y ~ x),ordered = TRUE)
Tukey multiple comparisons of means
95% family-wise confidence level
factor levels have been ordered

Fit: aov(formula = y ~ x)

$x
  diff  lwr    upr   p adj
B-A 1 -4.2294189 6.229419 0.9252929
C-A 3 -2.2294189 8.229419 0.3245304
D-A 5 -0.2294189 10.229419 0.0609499
C-B 2 -3.2294189 7.229419 0.6297636
D-B 4 -1.2294189 9.229419 0.1441838
D-C 2 -3.2294189 7.229419 0.6297636

有意水準0.05ではどのペアにも差はありませんでした.

等分散性の検定

準備としてカイ二乗分布とF分布が必要になります

カイ二乗分布
確率変数Z1,Z2,…Znが互いに独立で標準正規分布にしたがうとき
W=Z1^2+Z2^2+…+Zn^2
の従う分布を自由度nのカイ二乗分布とよび、χ^2(n)で表します.
期待値E[W]=n, 分散V[E]=2n

カイ二乗分布の性質
① 二つの確率変数W1,W2が独立にχ^2(n1), χ^2(n2)に従うとき
W1+ W2は、χ^2(n1+n2)に従います.(カイ二乗分布の再生性)
② 母集団N(μ,σ^2)からの無作為標本X(X1,X2,…X3)について
Xiを標準すると Zi= (Xi-μ)/ σとなる
Ziは互いに独立で標準正規分布に従うので
∑{(Xi-μ)/ σ}^2 はカイ二乗分布χ^2(n)に従うことになります
③ 標本不偏分散
s^2={∑{(Xi-X平均) ^2 }/(n-1) より
(n-1) ×s^2=∑(Xi-X平均) ^2
母集団N(μ,σ^2)からの無作為標本Xi=(X1,X2,…X3)は互いに独立でN(μ,σ^2)に従うので
W=∑{(Xi-X平均) ^2}/ σ^2 = {(n-1) ×s^2}/σ^2
母平均μをX平均で置き換えると自由度がn-1になります
したがって{(n-1) ×s^2}/σ^2は自由度n-1のカイ二乗分布χ^2(n-1)に従います.

f:id:yoshida931:20170222184024p:plain

http://tips-r.blogspot.jp/2015/01/r.htmlより)


2標本の母分散の比の検定
散らばりを表す分散は非負の値なので、その大きさを考える場合には差ではなく比で考える方が適切です.

例)単位はmm
A群の身長の標準偏差はσa=5 (分散25)
B群の身長の標準偏差はσb=10 (分散100)
σa/σb=0.5 A群の方がB群に比べると2倍散らばりの小さい集団である
現実的には標準偏差が便利で理解しやすいのですが、推定での取り扱いは上述のような分散を使います.

 

F分布
条件
・確率変数U:自由度mのカイ二乗分布に従う
・確率変数V:自由度nのカイ二乗分布に従う
・UとVは独立である

フィッシャーの分散比 F= (U/m)/(V/n)

Fが従う確率分布を自由度(m,n)のF分布F(m,n)という

2標本の標本分散s1^2, s2^2を考えます
標本不偏分散s^2={∑{(Xi-X平均) ^2 }/(n-1) 
標本s1^2:{(m-1) ×s1^2}/σ1^2は自由度m-1のカイ二乗分布χ^2(m-1)に従います
標本s2^2:{(n-1) ×s2^2}/σ2^2は自由度n-1のカイ二乗分布χ^2(n-1)に従います
s1^2とs2^2は独立です

F={(m-1) ×s1^2}/(σ1^2)*(m-1)÷{(n-1) ×s2^2}/(σ2^2)*(n-1)
= (s1^2/σ1^2)÷(s2^2/σ2^2)
= (s1^2/s2^2)*(σ2^2/σ1^2)
の従う分布は、自由度(m-1,n-1)のF分布F(m-1,n-1)となります

 

 

それでは始めます
赤文字のみRで実行

2標本の場合

母分散の比の検定になります
二つの母集団からの標本をXi~N(μx,σx^2),Yi~N(μy,σy^2)とします.
帰無仮説:σx^2=σy^2、(σx^2)/(σy^2)=1
対立仮設:σx^2≠σy^2、(σx^2)/(σy^2)≠1

サンプルは「石村卓夫;分散分析のはなし,東京図書,1998,p69」
Xi<-c(12.2,18.8,18.2)
Yi<-c(26.4,32.6,31.3)
σx二乗<-sum((Xi-mean(Xi) )^2)/2;σx二乗
σy二乗<-sum((Yi-mean(Yi) )^2)/2;σy二乗

T<-σx二乗/σy二乗
1.25
F分布を使用します

qf(0.025,2,2,lower=F) # F 分布の上側5%点
39
qf(0.975,2,2,lower=F) # F 分布の下側5%点
0.026

T=1.25は棄却域には入らないので帰無仮説は棄却できない
標本をXiとYiの母分散は等しいと仮定できることになります

 

SPSSではlenene検定が使用されていますので以下に算出しておきます
Levene's Test もF分布を使用していますが、上記方法とは異なります
leveneTest(x,group)

Levene's Test for Homogeneity of Variance (center = median)
    Df  F value  Pr(>F)
group  1  0.0031  0.9585
    4

 

 

 

 

3標本以上の場合

 

繰り返し数が等しいとき:Hartleyの検定
各水準の分散(Si^2 )について
Si^2の最大分散÷最小分散が1に近づくことで等分散性を仮定します
最大分散/最小分散=max{Si^2}/min{Si^2}

 

f:id:yoshida931:20170222184102p:plain

http://www2.vmas.kitasato-u.ac.jp/lecture0/statistics/hartley.pdf


検定表の見方:Fmax(k , f)(k) = Fmax(水準数,繰り返し数-1)(水準数)
帰無仮説:各水準の分散は等しい

サンプルは「石村卓夫;分散分析のはなし,東京図書,1998,p108」
A<-c(12.2,18.8,18.2)
B<-c(22.2,20.5,14.6)
C<-c(20.8,19.5,26.3)
D<-c(26.4,32.6,31.3)
E<-c(24.5,21.2,22.4)

x<-c(A,B,C,D,E)

SuppDistsをインストトールします
library(SuppDists)
factor()カテゴリーを要素としたベクトル
group=factor(rep(c("A", "B", "C", "D", "E"), c(3, 3, 3, 3, 3)))
levels()でグループ化をstr()で構成を確認できます
fmax=max( by(x, group, FUN=var) ) / min( by(x, group, FUN=var) )  
5.70

by関数の引数(data,INDICES,FUN=,simplify=)
data:data.frameやmatrixもOK
INDICES:factorをsubsetに指定
FUN :引数dataに適用する関数を指定。ここでは分散varを指定。
simplify :結果のデータ形式。TRUE:scalar、FALSE:list or array。

 

積分布pmaxFratio
pmaxFratio( fmax, 繰り返し数の自由度, 水準数)
1-pmaxFratio(fmax, 2, 5)  #p値
0.8154272

帰無仮説は棄却されず、等分散性を仮定することができる

 

 

繰り返し数が異なるとき:Bartlettの検定
Rヘルプより
繰り返し数は異なりませんが同じサンプルで試行してみます
以下実施方法のみ記載
A<-c(12.2,18.8,18.2)
B<-c(22.2,20.5,14.6)
C<-c(20.8,19.5,26.3)
D<-c(26.4,32.6,31.3)
E<-c(24.5,21.2,22.4)

x<-c(A,B,C,D,E)
group=factor(rep(c("A", "B", "C", "D", "E"), c(3, 3, 3, 3, 3)))

bartlett.test(x,group)


Bartlett test of homogeneity of variances
data: x and group
Bartlett's K-squared = 1.2291, df = 4, p-value = 0.8733

 


Leveneの検定
最もよくみかける等分散性の検定だと思われます
水準の繰り返し数に縛りはありません

A<-c(12.2,18.8,18.2)
B<-c(22.2,20.5,14.6)
C<-c(20.8,19.5,26.3)
D<-c(26.4,32.6,31.3)
E<-c(24.5,21.2,22.4)

統計量(Wikipediaより)
W={(N-k)/(k-1)}×{∑Ni(Zi.-Z..)^2/∑∑(Zij-Zi.)^2}

数値はsample
k:比較する群の数水準
水準数i  1,2,3,4,5
繰り返しj  1,2,3
N:総サンプル数
Ni:第i群のサンプル数
自由度N-k=15-5=10
自由度k-1=5-1=4
Z.. 全データの平均値
Zi. i番目の水準内の平均値 
Yij:第i群のj番目の変数
Zij:| Yij – Yiの平均|、| Yij – Yiの中央値|


carパッケージをインストールします
x<-c(A,B,C,D,E)
group=factor(rep(c("A", "B", "C", "D", "E"), c(3, 3, 3, 3, 3)))

leveneTest(x,group)

Levene's Test for Homogeneity of Variance (center = median)
Df F value Pr(>F)
group 4 0.1256 0.9698
10

p値=0.9698で帰無仮説は棄却されず、等分散性を仮定した検定が可能になります

 

引用・参考
バイオインフォマティクス入門 http://bio-info.biz/
石村卓夫;分散分析のはなし,東京図書,1998

一元配置分散分析 (対応なし) F値の算出方法

更新日2017.3.13
更新日2017.6.22

赤文字の部分だけRにペーストして実行

一元配置分散分析とは、一つの因子に基づく“3つ以上の母平均の差の検定”です
因子Factor:実験結果に影響を与えると考えられる要因
水準Lebel:因子をいくつかに分ける基準 (水準=群

例題)
治療因子に基づく、患者20名の満足度データを比較する(繰り返し数=5)
因子Factor=治療、水準Level=治療A, 治療B, 治療C, 治療D
因子<c(水準)としてベクトルを作成

まずイメージを一覧表にしておきます

A<-c(10,6,10,9,10)
B<-c(10,5,5,12,4)
C<-c(5,4,11,4,6)
D<-c(9,5,2,3,1)
xx<-data.frame(A,B,C,D)
xx<-t(xx)
colnames(xx)<-c(1:5)

f:id:yoshida931:20170622170844p:plain

繰り返し数(n)=5(順序が入れ替わっても問題はない)

 

ALL<-c(A,B,C,D)  #全データ 
mean(ALL)     #総平均
Levels<-c(rep("A",5),rep("B",5),rep("C",5),rep("D",5))    # 水準oneway.test(ALL~Levels,var.equal = T)         #一元配置分散

One-way analysis of means
data: AllD and AllL
F = 2.5789, num df = 3, denom df = 16, p-value = 0.08979

それではF = 2.5789を求めてみます

 

分散分析の理解
分散分析は、観測データにおける変動(誤差変動水準内変動)と各要因による変動(水準間変動)に分解し効果を判定する方法です。(2元配置では交互作用による変動も考えます)。
1、 データの分散成分の平方和(総平方和全変動全平方和)を分解し、誤差による変動と要因の効果による変動を分離します。
2、 平方和を自由度で割った平均平方を算出します。
3、 統計量はF値になります。「要因の効果によって説明される平均平方(水準間変動)」を分子、「誤差によって説明される平均平方(水準内変動)」を分母にとった値が分散分析の統計量です。
4、 F検定。F分布により各効果の有意性を判定します。

以下、水準間変動、水準内変動という表現を使用します

 

まず各水準の平均間に差がない場合を考えてみます
赤線は全平均
F = 0.30185, num df = 3, denom df = 16, p-value = 0.8236

f:id:yoshida931:20170219143818p:plain


残っている変動は水準内の変動になります。つまり、「全変動=水準内の変動」ということになります。偏差の2乗和(平方和)を全2乗和、または全変動といいます。

逆に考えると水準間変動>水準内変動の場合には、「各水準間に差がある」と考えることができます。したがって「水準間変動÷水準内変動」が分散分析の統計量になります

ここで一元配置分散分析の数学モデル(X=x1,x2,…)を考えてみます
X = 総平均 + (水準平均―総平均) + (X―水準平均)
(X-総平均) = (水準平均―総平均) + (X―水準平均)

f:id:yoshida931:20170313111752p:plain


ここで差の平方をとって考えていきます
(X-総平均)^2 、 (水準平均―総平均) ^2 、(X―水準平均) ^2
例) 20個       4個         20個 
平方和∑にすると次の式が成り立ちます
n=各水準に含まれている繰り返し数(例題ではn=5)

∑(X-総平均)^2=n∑(水準平均―総平均) ^2 + ∑(X―水準平均) ^2

f:id:yoshida931:20170313112921p:plain

f:id:yoshida931:20170622172124p:plain

 

水準間変動  (A~Dで4水準なので df=4-1=3 )
y<-5*{(mean(A)-mean(ALL))^2+
(mean(B)-mean(ALL))^2+
(mean(C)-mean(ALL))^2+
(mean(D)-mean(ALL))^2}
=66.15


群間平方和 (
水準間変動の平均平方和 )
水準間変動:n∑(水準平均―総平均) ^2
各水準平均は正規分布N( 総平均+(水準平均―総平均)、σ^2/n)に従います。
したがって∑(水準平均―総平均) ^2÷(σ^2/n)は自由度(水準数-1)のχ二乗分布に従います。
水準間平均平方=∑(水準平均―総平均) ^2 ÷ (水準数-1)
66.15÷(4-1) = 22.05


水準内変動   (A~Dで4水準, 繰り返い数は5回なので df=4×(5-1)=16 )

sum( (A-mean(A))^2)+
sum( (B-mean(B))^2)+
sum( (C-mean(C))^2)+
sum( (D-mean(D))^2)
=136.8

郡内平方和、残差平方和 ( 水準内変動の平均平方和 )
水準内平方和:∑(X―水準平均) ^2
それぞれの水準における∑(X―水準平均) ^2 ÷ σ^2 はn-1のχ二乗分布に従う統計量になります。したがって∑∑(X―水準平均) ^2 ÷ σ^2は自由度、水準数×(n-1)のχ二乗分布に従います。
水準内平均平方=∑∑(X―水準平均) ^2 ÷水準数×(n-1)
136.8÷{4☓(5-1)} = 8.55


ポイント
誤差を正規分布に従う確率変数と考えます。したがって例題の水準A標本(10,6,10,9,10)は、正規母集団(6.55+9.0、σ^2)からの大きさ5の標本ということになります。∑(A-6.55)^2/σ^2 は自由度4のχ二乗分布に従う統計量となります。

分散分析の統計量F
    水準間平均平方 ÷ 水準内平均平方
= 22.05 ÷ 8.55 = 2.5789


Rで分散分析表を求めてみます
anova(aov(ALL~Levels))

Analysis of Variance Table
Response: ALL
       Df Sum Sq Mean Sq F value Pr(>F)
Levels 3 66.15 22.05 2.5789 0.08979 .
Residuals 16 136.80 8.55
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

 

f:id:yoshida931:20170313123012p:plain

正規分布のグラフ

グラフ初歩の初歩

手探り状態でやってます

赤文字のみ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")

 

f:id:yoshida931:20170212153358p:plain

 

#四つのグラフを2行2列のまとめてみます
par(mfrow = c(2,2))
mynorm(30,"n=30")
mynorm(100,"n=100")
mynorm(300,"n=300")
mynorm(1000,"n=1000")

f:id:yoshida931:20170212153327p:plain

回帰係数の区間推定

赤文字をRで実行

 

サンプル)統計学入門(基礎統計学),東京大学出版会,1991,p258

次のデータをコピーしてください
東京 福岡
1019.4 1018.4
1005.7 1007.6
1002 1006.2
1006.7 1009.9
1005.1 1010.8
1010.1 1013.2
1016.7 1016.2
1011 1009.1
999.5 1003.1
1006.9 1012.5
1001.9 1006.4
1007.5 1006.3
1014.4 1012.2
1014.3 1015
1014.6 1017.4
1009 1016.5
1006.7 1012.1
1009.4 1008.7
1011.8 1009.2
1009.4 1009.2

x<-read.table("clipboard",header=T)

 

回帰分析
lm(x$東京~x$福岡)

 

Call:
lm(formula = x$東京 ~ x$福岡)

Coefficients:
(Intercept)       x$福岡 
    16.0876       0.9822 

標本回帰方程式
Y=16.09+0.98X

 

 

定義
yi   実測値   x$東京
xi   実測値   x$福岡
B 回帰係数(最小二乗法より) この信頼区間を求める
β 回帰係数の不偏推定量 
ei 誤差項
Ei   (回帰)残差
Txx    xの偏差平方和      ∑(xi-xi平均値)^2
Txy xyの偏差平方和    ∑(yi-yi平均値)(xi-xi平均値)
Tyy    yの偏差平方和      ∑(yi-yi平均値)^2
σ二乗 誤差の分散
s二乗    誤差の分散の不偏推定量

 


回帰残差
母回帰方程式 yi=α+βxi+Ei  Ei~(0,σ二乗)  (i=1,2,...,n)
回帰残差Ei=yi-α-βxi
Ei<-(x$東京-16.09-0.9822*x$福岡)

 

誤差の分散の不偏推定量
回帰残差の平方和を(n-2)で割ることで回帰方程式の当てはまりの良さがわかります
s二乗<-(sum(Ei^2))/(n-2)
s二乗<-(sum(Ei^2))/18

 

推定値の標準誤差
s<-sqrt(s二乗)

 

ほとんどのeiは0±2sの中に入入ります
残差の散布図で確認(言訳;グラフは初心者です)
plot(ei,type="p",ylim=c(-6,6))
abline(h=s, col='black' ,lty="twodash")
abline(h=-s, col='black' ,lty="twodash")
abline(h=2*s, col='black' ,lty="dashed")
abline(h=-2*s, col='black',lty="dashed")


graphics.off()

偏回帰係数の標本分布
σは未知なのでsに置き換えます
Z=(β-B)÷σ/√Txx ~N(0,σ二乗)
t=(β-B)÷s/√Txx ~t(n-2)

 

Bの信頼区間
β-t×s/√Txx ≦ B ≦ β+t×s/√Txx
Txx<-sum( (x$福岡-mean(x$福岡) )^2)
0.98-2.101*s/sqrt(Txx)
0.6189584
0.98+2.101*s/sqrt(Txx)
1.341042
     0.62 ≦ B ≦ 1.34

 

 

Rの関数で確認します
confint(lm(x$東京~x$福岡),level=0.95)

              
(Intercept) -348.9125481 381.087689
x$福岡         0.6211872   1.343239

 

単回帰分析

赤文字のみRで実行

 

線形回帰分析
説明変数と目的変数を直線関係で傾向を示す。
説明変数と目的変数との関係を直線でモデル化する回帰分析。

非線形回帰分析
非線形関係でモデル化する回帰分析。
説明変数と目的変数を非直線的関係で傾向を示す。

 

ここでは単回帰分析について解説します
以下のように定義します
観測値 ( xi , yi )
観測値の平均値 ( ya , xa )
予測式 y'= b0 + b1*xi

最小二乗法を用いて回帰係数を求めます
残差の二乗和(残差平方和)を最小にする回帰係数をもとめます

残差ei ( yi - y' )
残差平方和 Q( b0 , b1 )

Q(b0,b1) = ∑(yi-y')^2 = ∑{yi-(b0+b1*xi)}^2    最小二乗基準

Q(b0,b1) を最小にするb0とb1を求めることで回帰直線を得る

Q(b0,b1) はパラメータb0b1の二次関数で最小値が存在します
そこでb0b1偏微分し、0に等しいとします

 

f:id:yoshida931:20170317200530p:plain
Y=b0+b1*x1
この式をyのx上への回帰方程式、あるいは回帰直線といいます

b0 : y切片
b1 : 偏回帰係数

 

例)
統計学入門 ,東京大学出版会,p42,1991より
yi<-c(114,124,143,158,166)
xi<-c(35,45,55,65,75)

 

平均を求めます
(ya<-mean(yi))
141
(xa<-mean(xi))
55

 

偏回帰係数
B<-sum( (yi-ya)*(xi-xa) )/sum( (xi-xa)^2)
1.38

 

Y切片
A<-ya-1.38*55
65.1

y=65.1+1.38x

 

平方和の分解
次に平方和を分解してみます
観測値の平方和(または総平方和)をSyとします
Sy = ∑(yi-ya)^2
   = ∑{(yi-y')+(y'-ya)}^2
   = ∑(yi-y')^2+∑(y'-ya)^2+∑(yi-y')*(y'-ya)
となります
ここで∑(yi-y')*(y'-ya)を考えてみます
予測値と残差の相関は0になります(証明は省略)
ry'ei=Sy'ei/Sy'*Sei=0
Sy'ei=∑(yi-y')*(y'-ya)=0となります

したがって
     Sy     = ∑(yi-y')^2 + ∑(y'-ya)^2
     Sy     =    Sei      +    SR
観測値の平方和 =  残差平方和   +  回帰による平方和

 

決定係数
回帰直線で説明できる部分の割合
R^2 = SR / Sei
Se=0,R^2=1 完全に説明できる
SR=0,R^2=0 全く説明できない

さきほどの例の決定係数を求めてみます

決定係数<-sum( (65.1+1.38*xi-ya)^2)/sum( (yi-ya)^2)
0.9836777

 

Rの関数を使って確認してみます
summary(lm(formula=yi~xi) )

 

Call:
lm(formula = yi ~ xi)

Residuals:
   1    2    3    4    5
 0.6 -3.2  2.0  3.2 -2.6

Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept)  65.1000     5.8284   11.17 0.001538 **
xi            1.3800     0.1026   13.45 0.000889 ***
---
Signif. codes: 
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.246 on 3 degrees of freedom
Multiple R-squared:  0.9837,    Adjusted R-squared:  0.9782
F-statistic: 180.8 on 1 and 3 DF,  p-value: 0.0008894

 

散布図と回帰直線
plot(xi,yi)
abline(lm(yi~xi))

f:id:yoshida931:20170613173429p:plain

参考) 豊田秀樹 (著, 編集);回帰分析入門 (Rで学ぶ最新データ解析) ,東京図書 ,2012