ロジスティック回帰の基準(reference) 変更

ロジスティック解析を行うときには、オッズとオッズ比の基準 (reference) が非常に重要になります．特に日本語での記述の場合には、解析を行う前に基準を設定することをお勧めします．

df <- data.frame(c(rep("あり",15),rep("なし",15)),
               c(rep("有効",12),rep("無効",3),rep("有効",4),"無効","有効","有効",rep("無効",8)))
colnames(df) <- c("筋トレ","効果")
dat1 <- xtabs(~ 効果 + 筋トレ, data = df)  

dat1
      筋トレ
効果   あり なし
  無効    3    9
  有効   12    6

基準を確認します．オブジェクトの内容を情報付きで簡潔に表示 str( ) ．

str(dat1)

  ..$ 効果  : chr [1:2] "無効" "有効"
  ..$ 筋トレ: chr [1:2] "あり" "なし"

何も設定しない場合には、この先に書かれている方が基準として決定されます． 効果の基準は「無効」、筋トレの基準は「あり」．効果を目的変数、筋トレを説明変数としたロジスティック回帰を行います．

#データフレームを操作して、目的変数を0、1に置き換え名義変数にします
#0,1に置き換えていますが、基準は「0=無効」のままです
df$効果n <- ifelse(df$効果=="有効",1 ,0)
df$効果n <- as.factor(df$効果n)
glm(formula = 効果n ~ 筋トレ, family = binomial, data = df)

Coefficients:
(Intercept)   筋トレなし  
      1.386       -1.792  
#オッズ比
exp(-1.792) 
オッズ比＝0.1666266

オッズは効果が「無効」に対する「有効」の比率、オッズ比は 筋トレ「あり」のオッズに対する筋トレ「なし」のオッズの比率．つまり、筋トレしなかったら筋トレした場合と比べてオッズは0.167倍に下がることが推定できます．

変数を変換したり、層別を行ったりする過程で、何を基準にしているのか分からなくなる場合もあります．以下のように事前に基準を設定しておくことでエラーを防ぐことができます．

「筋トレなし」と「効果なし」を基準にする場合には以下のように設定します。

#データフレームから直接「効果なし」を基準に設定
df$効果n <- as.factor(df$効果n)
df$効果n <- relevel(df$効果n, ref="0") 
#データフレームから直接「筋トレ」なしを基準に設定
df$筋トレ <- as.factor(df$筋トレ)
df$筋トレ <- relevel(df$筋トレ, ref="なし") 
#これで解釈を間違うことはありません、後の記述は同じです
glm(formula = 効果n ~ 筋トレ, family = binomial, data=df)

(Intercept)   筋トレなし  
     -1.386        1.792 ``
exp(1.792) = 6.001443

筋トレしたらオッズが6倍になる・・・同じ結果なのですが、こちらの方が説得力はありそうです。

それでは、次のような例はどうでしょうか・・・

df2 <- data.frame(
        c(rep("多", 8), rep("中", 5), rep("無", 2),
        rep("多", 1), rep("中", 4), rep("無", 10)),
        c(rep("あり", 15), rep("なし", 15))
        )
colnames(df2) <- c("筋トレ", "効果")
dat2 <- xtabs(~ 筋トレ + 効果, data = df2)  

ロジスティック回帰を行ってみます

あり <- c(8, 5, 2)
なし <- c(1, 4, 10)
筋トレ <- c("多", "中", "無")
glm(cbind(あり, なし) ~ 筋トレ, family = binomial)

(Intercept)     筋トレ中     筋トレ無  
      2.079       -1.856       -3.689  

やはりここでも基準が理解できていないと解釈を間違うことになります．
オッズの基準を「効果なし」、筋トレ の基準を「無」に設定してみます．

筋トレ  <- as.factor(筋トレ ) 
筋トレ  <- relevel(筋トレ , ref=c("無")) 
#効果の基準「なし」を後に書きます
glm(cbind(あり,なし) ~ 筋トレ , family = binomial)

(Intercept)     筋トレ多     筋トレ中  
     -1.609        3.689        1.833  
     
筋トレ「多」のオッズ比
exp(3.689)＝40
筋トレ「無」に比べると「多」のオッズは40倍になることが推定される

効果「あり」を基準にしい場合

glm(cbind(なし, あり) ~ 筋トレ , family = binomial)

(Intercept)     筋トレ多     筋トレ中  
      1.609       -3.689       -1.833  

筋トレ「無」に比べると筋トレ「多」は、効果ありと比較した効果なしの割合は0.025倍に落ちることが推定されます．筋トレすると効果なしの比率が低くなるという結果です・・・研究目的次第だとは思いますが、やはり前の基準が説明しやすいと良いと思います．

基準値の変更

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)

標準誤差を挿入

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)

参考
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歳未満の人口が０％の場合には個人貯蓄合計が平均17.49660になる」と解釈できるのですが、15歳未満の人口が０％というのはあまり現実的な解釈ではありませんので、次のように書き換えます．

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

結果の切片は、15歳未満の人口が1０％の場合の平均的な個人貯蓄合計を意味することになります。切片が変わるだけすので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歳未満の人口が１％増加するごとに個人貯蓄の合計が平均して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

自由度 n=50, p=1 （単回帰なのでp=1となる）
Residual standard error（残差標準誤差）…でも、残差の「標準偏差」を意味している
  
Multiple R-squared （決定係数）
  R²=

決定係数は説明変数の数が増えるほど1に近づくので自由度で調整
Adjusted R-squared（自由度調整済み決定係数）
  AR²=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.4432 ～ 0.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分 運動後3分
1        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)

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

== 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君をお勧めします。