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

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

since2016 ときどきTEXのメモ

Tex 図を並べる

図を縦に二つ並べる

\begin{figure}[H]
  \centering
  \includegraphics[width=40mm]{ファイル01名.jpg} 
  \caption{図の説明01}
\end{figure}

\begin{figure}[H]
  \centering
  \includegraphics[width=40mm]{ファイル02名.jpg} 
  \caption{図の説明02}
\end{figure}

「図〇〇」をカットする方法

\begin{figure}[H]
  \centering
  \includegraphics[width=40mm]{ファイル01名.jpg} 
  \captionsetup{labelformat=empty,labelsep=none}\caption{図の説明01}
\end{figure}

横に並べる方法1

\begin{figure}[htbp]
  \begin{minipage}{0.5\hsize}
    \begin{center}
      \includegraphics[width=50mm]{ファイル01名.jpg}
    \end{center}
      \caption{図01の説明}\label{fig:one}
  \end{minipage}   %注意 minipageは、ここで改行したら2行になる
  \begin{minipage}{0.5\hsize}
    \begin{center}
      \includegraphics[width=20mm]{ファイル02名.jpg}
    \end{center}
      \caption{図02の説明}\label{fig:two}
  \end{minipage}
\end{figure}

横に並べる方法2
2段組を利用
\usepackage{multicol} : 2列表示するパッケージ 何列でもOK
\setlength{\columnsep}{10mm} : 列の間隔の設定  

\begin{multicols}{2}  %2列
  \begin{figure}[H]
    \begin{center}
      \includegraphics[width=10mm]{ファイル01.jpg}
    \end{center}
      \captionsetup{labelformat=empty,labelsep=none}\caption{図01の説明}
  \end{figure}
        
  \begin{figure}[H]
    \begin{center}
      \includegraphics[width=10mm]{ファイル02.jpg}
    \end{center}
      \captionsetup{labelformat=empty,labelsep=none} \caption{図02の説明}
    \end{figure}
\end{multicols}

カッパ係数

投稿日2018.3.12  最終更新日2020.10.13
個人的によく利用させていただいております以下のHPをもとに、今回はカッパ係数について少し勉強してみます
統計学入門−第5章
まずはHPに掲載してある次のサンプルデータを使用して、Rを使って処理してみます

分類数が2つの場合

rater01<-c(rep(1,40),rep(2,40),rep(1,10),rep(2,10)) 
rater02<-c(rep(1,40),rep(2,40),rep(2,10),rep(1,10))
tk <- table(rater01,rater02)
addmargins(tk)

       rater02
rater01   1   2 Sum
    1    40  10  50
    2    10  40  50
    Sum  50  50 100

次にRの関数を使用して求めてみます

install.packages('vcd')
library(vcd)
Kappa(tk)

           value  ASE   z  Pr(>|z|)
Unweighted   0.6 0.08 7.5 6.382e-14
Weighted     0.6 0.08 7.5 6.382e-14

Z検定のp値を算出

library(irr) # irrの読み込み
dk <- cbind(rater01,rater02)
kappa2(dk)

 Cohen's Kappa for 2 Raters (Weights: unweighted)
 Subjects = 100 
   Raters = 2 
    Kappa = 0.6 
        z = 6 
  p-value = 1.97e-09 = (1-pnorm(6))*2

分散を求めて95%信頼区間を算出

n <- length(rater01)
po <- (40+40)/100
pc <- (50/100)*(50/100)+(50/100)*(50/100)
vk <- (po*(1-po))/(n*(1-pc)^2)

as.numeric(kappa2(dk)[5])-1.96*sqrt(vk) 
as.numeric(kappa2(dk)[5])+1.96*sqrt(vk) 

0.44320.7568

参考HPと比較してください

分類数が3つの場合

rater02<-c(rep(1,19),rep(2,26),rep(3,4),rep(1,17),rep(1,7),rep(2,7),rep(2,5),rep(3,3),rep(3,12)) 
rater03<-c(rep(1,19),rep(2,26),rep(3,4),rep(2,17),rep(3,7),rep(1,7),rep(3,5),rep(1,3),rep(2,12))
tk2 <- table(rater02,rater03)
addmargins(tk2)

       rater03
rater02   1   2   3 Sum
    1    19  17   7  43
    2     7  26   5  38
    3     3  12   4  19
    Sum  29  55  16 100

d2 <-cbind(rater02, rater03)
kappa2(d2)

 Cohen's Kappa for 2 Raters (Weights: unweighted)

 Subjects = 100 
   Raters = 2 
    Kappa = 0.198 

        z = 2.8 
  p-value = 0.00508 

#95%信頼区間
n2 <- length(rater03)
po2 <- (19+26+4)/100
pc2 <- (43/100)*(29/100)+(38/100)*(55/100)+(19/100)*(16/100)
vk2 <- (po2*(1-po2))/(n2*(1-pc2)^2)

as.numeric(kappa2(d2)[5])-1.96*sqrt(vk2) 
as.numeric(kappa2(d2)[5])+1.96*sqrt(vk2)  

0.04390565~0.3520686

重み付けカッパ係数

上記の例で、解答1、2、3を順序尺度として考えているので、rankA=一致とした場合に、次のような順序が考えられます
rankA > rankB > rankC
1と3のペアよりも1と2のペアや2と3のペアの方が一致に近いという考え方です

     [,1]    [,2]    [,3]   
[1,] "rankA" "rankB" "rankC"
[2,] "rankB" "rankA" "rankB"
[3,] "rankC" "rankB" "rankA"

この重み付けでカッパ係数をRで算出してみます

kappa2(d2, "squared")   #自乗距離に応じた重み付け

 Cohen's Kappa for 2 Raters (Weights: squared)

 Subjects = 100 
   Raters = 2 
    Kappa = 0.196 

        z = 2 
  p-value = 0.0453 


kappa2(d2, "equal")   #均等に重み付け

 Cohen's Kappa for 2 Raters (Weights: equal)   #評価者間の不一致を均等に重み付け

 Subjects = 100 
   Raters = 2 
    Kappa = 0.197 

        z = 2.67 
  p-value = 0.0076 

ANOVA君 で分散分析

投稿日2017.8.25, 更新日2020.8.17

ANOVA君はフリーの統計ソフトウェア「R」で動作する分散分析関数です。ダウンロード、使い方、リリース情報などは「井関龍太のページ」を必ずご確認ください.

ANOVA君 - 井関龍太のページ

このブログではanovakun version 4.8.5を使用させていただきます 。このテキストファイルにプログラムが記載してあり、マニュアルにもなります。

まずはテキストファイル(コード)を読み込みます

source("C:\\Users\\-----------\\anovakun_485.txt")

A~Zが要因、被験者間要因は「s」の左側,被験者内要因は「s」の右側に記載、数字は各要因の水準の数。入れ子の最も外側にある要因から順番に「A」~「Z」のラベルが対応。

一元配置分散分析と多重比較

n=30, 繰り返し数10の例

testA <- c(108, 105, 98, 112, 110, 109, 115, 120, 107, 114)
testB <- c(91, 93, 89, 109, 98, 88, 102, 110, 97, 109)
testC <- c(79, 90, 87, 87, 78, 75, 96, 109, 89, 90)
dat1 <- data.frame(testA , testB , testC)

#long型に変換
library(tidyr)
df <- gather(dat1, test, data) 

    test data
1  testA  108
2  testA  105
3  testA   98
4  testA  112
5  testA  110
6  testA  109
7  testA  115
8  testA  120

anovakun(df, "As", 3)
#被験者間要因は「s」の左側のA,3は各要因の水準の数(A, B, C)

多重比較:デフォルトはShafferの方法(Bonferroniの方法以外が可能)

反復測定と球面性検定とイプシロン修正

n=10の例

運動後1<- c(108, 105, 98, 112, 110, 109, 115, 120, 107, 114)
運動後2<- c(98, 98, 95, 98, 100, 100, 105, 100, 98, 109)
運動後3<- c(95, 87, 90, 84, 88, 89, 90, 89, 90, 99)
dat1 <- data.frame(運動後1,運動後2, 運動後3) #反復測定はこの型のまま使用する

   運動後1分 運動後2分 運動後31        108        98        95
2        105        98        87
3         98        95        90
4        112        98        84
5        110       100        88
6        109       100        89
7        115       105        90
8        120       100        89
9        107        98        90
10       114       109        99

anovakun(dat1, "sA" ,3 , mau = T)  
#被験者内要因は「s」の右側のA,3は各要因の水準の数(運動後1分, 2分, 3分)
#mau = T:Mauchlyの球面性検定(デフォルトは、Mendozaの多標本球面性検定)  

== Mauchly's Sphericity Test and Epsilons ==

------------------------------------------------------------------------
 Effect      W  approx.Chi  df      p         LB     GG     HF     CM 
------------------------------------------------------------------------
      A 0.4396      6.5750   2 0.0373 *   0.5000 0.6409 0.7008 0.6081 

分散分析はそれぞれの群が独立していることが前提となります。検定の結果(p値=0.0373)から球面性の仮定が棄却されるので、被験者内効果をイプシロン修正します。さらに一般的な手法としてMANOVA混合モデルがあります。

anovakun(dat1, "sA" ,3 , gg = T)

イプシロン調整済p値が算出される

== Adjusted by Greenhouse-Geisser's Epsilon ==

-------------------------------------------------------------
 Source         SS    df         MS  F-ratio  p-value      
-------------------------------------------------------------
      s   374.6667     9    41.6296                        
-------------------------------------------------------------
      A  1940.6000  1.28  1514.0480  68.0293   0.0000 ***  
  s x A   256.7333 11.54    22.2558                        
-------------------------------------------------------------     

もし球面性の仮定が成立(p>0.05)しているのであれば、独立ではないのですが一元配置分散分析の結果を利用することになります。

二元配置分散分析と交互作用

n=30の例

性別 <- c("男", "男", "男", "男", "女", "女", "女", "女", "女", "女")
A群 <- c(108, 105, 98, 112, 110, 109, 115, 120, 107, 114)
B群 <- c(91, 93, 89, 109, 98, 88, 102, 110, 97, 109)
C群 <- c(79, 90, 87, 87, 78, 75, 96, 109, 89, 90)
dat2 <- data.frame(性別, A群, B群, C群)

#long型に変換
library(tidyr)
df2 <- gather(dat2, test, data, -性別) 

   性別 test data
1    男  A群  108
2    男  A群  105
3    男  A群   98
4    男  A群  112
5    女  A群  110
6    女  A群  109
7    女  A群  115
8    女  A群  120
9    女  A群  107
10   女  A群  114

anovakun(df2, "ABs" ,2 ,3)

Rを使用して分散分析を実施するのであれば、ANOVA君をお勧めします。

分散分析(繰り返し vs 反復)

登校日:2017.11.06, 更新日:2020.08.12

「繰り返し」と「反復」は区別なく混同して使用される場合もあります.要点のみ記載しておきます。

繰り返し

実験内容としては「水準の設定を最初から施して既定のサンプル数を繰り返してとる」ことになります. 「繰り返し」とは、同じ実験を繰り返すことを意味しています.実験Aを被験5名 に実施する場合、5回「繰り返す」ことになります。

一元配置分散分析の例
検査A、検査B、検査Cに被験者4名を無作為に割り当て、3回(検査A, 検査B, 検査C)の検査結果に差があるかを検定する。 この例では、N=12、繰り返し数=4回(3回ではないことに注意)。

No <- c(1:4)
検査A <- c(44.8, 55.3, 55.4, 38.4)
検査B <- c(50.8, 52.1, 62.2, 65.1)
検査C <- c(35.7, 31.8, 43.0, 31.2)
x <- data.frame(No, 検査A, 検査B, 検査C)

x <- read.table("clipboard", header=T)
#long形式に変換
library(tidyr)
df <- x %>% gather(test, data, -No) 
stripchart(data~test, data=df, vertical = T, method="jitter", pch=16, xlim=c(0.5,3.5))

f:id:yoshida931:20200812175921p:plain:w350

反復

反復繰り返しが混同されて使われている場合も多いので注意が必要です。ここでは対応のある測定方法を反復と呼びます。

例)1症例における漸増運動負荷試験による血圧の変化

data <- c(126, 147, 159, 165)
plot(data, type="b", xaxt="n", xlim=c(0.5,4.5), ylim=c(120, 170), xlab="", ylab="")
name<-c("運動前","20W","40W","60W")
axis(side=1,at=c(1:4),labels=name) 

f:id:yoshida931:20200812182256p:plain:w400

二元配置分散分析 vs 反復測定分散分析(repeated ANOVA)

例)4名の学生を対象とした実験。実習1ヶ月前、実習直前、実習直後の唾液アミラーゼ値に差があるか調べました。この場合、要因を「個人差」と「実験時期」の2つとして考えて二元配置分散分析を適応するのは不適切です。分散分析は無作為独立標本が大前提となっています。この例ではデータ同士が対応して関連性があるため、球面性の検定結果が有意であれば分散共分散構造を考慮した反復測定分散分析で解析すべきです。 この例ではN=3、反復回数5回となります。球面性検定でp値<0.05のときにはイプシロン修正で対応し、p値<0.0001の場合には分散分析を多変量に拡張した多変量分散分析(MANOVA:multivariate analysis of variance)で対応します(統計学入門−第18章)。

No <- c(1:4)
実習1ヶ月前 <- c(44.8, 58.4, 55.3, 38.4)
実習直前 <- c(50.8, 52.1, 62.2, 65.1)
実習直後 <- c(9.7, 31.8, 43.0, 22.2)
x <- data.frame(No, 実習1ヶ月前, 実習直前, 実習直後)

df <- x %>% gather(time, outcome, -No)
with(df,{ x1  = factor(time, levels=c("実習1ヶ月前", "実習直前", "実習直後")) #並び替え
interaction.plot( x1, No, outcome, col=1:3, xlab="", ylab="", legend=F, lwd=2.5)})

f:id:yoshida931:20200812173039p:plain

おまけ(乱塊法 Random Block method
Fisherの3原則を満たす実験方法

# 無作為抽出した繰り返し数3回(検査α、β、γ)の例
# 1日目=被験者1,2,3、2日目=被験者4,5,6、3日目=被験者7,8,9

検査α 被験者6→被験者8→被験者9
検査β 被験者4→被験者1→被験者2
検査γ 被験者5→被験者3→被験者7

# 次のように書き換えると偏りが見えてきます
# 次に1日に3名実施して場合、1日を1ブロック(同じ条件で実験できる1領域)としたら・・・
# 偏りのある採取方法になっています

検査順序
1日目(β → β → γ)
2日目(β → γ → α)
3日目(γ → α → α)

# 検査日による偏りを解消
1日目 (α → β → γ)
2日目 (β → γ → α)
3日目 (β → α → γ)

     1日目 2日目 3日目
      α   β   β
検査要因  β   γ   α
      γ   α   γ
#「時期」をブロックとした分散分析(繰り返しなし)・・・方法は二元配置分散分析

無作為に抽出されたサンプルに対して、制御因子(再現性あり)ブロック因子(再現性なし)の2因子を用いる特殊な分散分析が乱塊法と呼ばれています.検査項目が制御因子、時期(1日目、2日目、3日目)がブロック因子となります.

対応のある、対応のない

繰り返しがある場合には、別々の被験者からのデータ収集になるため「対応がない」と表現されています.また繰り返しがない場合には、上述の例のように被験者AのA1、A2、A3、A4は同一被験者のデータになるので、「対応がある」という表現が使用されているようです.

まとめ

結局色々な表現方法があるので、データ間の独立性やブロックなどを理解した上で、研究目的と実験方法に適応した解析をしなければなりません。

単回帰分析を一括で実行する方法

tokei.net 全人類が分かる統計学より
to-kei.net

上記サイトに自分なりの注釈をつけました。
sampleもまったく同じです。

番号 <- c(1:30)
年齢 <- c(22,23,24,25,27,28,28,29,30,31,32,32,33,33,34,36,37,37,38,39,40,42,46,49,50,53,56,58,64,65)
血圧 <- c(110,128,104,112,108,126,126,104,125,120,116,124,106,134,128,128,116,132,134,116,120,130,126,140,156,124,118,144,142,144)
肺活量 <- c(110,128,104,112,108,126,126,104,125,120,116,124,106,134,128,128,116,132,134,116,120,130,126,140,156,124,118,144,142,144)
性別 <- c("M","M","F","F","M","F","F","F","F","F","M","M","F","F","M","M","M","M","F","M","F","F","M","F","M","M","M","M","F","F")
病気 <- c(1,1,0,0,0,0,1,1,1,1,1,0,0,0,1,1,1,1,0,1,1,1,0,0,1,1,1,0,1,0)
体重 <- c(79,65,53,45,80,50,43,55,47,49,64,61,48,41,70,55,70,90,39,86,50,65,58,45,60,71,62,51,40,42)
data <- data.frame(番号, 年齢, 血圧, 肺活量, 性別, 病気, 体重)

head(data, n=3)

dat <- data
names(dat)[6] <- "Y" # 6列目をYに変更

#まず関数を作成して、どんな解析の、どこの結果を記載するのか決めます
#func01は、関数名なのでなんでもOK、他は特に変更する必要なし

func01 <- function(dat,filename,response) 
{
  ans <- lm(dat$Y~.,data=dat)    #重回帰
  s.ans=summary(ans)     #重回帰の要約
  coe <- s.ans$coef    #重回帰の結果(係数、標準誤差、t値、p値)
  CI.low <- coe[,1]-1.96*coe[,2]    #係数の検定結果 95%信頼区間下限
  CI.up <- coe[,1]+1.96*coe[,2]    #係数の検定結果 95%信頼区間上限
  aic <- AIC(ans)     #AIC
  N <- nrow(dat)    #セットの行数
  result <- cbind(coe,CI.low,CI.up,aic,N)    #95%信頼区間とAICを合体
  result[2:nrow(result),c(7,8)]=""   #2行目以降の7、8列目を空白にする・・・?
  res <- paste("Y = ",response)    #単回帰の結果の1行目の文字列になる
  #ファイルの構成
  write.table(res,filename,append=T,quote=F,sep=",",row.names=F,col.names=F)
  write.table(matrix(c("",colnames(result)),nrow=1),filename,append=T,quote=F,sep=",",row.names=F,col.names=F)
  write.table(result,filename,append=T,quote=F,sep=",",row.names=T,col.names=F)
  write.table("",filename,append=T,quote=F,sep=",",row.names=F,col.names=F)
}

filename <- "てすとファイル1.csv"    #ファイル名を作成

df1 <- dat[,-1]    #番号をカット(入れてても問題なし)
for(i in 1:ncol(df1))    #列数分繰り返す(6回)
{
  if(i != 5) #5回目(5行目)を回避
  {
    #単回帰を繰り返すために2列のデータフレームを作成(ここでは5列目が目的変数)
    dat2 <- df1[ , c(5,i)]
    colnames(dat2)[1] = "Y"      #2列になったdatの1列目が目的変数
    func01(dat2, filename, "病気")    #各解析結果の1行目Y =の後に目的変数名「病気」を挿入
  }
}

#if(i != 5)
#5回目(5行目)を回避しなくてもやれる!

dat <- data
names(dat)[6] <- "Y"     # 6列目をYに変更
filename <- "てすとファイル2.csv"
df2 <- dat[,-1]     #IDをカット 入れてても大丈夫

for(i in 1:ncol(df2))    #列数分繰り返す(6回)
{
  dat3 <- df2[,c(5,i)]
  colnames(dat3)[1] = "Y" #datが2列
  func01(dat3, filename, "病気") #ここで目的変数を挿入
}

#同じ変数同士で、fitし過ぎているので注意されます
#essentially perfect fit!!!