回帰直線と交互作用

パッケージ無しで描く方法

dat <- data.frame(y=c(1,2,3,4), x=c(4,5,8,12), a=c(1,1,0,0))
dat$a <- as.factor(dat$a)

plot(dat$x, dat$y)
fit <- summary(lm(y~x, data=dat))
lines(range(dat$x), fit$coef[1]+ fit$coef[2]*range(dat$x))

pchAB <- ifelse(dat$a == "1", 19, 21)
plot(dat$x, dat$y,  pch=pchAB, cex=1.5 )
legend("bottomright",
    legend=c("a=1", "a=0"),
    pch=c(19, 21),
    lty = c(1, 3)
    )

fit1 <- summary(lm(y ~ x, data=dat[dat$a==1,]))
lines(range(dat$x), fit1$coef[1]+ fit1$coef[2]*range(dat$x))
fit0 <- summary(lm(y ~ x, data=dat[dat$a==0,]))
lines(range(dat$x), fit0$coef[1]+ fit0$coef[2]*range(dat$x), lty=3)

共分散分析 ANCOVA

サンプル（勝手に作ったサンプルです。実存しません。）
治療Aと治療Bの降圧効果に関する検証

治療 <- c("B","A","B","A","A","B","B","A","A","B","A","B","B","B","A","B",
    "B","B","A","A")
治療前BP <- c(160,135,177,141,142,155,175,145,149,155,135,156,170,
    150,138,160,150,180,155,140)
治療後BP <- c(148,133,160,138,139,145,157,142,142,145,132,145,157,
    142,134,150,140,160,149,136)
前後差 <- 治療前BP - 治療後BP
dat1 <- data.frame(治療, 治療前BP, 治療後BP, 前後差)
dat1$治療 <- as.factor(dat1$治療)
head(dat1)

  治療 治療前BP 治療後BP 前後差
1    B      160      148     12
2    A      135      133      2
3    B      177      160     17
4    A      141      138      3
5    A      142      139      3
6    B      155      145     10

治療Bが治療Aよりも降圧効果があるようです！平均値の差の検定をやってみましょう。

t.test(前後差~治療, data=dat1)
t = -6.7187, df = 13.82, p-value = 1.047e-05
95％信頼区間: -11.543307  -5.951643
A群とB群の平均値
3.888889       12.636364 

差がありました。95％信頼区間から6～11程度の差があるようです。しかし、差が大きいのは治療前BPが高い人では・・・という疑問が残ります。

治療前BPと前後差の散布図と回帰直線

fitAll <- lm(前後差 ~  治療前BP,  data = dat1)
anova(fitAll)
fitAllhat <- fitAll$coef[1] + fitAll$coef[2] *dat1$治療前BP
plot(dat1$治療前BP,dat1$前後差, cex=1.5 ,xlab="治療前BP", ylab="前後差") 
lines(range(治療前BP), fitAll$coef[1] + fitAll$coef[2] * range(治療前BP))#全体の回帰

やはり、想定したように治療前の血圧が高い人は治療効果も高くなるようです。この散布図をA群・B群に色分けします。

fig1 <- function() 
    {
    pchAB <- ifelse(dat1$治療 == "A", 19, 21)
    plot(dat1$治療前BP,dat1$前後差, pch=pchAB, cex=1.5 , 
    xlab="治療前BP",  ylab="前後差")
    lines(range(治療前BP),fitAll$coef[1]+fitAll$coef[2]*range(治療前BP),
    col=1)#全体の回帰
    legend("topleft", legend = c("治療A", "治療B"), pch = c(19, 21))
    }
fig1()

治療B群にはもともと血圧の高い人が集まっていたようで、治療効果も大きかったようです。治療前BPが交絡変数として考えられます。この治療前BP因子を考慮したうえでA群とB群の差を考えなければいけません。それぞれの群の回帰直線を挿入してみます。

setA <- subset(dat1, dat1$治療 == "A")
setB <- subset(dat1, dat1$治療 == "B")
fitA <- lm(前後差~治療前BP, data=setA)
fitB <- lm(前後差~治療前BP, data=setB)
fig2 <- function() 
    {
    fig1()
    lines(range(治療前BP),fitA$coef[1]+fitA$coef[2]*range(治療前BP),col="blue")
    lines(range(治療前BP),fitB$coef[1]+fitB$coef[2]*range(治療前BP),col="green")
    }
fig2()

単純に平均値だけを考えても答えは出ないようです。治療前BPが影響していることは明らかなので、治療前BPを共変量として解析します。（今回は、例として治療前BPを共変量としましたが、実際の解析では多くの変数の中から選択することになります）

統計学入門−第8章 より引用
このように他のデータの影響を考慮して目的のデータを総合的に比較する手法を共分散分析(ANCOVA：analysis of covariance、アンコバ)といい、影響を考慮する変数のことを共変数(covariable)といいます。 この手法は比較したいデータを目的変数にし、影響を考慮するデータを説明変数にした回帰分析と、群を要因にした一元配置分散分析を組み合わせた手法に相当します。 共分散とは回帰分析で重要な役目をする値であり、この手法が共分散分析と呼ばれるのは共分散を利用した分散分析だからです。

単回帰分析の傾きの推定値（最小二乗法）

Rを使って解析やってみます。方法は至って簡単で、重回帰分析の説明変数の1つがカテゴリーになっているというだけです。

fit1 <- lm(前後差 ~  治療前BP + 治療, data = dat1)
summary(fit1)

             Estimate Std. Error t value Pr(>|t|)    
(Intercept) -39.79522    4.71524  -8.440 1.75e-07 ***
治療前BP      0.30715    0.03301   9.304 4.41e-08 ***
治療B         2.50511    0.89016   2.814   0.0119 *

共通の傾きは0.30715、２群の切片の差は2.50511。つまり、治療Bの前後差平均値は、治療Bより平均して2.50511高いと推定できます。平均値の差の検定では6～11の差を見積もっていたのですが、誤算であったことが理解できました。

fit1$coef[2]#共通の傾き
fit1$coef[3]+fit1$coef[1]#Bの切片
fit1$coef[1]#Aの切片

fig3 <- function() 
    {
    fig2()
    #共通回帰式
    lines(range(治療前BP), (fit1$coef[3]+fit1$coef[1]*2)/2 + 
    fit1$coef[2]*range(治療前BP), col="red")
    #Aの回帰式（傾きは共通回帰式）
    lines(range(治療前BP), fit1$coef[1]+ fit1$coef[2]*range(治療前BP), 
    lty = 2, col="red")
    #Bの回帰式（傾きは共通回帰式）
    lines(range(治療前BP), fit1$coef[3]+fit1$coef[1]+ 
    fit1$coef[2]*range(治療前BP), lty = 2, col="red")
    }
fig3()

黒実線：全体の共通回帰
青実践：Aの回帰直線
緑実線：Bの回帰直線
赤実践：全体の共通回帰
赤点線：各群の共通回帰

共通回帰の傾きの求め方

まず最初に全体の平均と各群のデータの平均値が重なるようにデータを移動します。
イメージ図

m治療前BP <- mean(dat1$治療前BP)#治療前BPの全体平均
m前後差 <- mean(dat1$前後差)#前後差の全体平均
ma治療前BP <- tapply(dat1$治療前BP, dat1$治療, mean)[1]#Aの治療前BP平均
mb治療前BP <- tapply(dat1$治療前BP, dat1$治療, mean)[2]#Bの治療前BP平均
ma前後差 <- tapply(dat1$前後差, dat1$治療, mean)[1]#Aの前後差平均
mb前後差 <- tapply(dat1$前後差, dat1$治療, mean)[2]#Bの前後差平均
#移動後の治療前BP、前後差
移動後の治療前BP <- ifelse(dat1$治療=="A", dat1$治療前BP + (m治療前BP - ma治療前BP), dat1$治療前BP + (m治療前BP - mb治療前BP))
移動後の前後差 <- ifelse(dat1$治療=="A", dat1$前後差 + (m前後差 - ma前後差), dat1$前後差 + (m前後差 - mb前後差))
#共通の傾き
lm(移動後の前後差 ~ 移動後の治療前BP)$coef[2]

#  傾き＝0.3071539 
# 各群の共通回帰の予測式 y = β0 + 0.3*x にそれぞれの群の平均値と共通回帰の傾きを当てはめて切片β0 を算出する

移動前後の散布図の比較

par(mfrow=c(1,2)) 
plot(治療前BP, 前後差, xlim=c(130,180), ylim=c(0, 20))
lines(range(130:180), fitAll$coef[1]+ fitAll$coef[2]*range(治療前BP))
plot(移動後の治療前BP, 移動後の前後差, xlim=c(130,180), ylim=c(0, 20))
f <- lm(移動後の前後差~移動後の治療前BP)
lines(range(130:180), f$coef[1]+ f$coef[2]*range(130:180))
par(mfrow=c(1,1)) 

共通回帰の検定

共通回帰直線の傾きが各群のもつ傾きと比較して意味を持つかどうかを検定します。この結果から「治療前BPを共変量として、前後差に影響する因子として考慮しなければならない」という解析の方向性を示すことになります．

各群の共通回帰から得られる推定値と各群の平均値との差の平均平方和を残差の平均平方和で除したF値で検定します。共通回帰のF値が大きければ共通回帰が意味を持つことになる。小さい場合には、共通回帰の傾きが0に近いことを意味します。
F値= （AB群の共通回帰の推定値の平均平方和ー交互作用の平均平方和）÷ 残差平方和

setA <- subset(dat1, dat1$治療 == "A")#A群のみのセット
setB <- subset(dat1, dat1$治療 == "B")#B群のみのセット
fitA <- lm(前後差~治療前BP, data=setA)#A群のみの単回帰
fitB <- lm(前後差~治療前BP, data=setB)#B群のみの単回帰
fitAB <- lm(前後差~治療前BP*治療, data=dat1)#全データの回帰、交互作用あり
S1 <- anova(fitA)$Mean[1] + anova(fitA)$Mean[1]#各群の平方和
S2 <- anova(fitAB)$Mean[3]#交互作用の平均平方和
S3 <- anova(fitAB)$Mean[4]#残差の平方和
Fvalue <- (S1 - S2) / S3
#F検定
#帰無仮説：（AB群の共通回帰の推定値の平均平方和ー交互作用の平均平方和）の分散と残差平方和の分散が等しい
pf(Fvalue, 1, 16, lower.tail=F) 
# 0.00198736

非並行性の検定（交互性の検定）

共通回帰のF値が大きく、非平行性のF値が大きい場合には、両群の回帰直線の傾きが非並行ということになり、両群の共通回帰直線が意味を持つことになります。 共通回帰のF値が小さく、非平行性のF値も小さい場合には、共変量の影響を考慮する必要はなく分散分析で解析します。 ​

#交互作用の平均平方和 ÷ 残差平均平方和
f <- S2 / S3
#F検定
#帰無仮説：交互作用の平均平方和の分散と残差平方和の分散が等しい
pf(f, 1, 16, lower.tail=F) 
# 0.05857521

P=0.06ですので、有意水準をどのように設定するかで、A群とB群の非平行性の検定結果は異なります。有意水準は、検定の前に設定しなければなりません。ｐ値から、どのような解析手法にするのか吟味しなければなりません。　　

混合モデル（マルチレベルモデル）の基礎

投稿日:2021.5.2

忘れないように基本的な部分のみ貼っておきますが、名称だけでも統一していただきたいです・・・

• 線形混合モデル ()
• 線形混合効果モデル ()
• 階層線形モデル ()
• マルチレベルモデル ()・・・etc

サンプルはパッケージlme4の「cbpp」を使用します
概要：牛肺疫（CBPP:contagious bovine pleuropneumonia）は、マイコプラズマを原因とするウシおよびスイギュウの感染症。本データセットは、エチオピアのBoji地区にある商業用の牛群（15グループ）における血清学的なCBPP発生状況を示したもの。これら15グループすべての牛から四半期ごとに血液サンプルを採取した追跡調査。 ​

library(lme4)
data(package="lme4") #パッケージのサンプルを見る

head(cbpp)
#    herd incidence size period
#  1    1         2   14      1
#  2    1         3   12      2
#  3    1         4    9      3
#  4    1         0    5      4
#  5    2         3   22      1
#  6    2         1   18      2

herd:牛のグループ名（15グループ）
incidence:各グループの各期間におけるCBPP発症数。
size:各グループに飼育されている牛飼養頭数
period:調査時期

サンプルの目的とは異なりますが「牛飼養頭数とCBPP発症頭数の関連性について」という設定で進めます。

全体の散布図と回帰直線

library(lattice)
xyplot(incidence  ~ size,  data = cbpp, col=1, type = c("p","r"))

fit1 <- lm(incidence  ~ size,  data = cbpp)
summary(fit1)

#             Estimate Std. Error t value Pr(>|t|)   
# (Intercept) -0.31106    0.73554  -0.423  0.67405   
# size         0.13827    0.04389   3.150  0.00266 **

牛の数が1頭増えるごとに、発病する牛は0.14頭増加する傾向にある（p=0.0027）。つまり牛の飼養頭数が10頭増すとCBPP感染牛が1.4頭増える傾向にあることが示唆されている。

調査時期を考慮した場合はどうなるか？

xyplot(incidence  ~ size, data = cbpp, group = period, pch = 19, type = c("p","r"), auto.key = list(pch = 19, corner=c(0,1), border = T, padding.text = 1.5, columns = 1), par.settings = simpleTheme(pch = 16))

調査時期ごとの単回帰分析の結果

md <- lmList(incidence ~ size|period, data = cbpp)
summary(md)  
   
# Coefficients:
#    (Intercept) 
#     Estimate Std. Error    t value  Pr(>|t|)
# 1 -0.5048610   1.731796 -0.2915246 0.7719063
# 2 -0.1375479   1.323172 -0.1039531 0.9176397
# 3  1.3450519   1.390118  0.9675809 0.3381039
# 4 -0.1213868   1.066066 -0.1138642 0.9098204
#    size 
#      Estimate Std. Error    t value    Pr(>|t|)
# 1  0.24666516 0.08851373  2.7867445 0.007603578
# 2  0.08927203 0.07871405  1.1341308 0.262372119
# 3 -0.02452146 0.08995962 -0.2725829 0.786343078
# 4  0.05534211 0.07412426  0.7466127 0.458938830

両変数の関連に有意な傾向を示したのは、時期1のみである。時期3は有意ではないが、やや下降傾向にある。各調査の回帰直線を考慮した上で、モデルを構築する必要があります。

単回帰モデル

fit <- lm(incidence  ~ size,  data = cbpp)
summary(fit)
#            Estimate Std. Error t value Pr(>|t|)   
#(Intercept) -0.31106    0.73554  -0.423  0.67405   
#size         0.13827    0.04389   3.150  0.00266 **

期間をランダム効果として取り入れたモデル
ただし傾きと切片に共分散を仮定（共分散≠0）

fitm1 <-  lmer(incidence  ~ size + (size|period),  data = cbpp)
summary(fitm1)
# Random effects:
#  Groups   Name        Variance Std.Dev. Corr 
#  period   (Intercept) 0.24363  0.4936        
#           size        0.01047  0.1023   -1.00
#  Residual             4.40888  2.0997        
# Number of obs: 56, groups:  period, 4
# Fixed effects:
#             Estimate Std. Error t value
# (Intercept) -0.01094    0.68611  -0.016
# size         0.10386    0.06422   1.617

sizeは全期間を通して0.1程度のバラツキがある. 各時期（period）の傾向を考慮したことで、傾きは0.13827から0.10386に減少しているが各グループの牛の数に応じて発症数は増加傾向にあることが理解できる.  時期（period）を考慮した場合、牛の数が1頭増えるごとに、発病する牛は0.1頭増加する傾向にある.

それぞれの時期における推定値

ranef(fitm1)$period

#   (Intercept)        size
# 1  -0.6635655  0.13756593
# 2   0.1108923 -0.02298944
# 3   0.2900289 -0.06012683
# 4   0.2626443 -0.05444966

# 例）時期１の場合の発症数の推定値
# incidence  = -0.01094 + (-0.01094)*size + (-0.6635655) + 0.13756593*size
#              =(-0.01094 - 0.6635655)+(-0.01094 + 0.13756593)*size

他の混合モデルの例
切片を固定、傾きをランダムにしたモデル

(　 )の中：0＝切片固定、1＝切片ランダム
fitm2 <-  lmer(incidence  ~ size + (0 + size|period),  data = cbpp)
summary(fitm2)

# Random effects:
#  Groups   Name Variance Std.Dev.
#  period   size 0.006005 0.0775  
#  Residual      4.430191 2.1048  
# Number of obs: 56, groups:  period, 4
# 
# Fixed effects:
#             Estimate Std. Error t value
# (Intercept)  0.08343    0.64411   0.130
# size         0.10009    0.05508   1.817

傾きを固定、切片をランダムにしたモデル

fitm3 <-  lmer(incidence  ~ size + (1|period),  data = cbpp)
summary(fitm3)

# Random effects:
#  Groups   Name        Variance Std.Dev.
#  period   (Intercept) 1.506    1.227   
#  Residual             4.800    2.191   
# Number of obs: 56, groups:  period, 4
# 
# Fixed effects:
#             Estimate Std. Error t value
# (Intercept)  0.24068    0.91873   0.262
# size         0.09878    0.04136   2.389

傾き、切片をランダムにしたモデル
最初に作ったfitm1と同じモデル

fitm4 <-  lmer(incidence  ~ size + (1 + size|period),  data = cbpp)
summary(fitm4)

#  Groups   Name        Variance Std.Dev. Corr 
#  period   (Intercept) 0.24363  0.4936        
#           size        0.01047  0.1023   -1.00
#  Residual             4.40888  2.0997        
# Number of obs: 56, groups:  period, 4
# 
# Fixed effects:
#             Estimate Std. Error t value
# (Intercept) -0.01094    0.68611  -0.016
# size         0.10386    0.06422   1.617

null model（period切片をランダムにした切片固定モデル）

fitm5 <-  lmer(incidence  ~ 1 + (1|period),  data = cbpp)
summary(fitm5)

# Random effects:
#  Groups   Name        Variance Std.Dev.
#  period   (Intercept) 2.219    1.490   
#  Residual             5.124    2.264   
# Number of obs: 56, groups:  period, 4
# 
# Fixed effects:
#             Estimate Std. Error t value
# (Intercept)    1.714      0.804   2.131

null modelと比較

lm(incidence  ~ 1 , data = cbpp)

傾きと切片をそれぞれ独立したランダム変数としたモデル

lmer(incidence ~ size + (1|period) + (0 + size|period), data = cbpp)

グループを考慮した場合はどうなるか？

各グループ最大4回調査を実施していますので、各時期の頭数と発症数の関連性を見てみます。

xyplot(incidence  ~ size, data = cbpp, group = herd, pch = 19, type = c("p","r"), auto.key = list(pch = 19, corner=c(0,1), border = T, padding.text = 1.5, columns = 1), par.settings = simpleTheme(pch = 16))

15グループあるので非常に複雑で、しかも切片・傾きともにランダム（バラバラ）であることには間違いないようです・・・といった具合にマルチレベルの分析は果てしなく続きます。

級内相関係数 (ICC：Intraclass Correlation Coefficient)

投稿日:2020.7.30
最終更新日: 2021.4.18

下記の文献をRを使って解析します
Shrout PE, Fleiss JL. Intraclass correlations: uses in assessing rater reliability. Psychol Bull. 1979;86(2):420-8.

: between-targets mean square
: within-targets mean square
: between-judges mean square
: residual (error) mean square

sampleは、下記のページより
統計学入門−第5章統計学入門−第5章

ICC(1,1)≒0.942、95％信頼区間：0.848  0.984
ICC(1,3)≒0.980、95％信頼区間：0.943  0.995

ICC(2,1)≒0.942、95％信頼区間：0.845  0.984
ICC(2,3)≒0.980、95％信頼区間：0.942  0.995

ICC(3,1)≒0.938、95％信頼区間：0.831  0.983
ICC(3,3)≒0.978、95％信頼区間：0.936  0.994

iccを実行するパッケージ

library(irr)

下記のような引数が用意されています

icc(ratings, model = c("oneway", "twoway"), 
    type = c("consistency", "agreement"), 
    unit = c("single", "average"), r0 = 0, conf.level = 0.95)

データセットをRにコピーします

ID <- 1:10
A <- c(15,30,34,52,58,69,76,88,91,95)
B <- c(10,14,42,38,51,78,88,90,94,87)
C <- c(21,38,36,40,42,63,72,84,98,96)
dat1 <- data.frame(ID, A, B, C)

分散分析のために縦のデータに変換
パッケージtidyr

library(tidyr)
dat2 <- dat1 %>% gather(group, data, -ID) 

一要因モデル

1人の評価者が被験者 ( n = 10 ) に対して複数回 ( k = 3回 ) 評価を実施した時の信頼性

icc(dat1[,-1], model ="oneway", type ="consistency", unit = "single")  
# ICC(1) = 0.942
# 95%-Confidence Interval for ICC Population Values:
# 0.848 < ICC < 0.984

分散分析の前に線形回帰について

fit1 <- lm(data ~ factor(ID), data=dat2)
summary(fit1)

この回帰分析のモデルは以下のようになる






以下続く・・・

ここでは分散分析表を使います

anova(fit1)
#            Df Sum Sq Mean Sq F value    Pr(>F)    
# factor(ID)  9  22163  2462.5  50.153 9.014e-12 ***
# Residuals  20    982    49.1     

上記の分散分析のResidualsがWithin targetsとなります

	自由度	平方和	平均平方和	F値
Between targets	n - 1	BSS	BMS	BMS / EMS
Within targets	n(k - 1)	WSS	WMS	WMS / EMS

	自由度	平方和	平均平方和	F値
Between targets	9	22163	2462.5	50.153
Within targets	20	982	49.1

算出方法

n = 10 #(被験者数)      
k = 3 #(評価回数 or 評価者数)          
BMS=2462.5    
WMS=49.1 

分散分析モデル

番目の被験者の効果
とは、全体の分散に対するの分散の割合
の分散を、 の分散を とした場合、

とは分散分析よりすでに算出済み

；k回(3回)評価しているのでｋをかける

( ICC1.1 <- (BMS - WMS) / (BMS + (k-1) * WMS) )

ICC(1,1)の95%信頼区間の求め方 (分散比の信頼区間より)

F1 <- BMS / WMS
FL1 <- F1 / qf(0.975, n-1, n*(k-1))
FU1 <- F1 / qf(0.025, n-1, n*(k-1))
( ICC_1.1_L <- (FL1 - 1) / (FL1 + (k-1)) )#下限
( ICC_1.1_U <- (FU1 - 1) / (FU1 + (k-1)) )#上限

一要因モデル 　　

One-way random effects for Case1
1人の評価者が被験者 ( n = 10 ) に対して複数回 ( k = 3回 ) 評価を実施した時の評価平均値の信頼性に関する指標で、 の分散をkで割った値を使用する
は、 に対するの分散

# unit = averageを使用する
icc(dat1[,-1], model ="oneway", type ="consistency", unit = "average")
# ICC(3) = 0.98     
# 95%-Confidence Interval for ICC Population Values:       
# 0.943 < ICC < 0.995      

算出方法
ICC(1.1)と同様に


よりを求める

# A, B, Cの平均になるので kで除した値を使用する
(ICC_1.k <- ( BMS - WMS ) / BMS)
# ICC(1,k)の95%信頼区間の求め方 (分散比の信頼区間より)
(ICC_1.k_L <- (FL1 - 1) / FL1)#下限
(ICC_1.k_U <- (FU1 - 1) / FU1)#上限

二要因モデル 　

Two-way random effects for Case2

評価者のA, B, Cは、たまたま選ばれた3名（変量モデル）
同じ評価を実施したときに、いつも同じ評価者ではないことが前提となっている。
評価を実施するたびに評価者が異なるので、評価者を変数扱いとなる。
複数の評価者 ( k=3; A, B, C ) が複数の被験者 ( n = 10 ) に評価したときの評価者間の信頼性

fit2 <- lm(data ~ group + factor(ID), data = dat2)    
anova(fit2)

#           Df  Sum Sq Mean Sq F value   Pr(>F)    
#group       2    19.5    9.73   0.182   0.8351    
#factor(ID)  9 22162.7 2462.52  46.051 1.33e-10 ***
#Residuals  18   962.5   53.47                     

icc(dat1[,-1], model ="twoway", type ="agreement", unit = "single")
# ICC(A,1) = 0.942
# 95%-Confidence Interval for ICC Population Values:
# 0.844 < ICC < 0.984

算出方法
分散分析モデル

；評価者の効果   randam variable
；被験者の効果
；被験者と評価者の交互作用

とは、全体の分散に対するの分散の割合
の分散=
の分散=
の分散=
の分散=

上記の分散分析のResidualsの平均平方和が となります

分散分析表より
JMS = 9.73
BMS = 2462.52
EMS = 53.47
# ICC_2.1 = ((BMS - EMS)/k) / ((BMS - EMS)/k + (JMS - EMS)/n + EMS)  より
(ICC_2.1 <- ( BMS - EMS ) / ( BMS + (k-1)*EMS + k*(JMS - EMS)/n) )

95%信頼区間

Fj <- JMS / EMS
c <- (n-1)*(k-1)*(k*ICC_2.1*Fj + n*(1+(k-1)*ICC_2.1) - k*ICC_2.1)^2
d <- (n-1)*k^2*ICC_2.1^2*Fj^2+(n*(1+(k-1)*ICC_2.1)-k*ICC_2.1)^2
(FL2 <- qf(0.975, n-1, round(c/d,0)))
(FU2 <- qf(0.975, round(c/d,0), n-1))
(ICC_2.1_L <- (n*(BMS - FL2*EMS)) / (FL2*(k*JMS+(n*k-n-k)*EMS)+n*BMS))#下限
(ICC_2.1_U <- n*(FU2*BMS - EMS) / ((k*JMS + (n*k - k- n)*EMS) + n*FU2*BMS))#上限

二要因モデル 　

複数の評価者 ( k=3; A, B, C ) が複数の被験者 ( n = 10 ) に評価したときの平均値の信頼性

icc(dat1[,-1], model ="twoway", type ="agreement", unit = "average")
# ICC(A,3) = 0.98
#  95%-Confidence Interval for ICC Population Values:
#  0.942 < ICC < 0.995

算出方法
は、 に対するの割合

# ICC_2.k = ((BMS - EMS)/k) / ((BMS - EMS)/k + ((JMS - EMS)/n + EMS)/k)  より
(ICC_2.k <- ( BMS - EMS ) / ( BMS + (JMS - EMS)/n ))

95%信頼区間

(ICC_2.k_L <- (k*ICC_2.1_L / (1+(k-1)*ICC_2.1_L)))#下限
(ICC_2.k_U <- (k*ICC_2.1_U / (1+(k-1)*ICC_2.1_U)))#上限

二要因モデル 　

Two-way mixed model for Case3

特定の評価者の信頼性を検討したいときに使用する。同じ試験を何度も実施したときに、評価者は常に同じであるため定数扱いとなる。被験者については変量モデルなので、混合モデルと呼ばれる場合もある。

icc(dat1[,-1],  model ="twoway", , type ="consistency", unit = "single")

# ICC(C,1) = 0.938
# 95%-Confidence Interval for ICC Population Values:
# 0.831 < ICC < 0.983

算出方法
分散分析モデルはICC2.1と同じだが、評価者の効果は定数扱いとなる

；評価者の効果   fixed effect
；被験者の効果
；被験者と評価者の交互作用
の分散＝0

全体の分散  

評価者の効果は定数扱いとなるので、
ICC(3, 1)は、 から を引いた値に対するの割合

anova(fit2)
BMS <- 2462.52
EMS <- 53.47
# ICC_3.1 <- ((BMS - EMS)/k) / (((BMS - EMS)/k) + EMS) より
(ICC_3.1 <- (BMS - EMS) / (BMS + (k - 1)*EMS))

95%信頼区間

FL3 <- (BMS/EMS) / ( qf(0.975, n-1, (n-1)*(k-1)) )
FU3 <- (BMS/EMS) * ( qf(0.975, (n-1)*(k-1), n-1) )
(ICC_3.1_L <- (FL3 - 1) / (FL3 + (k - 1)))#下限
(ICC_3.1_U <- (FU3 - 1) / (FU3 + (k - 1)))#上限

二要因モデル 　

クロンバックのα係数、エーベルの級内相関係数 r11
「特定の評価者(k=3人)」が1回評価したときの「評価平均値」の信頼性

icc(dat1[,-1],  model ="twoway", , type ="consistency", unit = "average")

# ICC(C,3) = 0.978
# 95%-Confidence Interval for ICC Population Values:
# 0.936 < ICC < 0.994

算出方法

；評価者の効果   fixed effect
；被験者の効果
；被験者と評価者の交互作用
の分散＝0

全体の分散（ 評価平均値なので、残差の効果はをで除した値となる）

(ICC_3.k  <- (BMS - EMS) /BMS )

95%信頼区間

(ICC_3.k_L <- 1 - (1/FL3))#下限
(ICC_3.k_U <- 1 - (1/FU3))#上限

ロジスティック回帰の基準(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}