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

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

since2016 ときどきTEXのメモ

基準値の変更

df <- data.frame(c(rep("あり",20),rep("なし",10)),
           c(rep("あり",9),rep("なし",6),rep("あり",5),rep("なし",10)),
           c(30, 35, 29, 22, 32, 29, 30, 33, 35, 25, 
             30, 35, 29, 22, 32, 29, 30, 31, 25, 25, 
             16, 16, 18, 22, 10, 29, 20, 18, 12, 16))
colnames(df) <- c("筋トレ","歩行練習","効果")
library(tableone)
vars <- c("筋トレ", "歩行練習", "効果")
facvars <- c("筋トレ", "歩行練習")
table <- CreateTableOne(
    vars = vars, 
    #strata =
    data = df, 
    factorVars = facvars) 
print(table,quote = TRUE)

筋トレの基準を「なし」に変更

df$筋トレ <- as.factor(df$筋トレ)
df$筋トレ <- relevel(df$筋トレ, ref="なし")

vars <- c("筋トレ", "歩行練習", "効果")
facvars <- c("筋トレ", "歩行練習")
table <- CreateTableOne(
    vars = vars, 
    #strata =
    data = df, 
    factorVars = facvars) 
print(table,quote = TRUE)

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}

stripchartの横にエラーバー

サンプル

ID <- c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20)
治療 <- c("A","B","A","B","B","A","A","A","B","B","A","B","A","A","A","B","A","B","B","B")
BP変化量 <- c(治療後BP - 治療前BP)
dat1 <- data.frame(id, 治療, 治療前BP, 治療後BP, BP変化量)

標準偏差を挿入

stripchart(BP変化量~治療, data = dat1, method="jitter", 
    jitter = 0.1, vertical=TRUE, pch=20, cex=1.5, xlim=c(0.5,2.5), ylim=c(-20, -0),
    xlab="治療", ylab="BP変化量")

par(new=T) #図を追加

Aの平均 = meanAB[[1]] #ここに平均を代入します
Bの平均 = meanAB[[2]]

Aの標準偏差 = sdAB[[1]]
Bの標準偏差= sdAB[[2]]

dplot <- plot(c(1.15,2.15), c(Aの平均, Bの平均), ylim=c(-20, -0), xlim = c(0.5, 2.5), xaxt="n",xlab="", ylab="")

arrows(1.15:2.15, c(Aの平均, Bの平均) - c(Aの標準偏差, Bの標準偏差), 1.15:2.15,  c(Aの平均, Bの平均) + c(Aの標準偏差, Bの標準偏差), lwd=1.0,angle=90,length=0.1,code=3)

f:id:yoshida931:20201218162617p:plain:w400

標準誤差を挿入

stripchart(BP変化量~治療, data = dat1, method="jitter", 
    jitter = 0.1, vertical=TRUE, pch=20, cex=1.5, xlim=c(0.5,2.5), ylim=c(-20, -0),
    xlab="治療", ylab="BP変化量")

par(new=T) #図を追加

Aの平均 = meanAB[[1]] #ここに平均を代入します
Bの平均 = meanAB[[2]]

dA <- dim(subset(dat1, dat1$治療=="A"))
dB <- dim(subset(dat1, dat1$治療=="B"))
Aの標準誤差 = sdAB[[1]]/sqrt(dA[1])
Bの標準誤差 = sdAB[[2]]/sqrt(dB[1])

dplot <- plot(c(1.15,2.15), c(Aの平均, Bの平均), ylim=c(-20, -0), xlim = c(0.5, 2.5), xaxt="n",xlab="", ylab="")

arrows(1.15:2.15, c(Aの平均, Bの平均) - c(Aの標準誤差, Bの標準誤差), 1.15:2.15,  c(Aの平均, Bの平均) + c(Aの標準誤差, Bの標準誤差), lwd=1.0,angle=90,length=0.1,code=3)

f:id:yoshida931:20201218163524p:plain:w400

参考
yoshida931.hatenablog.com

回帰分析の結果の解釈

単回帰分析

dat <- LifeCycleSavings      

# sr(個人貯蓄の合計), pop15(15歳未満の人口の数値%), pop75(75歳未満の人口の数値%), dpi(可処分所得、自由に使える収入), ddpi(dpiの成長率)

reg1 <- lm(sr ~ pop15, data=dat)    
summary(reg1) 

次のような結果が出力されます

            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 17.49660    2.27972   7.675 6.85e-10 ***
pop15       -0.22302    0.06291  -3.545 0.000887 ***

結果の解釈について
Intercept=17.4966
「15歳未満の人口が0%の場合には個人貯蓄合計が平均17.49660になる」と解釈できるのですが、15歳未満の人口が0%というのはあまり現実的な解釈ではありませんので、次のように書き換えます.

reg2 <- lm(sr ~ I(pop15-10), data=dat)
summary(reg2)

結果の切片は、15歳未満の人口が10%の場合の平均的な個人貯蓄合計を意味することになります。切片が変わるだけすのでpop15の係数の推定値は変化していません。

              Estimate Std. Error t value Pr(>|t|)    
(Intercept)   15.26642    1.67802   9.098 5.09e-12 ***
I(pop15 - 10) -0.22302    0.06291  -3.545 0.000887 ***

結果の解釈について
Intercept=15.26642
15歳未満の人口が10%の場合には個人貯蓄の合計が約15.3になり、さらに15歳未満の人口が1%増加するごとに個人貯蓄の合計が平均して0.22減少することが推定されます。
<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<

Residual standard error: 4.03 on 48 degrees of freedom
Multiple R-squared:  0.2075,    Adjusted R-squared:  0.191 
F-statistic: 12.57 on 1 and 48 DF,  p-value: 0.0008866

f:id:yoshida931:20201216181235p:plain

自由度 n=50, p=1 (単回帰なのでp=1となる)
Residual standard error(残差標準誤差)…でも、残差の「標準偏差」を意味している
  \sqrt{\frac{SSE}{自由度}}
Multiple R-squared (決定係数)
  R2= \frac{SSR}{SST}\

決定係数は説明変数の数が増えるほど1に近づくので自由度で調整
Adjusted R-squared(自由度調整済み決定係数)
  AR2=1- \frac{\frac{SSE}{n-p-1}}{\frac{SST}{n-1}}

F-statistic: 13.92 on 2 and 12 DF, p-value: 0.0007471
  MSR <- SSR/p
  MSE <- SSE/(n-p-1)
  F <- MSR/MSE    p-value=1-pf(F, p, n-p-1) ---> F分布から算出

分散分析表から求めることも可能

anova(reg2)

              Df Sum Sq Mean Sq F value    Pr(>F)    
I(pop15 - 10)  1 204.12  204.12  12.569 0.0008866 ***
Residuals     48 779.51   16.24          

SSR =  204.12     
MSR =  204.12 / 1 =  204.12   
SSE =  779.51
MSE =  779.51 / (50-1-1) = 16.24

回帰式から直接求める方法
注) reg1の係数を使用します

b0 <- reg1$coeff[1] #切片
b1 <- reg1$coeff[2] #傾き
yhat <- b0 + b1*dat$pop15 #推定値(srの推定値)
heikin <- mean(dat$sr) #=推定値の平均
SSE <- sum((dat$sr - yhat)^2)
SSR <- sum((yhat - heikin)^2)
SST <- sum((dat$sr - heikin)^2)

カッパ係数

投稿日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は同一被験者のデータになるので、「対応がある」という表現が使用されているようです.

まとめ

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