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

読者です 読者をやめる 読者になる 読者になる

統計学備忘録 since2016

Rを使用した統計学

異なる周期の三角関数の合成いろいろ

フーリエ変換の予備知識として
色々な三角関数の合成を可視化しておきます
あくまで準備段階ですので…

弧度法を度数法に変換します
x<-pi/180
y<-0:360
x<-x*y 

 

sin(x)、sin(2x) の合成
plot(sin(x)+sin(2*x),type="l",ylim = c(-2,2),ann = F)
lines(sin(x),type="l",col=2)  #赤
lines(sin(2*x),type="l",col=3) #緑

f:id:yoshida931:20170515150827p:plain


sin(x)、sin(2x)、sin(3x)の合成

plot(sin(x)+sin(2*x)+sin(3*x),type="l",ylim = c(-2.5,2.5),ann = F)
lines(sin(x),type="l",col=2)  #赤
lines(sin(2*x),type="l",col=3)  #緑
lines(sin(3*x),type="l",col=4)  #青

f:id:yoshida931:20170515150740p:plain


cos(x)、sin(3x)、-0.3cos(2x)の合成

plot(cos(x)+sin(3*x)-0.3*cos(2*x),type="l",ylim = c(-2.5,2.5),ann = F)
lines(cos(x),type="l",col=2)  #赤
lines(sin(3*x),type="l",col=3)  #緑
lines(-0.3*cos(2*x),type="l",col=4)  #青

f:id:yoshida931:20170515150715p:plain

 

 

cos(nx)とsin(nx)の合成

グラフでイメージするa*cos(x)とb*sin(x)の合成
ベクトル(振幅)は√(a^2+b^2)
sun(nx)、cos(nx)を合成すると位相が変化するが、周期は変化なし


振幅 √(2^2+2^2)=2.828427 #sqrt(8)
1周期
plot(2*cos(y*x)+2*sin(y*x),type="l",ylim = c(-3,3),ann = F)
lines(2*cos(y*x),type="l",col=2) #赤
lines(2*sin(y*x),type="l",col=3) #緑

f:id:yoshida931:20170515145631p:plain


2周期
plot(2*cos(2*y*x)+2*sin(2*y*x),type="l",ylim = c(-3,3),ann = F)
lines(2*cos(2*y*x),type="l",col=2) #赤
lines(2*sin(2*y*x),type="l",col=3) #緑

f:id:yoshida931:20170515145507p:plain

         振幅は変化なし、周期のみが変化

 


振幅√(2^2+5^2)=5.385165 #sqrt(29)
2周期
plot(2*cos(2*y*x)+5*sin(2*y*x),type="l",ylim = c(-5.5,5.5),ann = F)
lines(2*cos(2*y*x),type="l",col=2)
lines(5*sin(2*y*x),type="l",col=3)

f:id:yoshida931:20170515145436p:plain

          周期は変化なし、振幅のみが変化

三角関数のグラフ

角度の概念を復習します

弧度法
ラジアン(弧度)とは、「弧の長さ=半径」となる角度を「1」とする角度の単位のこと。
単位記号は「rad」
R言語ではこの「弧度法」を使用しています.

度数法
円周を360等分したものを1とする角度の単位のこと。
単位記号は「°」

つまり2πr÷r が360度に値します
2πが360度

グラフを理解するために弧度法を度数法に変換します
x<-pi/180
90度はx*90 となります。

sin(90度)=sin(90*x)=sin(pi/2)=1


したがって0°から360°までは次のようになります
y=0:360 
x*y=0°から360°まで

0°から360°までのsin波(黒),cos波(赤),sin波×cos波(緑),sin波×cos波(青)のグラフを重ねて書いてみます
plot(sin(y*x),type="l",ylim = c(-1.5,1.5),ann = F)    #タイトルを全て消去 ann = F
lines(cos(y*x),type="l",col=2)
lines(cos(y*x)+sin(y*x),type="l",col=3)
lines(cos(y*x)*sin(y*x),type="l",col=4)
na<-c("sin波","cos波","sin波+cos波","sin波×cos波")  #波に名前をつけます
legend(230,1.5,bty="n",na,text.col=1:4)       #bty="n" 凡例の枠無し
abline(h = 0,lty=2)                 #y軸のライン (破線)

 

f:id:yoshida931:20170512162035p:plain

RでBland Altman Plotを描く

 

データ
国際医療福祉大学 下井研究室
信頼性の検討方法:相対信頼性と絶対信頼性
例1(http://shimoi.iuhw.ac.jp/reliability_3_systematicbias.html
を使用させていただきました

x1
25.1 23.3 33.2 20.9 22.8 34.9 26.5 26.1 29.4 29.8 29.2 26.5 14.3 14.6 17.3 20.2 17.4 19.6 30.2 30.3 23.1 17.7 18.9 28.0 15.2 16.8 11.7 10.5 25.2 15.9

y1
17.6 21.5 25.6 10.5 16.4 31.2 20.8 27.0 32.8 30.9 22.9 24.0 18.2 14.9 13.2 17.2 19.5 15.8 21.2 26.1 14.1 14.8 12.3 17.6 14.3 8.2 11.3 11.8 23.4 10.6

 

まずはインストールを忘れずに

install.packages("BlandAltmanLeh")
library(BlandAltmanLeh)

 

ba <- bland.altman.stats(x1, y1)  
str(ba)            #データを概観して必要な部分だけ抜き取ります
plot(ba$means, ba$diffs,ylim=c(-8,15))
abline(h=ba$mean.diffs, lty=2, col=2)
abline(h=ba$CI.lines[2], lty=2) #誤差の許容範囲として最もあまい(optimisticな)範囲
abline(h=ba$CI.lines[5], lty=2) #誤差の許容範囲として最もあまい(optimisticな)範囲

f:id:yoshida931:20170509192237p:plain

 

bland.altman.plot(x1,y1)f:id:yoshida931:20170509192812p:plain

bland.altman.plot(x1,y1,conf.int = .95)

f:id:yoshida931:20170509192854p:plain

 

 

 

1標本t検定のスクリプト

只今練習中~

2標本xとyについて
sumxy <- function(x,y,sxy)  
{
xs<-summary(x)
ys<-summary(y)
xt<-t.test(x)
yt<-t.test(y)
print(xs);print(ys);print(xt);print(yt)
}
sumxy (x,y,sxy)
より

 

x<-c( );y<-c( )にデータを入れると、
各標本の要約・t値・p値・95%CIのみを出力させるプログラム

 

sumxy <- function(x,y,sxy)  

{
cat1<-c("Min.","1st Qu","Median","Mean", "3rd Qu.","Max.")
xs<-summary(x)
ys<-summary(y)
cat(cat1,"\n","x",xs,"\n","y",ys,"\n\n",file = "TEST001")  #各標本のまとめ

tx<-t.test(x)
cat2<-c("tvalue","pvalue","95%CI")
txt<- tx$statistic
txp<- tx$p.value
txc1<- tx$conf.int[1]
txc2<- tx$conf.int[2]
x1 <- sprintf("%.3f\t%.3f\t%.3f\t%.3f",txt,txp,txc1,txc2)

ty<-t.test(y)
tyt<- ty$statistic
typ<- ty$p.value
tyc1<- ty$conf.int[1]
tyc2<- ty$conf.int[2]
y1 <- sprintf("%.3f\t%.3f\t%.3f\t%.3f",tyt,typ,tyc1,tyc2)
cat(cat2,"\n","x",x1,"\n","y",y1,"\n",file = "TEST001", append =T)  

}
sumxy (x,y,sxy)

繰り返し for

いつも忘れるので書いておこう
#x=3に2を3回足す式
#3+2+2+2=9

myf<-function(){
x<-3
for(i in 1:3){
x<-x+2}
return(x)
}
myf()

 

#Xに2を3回加える式
#x+2+2+2
#xに5を代入すると5+2+2+2=11


myf<-function(x){
a<-x
for(i in 1:3){
a<-a+2}
return(a)
}
myf(5)


#Xに3を掛けることを3回繰り返す
#X*3*3*3
#xに5を代入すると
5*3*3*3

myf<-function(x){
for(i in 1:3){
x<-x*3}
return(x)
}
myf(5)

 

#10にxを3回加える式
#10+x+x+x
#x=2のとき10+2+2+2=16

myf<-function(x){
a<-10
for(i in 1:3){
a<-a+x}
return(a)
}
myf(2)

 

級内相関係数 ICC(1,1), ICC(1,k)

1人の検者により検査結果のばらつきを求める

ICC(1,1)

= ( (被験者間の平均平方-被験者内の平均平方)÷繰り返し数 ) ÷

 ( ( (被験者間の平均平方-被験者内の平均平方)÷繰り返し数) - 被験者内の平均平方)

 

被験者=水準
被験者間の平均平方
 =∑(被験者平均―総平均) ^2 ÷ (被験者数-1)

被験者内の平均平方

     =∑∑(X―被験者平均) ^2 ÷被験者数×(繰り返し数-1)

 

「一元配置分散分析 (対応なし) F値の算出方法」より
http://yoshida931.hatenablog.com/entry/2017/02/19/144239

 

仮想データ

一人の理学療法士が被験者A~Dに対してそれぞれ5回づつ検査を実施した結果

f:id:yoshida931:20170313152805p:plain

Rで算出してみます
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)
x<-data.frame(A,B,C,D)

パッケージirrのダウンロードが必要です
library(irr)
icc(t(x))   #行に各被験者のデータが並ぶように設定します

Single Score Intraclass Correlation
Model: oneway
Type : consistency
Subjects = 4
Raters = 5
ICC(1) = 0.24
F-Test, H0: r0 = 0 ; H1: r0 > 0
F(3,16) = 2.58 , p = 0.0898
95%-Confidence Interval for ICC Population Values:
-0.079 < ICC < 0.877

 

ICCの値と95%信頼区間の求め方
被験者数4、繰り返し数5
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)
ALL<-c(A,B,C,D)  #全データ 
mean(ALL)     #総平均


被験者間の平方和
y<-5*{(mean(A)-mean(ALL))^2+
(mean(B)-mean(ALL))^2+
(mean(C)-mean(ALL))^2+
(mean(D)-mean(ALL))^2}
=66.15

被験者内の平均平方
水準間平均平方=∑(水準平均―総平均) ^2 ÷ (水準数-1)
66.15÷(4-1) = 22.05


被験者内の平方和
sum( (A-mean(A))^2)+
sum( (B-mean(B))^2)+
sum( (C-mean(C))^2)+
sum( (D-mean(D))^2)
=136.8

被験者内
平均平方=∑∑(X―水準平均) ^2 ÷水準数×(n-1)
136.8÷{4☓(5-1)} = 8.55

 

ICC(1,1) 
= ( (被験者間の平均平方-被験者内の平均平方)÷繰り返し数 ) ÷
( ( (被験者間の平均平方-被験者内の平均平方)÷繰り返し数) - 被験者内の平均平方)
= ( (22.05-8.55)/5)/( ( (22.05-8.55)/5)+ 8.55)
=0.24

 

ICC(1,1)の95%信頼区間の求め方
被験者間の平均平方÷被験者内の平均平方
= 2.578947
qf(0.975,16,3)
= 14.23152
qf(0.975,3,16)
= 4.076823
( (2.578947/4.076823)-1)/ {(2.578947/4.076823)+(5-1)}
= -0.07931044
( (2.578947*14.23152)-1)/( (2.578947*14.23152)+(5-1))
=0.877157
95CI = [ -0.07931044 0.877157 ]

 


1人の検者により検査の信頼性を調べる
ICC(1,k)
= (被験者間の平均平方-被験者内の平均平方)÷ 被験者間の平均平方

 

Rで算出してみます
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)
x<-data.frame(A,B,C,D)
icc(t(x),"oneway","consistency", "average")

Average Score Intraclass Correlation
Model: oneway
Type : consistency (整合性、一貫性)
Subjects = 4
Raters = 5
ICC(5) = 0.612

F-Test, H0: r0 = 0 ; H1: r0 > 0
F(3,16) = 2.58 , p = 0.0898

95%-Confidence Interval for ICC Population Values:
-0.581 < ICC < 0.973

 

ICCの値と95%信頼区間の求め方
ICC(1,5)
= (被験者間の平均平方-被験者内の平均平方)÷ 被験者間の平均平方
= (22.05-8.55) / 22.05
= 0.6122449

ICC(1,5)の95%信頼区間の求め方
1-(1/(2.578947/4.076823))= -0.5808091
1-(1/(2.578947*14.23152))= 0.9727538
95CI = [-0.5808091 0.9727538]

 


次のような関数で全てのICCが算出可能になります
パッケージpsychを使用します
library(psych)
ICC(t(x))

Call: ICC(x = t(x))
Intraclass correlation coefficients
type ICC F df1 df2 p lower bound upper bound
Single_raters_absolute ICC1 0.24 2.6 3 16 0.09 -0.079 0.88
Single_random_raters ICC2 0.24 2.6 3 12 0.10 -0.084 0.88
Single_fixed_raters ICC3 0.24 2.6 3 12 0.10 -0.094 0.88
Average_raters_absolute ICC1k 0.61 2.6 3 16 0.09 -0.581 0.97
Average_random_raters ICC2k 0.61 2.6 3 12 0.10 -0.631 0.97
Average_fixed_raters ICC3k 0.61 2.6 3 12 0.10 -0.752 0.97

Number of subjects = 4 Number of Judges = 5

 

f:id:yoshida931:20170313152828p:plain

今回使用した結果は
Single_raters_absolute、Average_raters_absoluteの部分になります

 

対馬栄輝;頼性指標としての級内相関係数

http://www.hs.hirosaki-u.ac.jp/~pteiki/research/stat/icc.pdf