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

統計学備忘録(R言語のメモ)

since2016 ときどきTEXのメモ

トービット回帰直線

トービット回帰直線

投稿日2017.8.4

打ち切りデータの場合、つまり天井効果や床効果が生じているデータの場合には、トービット回帰直線で分析します.

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

 

Rのサンプルanscombeを使います.( n=11の4セット )

x1<-anscombe$x1

x4<-anscombe$x4
y1<-anscombe$y1
y4<-anscombe$y4

par(mfrow = c(1,2))
plot(x1,y1,xlab = "図A",ylab = "")
plot(x4,y4,xlab = "図B",ylab = "")

f:id:yoshida931:20170804195919p:plain



dev.off()

図A 通常の散布図
lm1<-lm(y1~ x1,data= anscombe)
plot(x1,y1)
abline(lm1$coefficients[1],lm1$coefficients[2])

f:id:yoshida931:20170804195838p:plain

視覚的にも何となくこの直線でフィットしているように見えます
x軸を予測値にして、残差を標準化して更に見やすくしてみます
ei<- y1-( lm1$coefficients[1]+ lm1$coefficients[2]* x1)
eis<-ei/sqrt( (sum(ei^2) )/(11 - 2) )  #標準化残差
plot(lm1$coefficients[1]+ lm1$coefficients[2]* x1, eis,xlab="図A(x軸:予測値)",ylab=" ", xaxt="n")
abline(h=0, lty=2, col=2)

f:id:yoshida931:20170804200055p:plain
図B 打ち切りデータがある場合
トービット回帰直線( 関数vglmを使用 )
通常の回帰直線(実践)
y4x4 <- lm(y4 ~ x4); summary(y4x4)
plot(x4,y4)
abline(summary(y4x4)$coefficients1,summary(y4x4)$coefficients2)
関数vglmを使用してトービット回帰直線を破線で加えます
y4x4t <- vglm(y4 ~ x4, tobit(Lower=8, Upper=19), trace=TRUE); summary(y4x4t)
abline(summary(y4x4t)@coefficients1,summary(y4x4t)@coefficients3,lty=2)

f:id:yoshida931:20170804200433p:plain

床効果で下に引っ張られていた回帰直線をやや上方修正できました.

ユールのQ

ユール ( Yule ) のQ
ユールの関連係数、ユールの連関係数とも呼ばれています.関連が強いほど1または-1に近い値をとります。

f:id:yoshida931:20170718215027p:plain

Q= ( a*d - b*c ) / ( a*d + b*c ) =  ( オッズ比 - 1 ) / ( オッズ比 + 1 )

例)2つの質問 ( q1 , q2 ) に関する答え ( yes , no ) のベクトルを考えます.
q1<-c("yes","yes","yes","yes","yes","yes","yes","no","no","no")
q2<-c("yes","yes","yes","no","no","no","no","yes","yes","no")

tq<- table ( q1,q2 )    

f:id:yoshida931:20170718172336p:plain

Q = ( 1*3 - 2*4 ) / ( 1*3 + 2*4 )
    = - 0.4545455

Rの関数で確認してみます
install.packages("psych")
library(psych)
Yule(tq)
- 0.4545455

 

Rstudioでの data.frame 表示の不具合

データフレームを表示したときに、ラベルが見えなくなりました

f:id:yoshida931:20170720121311p:plain

 

RStudio-1.0.143 から RStudio-1.0.153 へバージョンアップすることで不具合は解消しました.
https://www.rstudio.com/products/rstudio/download/preview/

 

f:id:yoshida931:20170720121437p:plain

 

Rstudioのサポートページの掲示板で助言を受けました!

https://www.rstudio.com/

f:id:yoshida931:20170720121521p:plain

ファイ係数

連関係数
 クロス集計表における2つの変数間の関連性の程度を表す指標として,連関係数が提案されています.

連関係数<|0.2| 連関はほとんどない
|0.2|≦連関係数<|0.4| やや連関がある
|0.4|≦連関係数<|0.7| 強い連関がある
|0.7|≦連関係数 かなり強い連関がある
   対馬栄輝; よくわかる医療統計 -「なぜ?」にこたえる道しるべ, 東京図書, 2015

 

ファイ係数 (φ係数)
 クロス集計表における行要素と列要素の関連の強さを示す指標. 1と0の2つの値からなる変数に対して計算される相関係数です。つまり、2つの2値変数をそれぞれ1と0で数値化した上で算出した積率相関係数です。1と0の2つの値とは、「好き、嫌い」、「はい、いいえ」など、2択で表現されるデータです. 選択肢Aに0を、選択肢Bに1を割り付けることで、ファイ係数を求めることができます.


例)2つの質問 ( q1 , q2 ) に関する答え ( yes , no ) のベクトルを考えます.
q1<-c("yes","yes","yes","yes","yes","yes","yes","no","no","no")
q2<-c("yes","yes","yes","no","no","no","no","yes","yes","no")
q101<-ifelse(q1=="yes",1,0)        # 0,1のベクトルに変更
  q101 <- c(1, 1, 1, 1, 1, 1, 1, 0, 0, 0)
q201<-ifelse(q2=="yes",1,0)       # 0,1のベクトルに変更
  q201 <- c(1, 1, 1, 0, 0, 0, 0, 1, 1, 0)


ピアソンの積率相関係数を2行×2列のクロス集計表に適用したもの
   ピアソンの積率相関係数

           f:id:yoshida931:20170718173706p:plain

mean ( q101 )
=0.7
mean ( q201 )
=0.5
sum ( ( q101-0.7 ) * ( q201 - 0.5 ) )
= - 0.5
sum ( ( q101 - 0.7 )^2 )
=2.1
sum ( ( q201 - 0.5 )^2 )
=2.5
r = - 0.5 / sqrt ( 2.1 * 2.5 )
  = - 0.2182179


Rの関数で確認してみます

cor(q101,q201)    
           
#相関係数をcor( )で算出
= - 0.2182179


次に例題をテーブルにして眺めてみます
tq<- table ( q1,q2 )
 
   

f:id:yoshida931:20170718172336p:plain

t<- tq [ c(2,1) , c(2,1) ]   
#これは不要ですがラベルの並べ替えです.yesから始めたいので挿入しておきます(1番目と2番目の入れ替え)

f:id:yoshida931:20170718172305p:plain

理解のためにデータフレームに変換します
t< - data.frame ( t ) 
     

f:id:yoshida931:20170720121016p:plain

 
     q1   q2    Freq
1   yes  yes   3         #  t $ Freq[1]
2   no    yes   2        #  t $ Freq[2]
3   yes   no    4        #  t $ Freq[3]
4   no    no     1        #  t $ Freq[4]

ファイ係数φその1 カイ二乗値より
φ=√(カイ二乗値/n)
chisq.test ( tq , correct=F )
X-squared = 0.47619
φ=sqrt ( 0.47619 / 10 )
= 0.2182178

ファイ係数φその2  観測値と周辺度数より

f:id:yoshida931:20170718213849p:plain

φ =  ( a*d - b*c ) / √ ( a+b )*( c+d )*( a+c )*( b+d )*
x <- t$Freq[1] * t$Freq[4] - t$Freq[2] * t$Freq[3]   # x = 対角線上の観測値を掛けたものの差
y <- ( t$Freq[1] + t$Freq[3] )*                                 # y = 全ての周辺度数を掛けたもの
  ( t$Freq[2] + t$Freq[4] )*
  ( t$Freq[1] + t$Freq[2] )*
  ( t$Freq[3] + t$Freq[4] )
φ = x / sqrt(y)
    = - 0.2182179

 

Rの関数で確認してみます
install.packages ( "psych" )
library ( psych ) 
phi ( tq , digits=8 )
= - 0.2182179