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

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

since2016 ときどきTEXのメモ

t分布からの信頼区間

投稿: 2017-06-12
更新: 2022-1-16

例)p=0.025、自由度5の場合

# 分位点関数 quantile
(q <- qt(0.025, 5, lower.tail = F))
# 確率関数 probability
(p <- pt(q, 5, lower.tail = F))
# 密度関数 density
(d <- dt(q, 5))

サンプルサイズ6の95%信頼区間を求める

#データサンプル
DBP<-c(26,19,28,24,29,32)
vDBP<-var(DBP)

統計量T値の求め方

(t <- mean(DBP) /sqrt(vDBP/6))

95%信頼区間の求め方

mean(DBP)-qt(0.025,5,lower.tail = F)*sqrt(vDBP/6)
mean(DBP)+qt(0.025,5,lower.tail = F)*sqrt(vDBP/6)

Rの関数で確認

t.test(DBP)

One Sample t-test
data:  DBP
t = 14.328, df = 5, p-value = 2.985e-05
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
 21.60893 31.05774
sample estimates:
mean of x 
 26.33333 

t分布のグラフ

x <- seq(-4, 4, 0.1)
plot(x, dt(x, 20), type="l")

for(i in c(2,3,4)){
for(i2 in 2:4){
curve(dt(x, 6-i), col=i, type="l",add=T, )
}
}
labels<-c("自由度","黒:20","赤:4","緑:3","青:2")
legend("topleft", legend = labels)

f:id:yoshida931:20220116223222p:plain:w500

対応のあるデータのグラフ

id <- 1:5
x <- c(rnorm(5, 20, 2))
y <- c(rnorm(5, 25, 1))
data <- data.frame(id, x, y)

  id        x        y
1  1 22.83957 23.13902
2  2 15.05851 23.15993
3  3 19.49754 24.70116
4  4 17.43654 26.60801
5  5 17.01846 25.09940

x <- c(1,2)
t_data <- t(data[,2:3])

      [,1]     [,2]     [,3]     [,4]     [,5]
x 22.83957 15.05851 19.49754 17.43654 17.01846
y 23.13902 23.15993 24.70116 26.60801 25.09940

matplot(x, t_data, type = "l", lty = 1)
legend("topleft", legend=id, col = id,  lty = 1)

f:id:yoshida931:20220125222246p:plain:w500

id <- 1:20
group <- c(rep("A",10),rep("B",10))
x <- c(rnorm(10,20,2), rnorm(10,15,3))
y <- c(rnorm(10,25,1), rnorm(10,20,0.5))
data <- data.frame(id, group, x, y)

x <- c(1,2)
t_data <- t(data[,3:4])
mg <- ifelse(data$group == "A", 1, 2)
matplot(x, t_data, type = "l", lty = 1, pch = 1, col = mg, xlab = "", ylab = "",xaxt="n", xlim = c(0.9,2.1))
name<-c("前","後")
axis(side=1,at=c(1, 2),labels=name)  #指定した場所に挿入
legend(1.8, 15, c("治療A","治療B"), col = c(1, 2), lty = 1)

f:id:yoshida931:20190711140353p:plain:w500

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

 ICC(2, 1) 二要因モデル  

Two-way random effects for Case2

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

下記の文献のICC(2, 1)、ICC(2, k)をRでやってみます
Shrout PE, Fleiss JL. Intraclass correlations: uses in assessing rater reliability. Psychol Bull. 1979;86(2):420-8.

平均平方和
 MSA : mean of square(評価者)
 MSB : mean of square (被験者)
 MSE : mean of square

 N : 全体のサンプル数
 n : 被験者数
 k : 繰り返し数(評価者A、B、C)

使用するパッケージ

library(irr)
library(tidyr)      

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

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

sampleは、下記のページより
統計学入門−第5章統計学入門−第5章
このページでもとめる値

ICC(2,1)0.94295%信頼区間:0.845  0.984
ICC(2,3)0.98095%信頼区間:0.942  0.995

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

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

分散分析のために縦のデータに変換

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

変数を名義変数へ変換

dat2$group <- as.factor(dat2$group)
dat2$ID <- as.factor(dat2$ID)

以下のようなグラフをイメージして進めます

plot(as.integer(dat2$ID), dat2$data, pch=19, col=c("red", "green", "blue"), xlab="被験者要因(i:1~10)", ylab="評価 (A,B,C) ")
abline(h=mean(dat2$data))
text(x=1.8,y=63,"全体平均=59.7")
legend("topleft", legend = c("A", "B", "C"), col = c("red", "green", "blue"), pch =19)

f:id:yoshida931:20210912221407p:plain:w500

要因:評価者被験者
評価者は変量モデル:多くの評価者からたまたま選ばれた3人
複数の評価者 (k = A, B, C ) が複数の被験者 ( n = 1, 2, … , 10 ) に評価したときの評価者間( A, B, C ) の信頼性

分散分析の前に線形回帰について
(要因を被験者としていることに注意)

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

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   16.467      4.625   3.560 0.002235 ** 
groupT2       -1.600      3.270  -0.489 0.630571    
groupT3       -1.800      3.270  -0.550 0.588805    
ID2           12.000      5.971   2.010 0.059684 .  
ID3           22.000      5.971   3.685 0.001696 ** 
ID4           28.000      5.971   4.690 0.000183 ***
ID5           35.000      5.971   5.862 1.50e-05 ***
ID6           54.667      5.971   9.156 3.41e-08 ***
ID7           63.333      5.971  10.607 3.58e-09 ***
ID8           72.000      5.971  12.059 4.66e-10 ***
ID9           79.000      5.971  13.231 1.03e-10 ***
ID10          77.333      5.971  12.952 1.46e-10 ***

Intercept は当てはめ値(予測値)より
fitted(fit2)[1]

分散分析

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

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

group = 評価者 (k:A, B, C )
ID = 被験者 ( n:1, 2, … , 10 )

自由度 平方和 平均平方和 F値
評価者 k - 1 SSA MSA=SSA/(k - 1 ) MSA/ MSE
被験者 n - 1 SSB MSB=SSB/(n - 1) MSB/ MSE
残差 (k - 1)(n - 1) SSE MSE=SSE/(k - 1)(n - 1)
全変動 k*n - 1 SST  
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

算出方法
分散分析モデル
 y _ {ij}=μ+a _ {i}+b _ {j}+(ab) _ {ij} + e _ {ij}  \   \   \   (i=1, 2, 3 \   \   \   j=1,2,\cdots,10)
 a _ {i} ;評価者の効果
 b _ {j} ;被験者の効果
 (ab) _ {ij} ;被験者 jと評価者 iの交互作用

 ICC(2.1) とは、全体の分散に対する b _ {j}の分散の割合
この値が大きくなれば、評価者間の分散が小さい(信頼性がある)と判断する

 α _ {i}の分散=  V(a _ {i})
 b _ {i}の分散= V(b _ {i})
 (ab) _ {ij} e _ {ij} の分散= V(e _ {ij})

 ICC(2.1) =\frac{V(b_ {j})}{V(a _ {i}) + V(b _ {j}) + V(e _ {ij})}

 MSA ≒ k×V(a _ {i}) + MSE
 MSB ≒ n×V(b _ {i}) + MSE

 V(a _ {i}) = \frac{MSA - MSE}{k}
 V(b _ {i}) = \frac{MSB - MSE}{n}

  MSE = V(e _ {ij})

分散分析表より
k = 3
n = 10
MSA = 9.73
MSB = 2462.52
MSE = 53.47

va <- (MSA - MSE) / 10 
vb <- (MSB - MSE) / 3
(ICC_2.1 <- vb / (va + vb +MSE))

95%信頼区間

Fj <- MSA / MSE
c <- (k - 1)*(n - 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*(MSB - FL2*MSE)) / (FL2*(k*MSA+(n*k - n - k)*MSE)+n*MSB))#下限
(ICC_2.1_U <- n*(FU2*MSB - MSE) / ((k*MSA + (n*k - n- k)*MSE) + n*FU2*MSB))#上限

 ICC(2, k) 二要因モデル  

k人の評価者が複数の被験者に対して1回実施した時の平均値の信頼性に関する指標

k人の平均値が対象となるため評価者と誤差の分散が1/k倍となる

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

算出方法

k = 3 #評価者数K1  
n = 10 #被験者数K2   
MSA =  9.73
MSB = 2462.52
MSE = 53.47 

ICC(2,k)の分散

va_k <- (MSA - MSE/3) / (10*3)  #評価者の分散
vb_k <- (MSB - MSE/3) / 3  #被験者の分散
(ICC_2.k <- vb_k / (va_k + vb_k +MSE/3))

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)))   #上限

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

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

下記の文献のICC(1, 1)、ICC(1, k)をRでやってみます
Shrout PE, Fleiss JL. Intraclass correlations: uses in assessing rater reliability. Psychol Bull. 1979;86(2):420-8.

 SST : total sum of squares, sum of squares total
 SSA : sum of squares among groups
 SSE : sum of squares error
 MSA : mean of square
 MSE :mean of square

 N : 全体のサンプル数
 n : 被験者数
 k : 繰り返し数(1回~3回)

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

ICC(1,1)0.94295%信頼区間:0.848  0.984
ICC(1,3)0.98095%信頼区間:0.943  0.995

使用するパッケージ

library(irr)
library(tidyr)

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

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

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

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

分散分析のために縦のデータに変換

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

変数を名義変数へ変換

dat2$group <- as.factor(dat2$group)
dat2$ID <- as.factor(dat2$ID)

 ICC(1, 1) 一要因モデル

stripchart(data ~ group, data = dat2, method="jitter",  jitter = 0.05, vertical=TRUE, pch=20, cex=1.5, xlim=c(0.5,3.5), ylim=c(0, 100), xlab="評価3回実施", ylab="被験者10名")

f:id:yoshida931:20210913230228p:plain:w500

一般的な分散分析の散布図、この図では評価(T1, T2, T3)が要因となります。
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)

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)    15.333      4.046   3.790  0.00115 ** 
factor(ID)2    12.000      5.721   2.097  0.04887 *  
factor(ID)3    22.000      5.721   3.845  0.00101 ** 
factor(ID)4    28.000      5.721   4.894 8.77e-05 ***
factor(ID)5    35.000      5.721   6.117 5.60e-06 ***
factor(ID)6    54.667      5.721   9.555 6.76e-09 ***
factor(ID)7    63.333      5.721  11.070 5.58e-10 ***
factor(ID)8    72.000      5.721  12.585 5.84e-11 ***
factor(ID)9    79.000      5.721  13.808 1.10e-11 ***
factor(ID)10   77.333      5.721  13.517 1.61e-11 ***

Intercept は当てはめ値(予測値)より
fitted(fit1)[1]
ID=1 の平均値と同値

この回帰分析のモデルは以下のようになる
 \boldsymbol{B ^ T}= (b1, b2, b3, \cdots , b10)
 \boldsymbol{X ^ T}=(1, x2, x3, \cdots ,x10)  ダミー変数
 Y =  \boldsymbol{B} ^ T*\boldsymbol{X} + ei
 ID1: x2~x10=0
 ID2: x2=1, x3=\dots=x10=0
 ID3: x2=0, x3=1, x3=\dots=x10=0
以下続く・・・

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

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     

要因が被験者となるので繰り返し数は評価の数
N = 30  #(総サンプル数)    
k = 3   #繰り返し3回(T1, T2, T3)
n = 10  
SSAの自由度
n - 1 = 9    
SSEの自由度
N - n = 20 

★ここでは被験者が要因になっている

自由度 平方和 平均平方和 F値
群間 n - 1 SSA MSA=SSA/(n - 1) MSA/ MSE
郡内 N - n SSE MSE=SSE/(N - n)
自由度 平方和 平均平方和 F値
群間 9 22163 2462.5 50.153
郡内 20 982 49.1

注意)一般的な分散分析と異なり被験者が要因となっている。以下の図をイメージしておく。

plot(as.integer(dat2$ID), dat2$data, pch = 20,  xlab = "被験者要因(n : 1~10)", ylab = "k = 3 (T1, T2 ,T3) ")
# 注意)IDが名義のままだとボックスプロットになります

f:id:yoshida931:20210924214059p:plain:w500

ICC1.1の算出方法

MSA = 22163/9    
MSE = 982/20

分散分析モデル
 y _ {ij}=μ+a _ {i}+ e _ {ij}  \   \   \   ( i=1,2,\cdots,10\   \   \   j=1,2,3 )
 a _ {i}=i番目の被験者の効果

全体の分散 V( y _ {ij})
 a _ {i}の分散 V(a _ {i})
 e _ {ij} の分散を V(e _ {ij})

 V( y _ {ij}) = V(a _ {i}) + V(e _ {ij})

 ICC(1.1) とは、全体の分散 V( y _ {ij}) に対する分散 V(a _ {i}) の割合
 ICC(1.1) =\frac{ V(a _ {i})  }{ V(a _ {i})  + V(e _ {ij}) }

 MSA ≒ k × V(a _ {i}) + V(e _ {ij})
 V(e _ {ij}) = MSE
 V(a _ {i}) =  \frac{MSA - MSE}{k}
 ICC(1.1)  = \frac{V(a _ {i})}{V(a _ {i}) +V(e _ {ij})} = \frac{MSA - MSE}{MSA + (k - 1)MSE}

( ICC1.1 <- ( MSA - MSE ) / (MSA + (k - 1) * MSE) )

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

F1 <- MSA / MSE
FL1 <- F1 / qf(0.975, n - 1, N - n)
FU1 <- F1 / qf(0.025, n - 1, N - n)
#繰り返し数 : 3 ( k = 3 )
( ICC_1.1_L <- (FL1 - 1) / (FL1 + (k - 1)) )#下限
( ICC_1.1_U <- (FU1 - 1) / (FU1 + (k - 1)) )#上限

 ICC(1, k) 一要因モデル   

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

k回の平均値が対象となるため誤差の分散が1/k倍となる

ICC(1, 1) より
 V(a _ {i}) =  \frac{MSA - MSE}{k}
誤差の分散
 V(e _ {ij}) =  \frac{MSE}{k}

 ICC(1.1)  = \frac{V(a _ {i})}{V(a _ {i}) +V(e _ {ij})} = \frac{ \frac{MSA - MSE}{k}} {\frac{MSA - MSE}{k} + \frac{MSE}{k} }
  = \frac{MSA - MSE}{MSA}

# 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.k <- ( MSA - MSE ) / MSA)
# ICC(1,k)の95%信頼区間の求め方 (分散比の信頼区間より)
(ICC_1.k_L <- (FL1 - 1) / FL1)#下限
(ICC_1.k_U <- (FU1 - 1) / FU1)#上限

回帰直線と交互作用

f:id:yoshida931:20210611162531p:plain
パッケージ無しで描く方法

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

f:id:yoshida931:20210529081839p:plain:w500
治療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))#全体の回帰

f:id:yoshida931:20210530104037p:plain

やはり、想定したように治療前の血圧が高い人は治療効果も高くなるようです。この散布図を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()

f:id:yoshida931:20210530103932p:plain
治療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()

f:id:yoshida931:20210530104336p:plain

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

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

単回帰分析の傾きの推定値(最小二乗法)
\frac{\sum{(x_i-\bar{x})(y_i-\bar{y})}}{\sum{(x_i-\bar{x})}^2}

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群の切片の差は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()

f:id:yoshida931:20210530104724p:plain 黒実線:全体の共通回帰
青実践:Aの回帰直線
緑実線:Bの回帰直線
赤実践:全体の共通回帰
赤点線:各群の共通回帰

共通回帰の傾きの求め方

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

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)) 

f:id:yoshida931:20210529224051p:plain

共通回帰の検定

共通回帰直線の傾きが各群のもつ傾きと比較して意味を持つかどうかを検定します。この結果から「治療前BPを共変量として、前後差に影響する因子として考慮しなければならない」という解析の方向性を示すことになります.
f:id:yoshida931:20210529105002p:plain:w300
各群の共通回帰から得られる推定値と各群の平均値との差の平均平方和を残差の平均平方和で除した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:id:yoshida931:20210530090718p:plain

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

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

混合モデル(マルチレベルモデル)の基礎

投稿日:2021.5.2

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

  • 線形混合モデル ( linear\ mixed\ model)
  • 線形混合効果モデル ( linear\ mixed\ effect\ model)
  • 階層線形モデル ( hierarchical\  linear\ model )
  • マルチレベルモデル ( multilevel \ model )・・・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"))

f:id:yoshida931:20210502225944p:plain:w400

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.13827頭増加する傾向にある(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))

f:id:yoshida931:20210502230145p:plain:w400

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

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.1023程度のバラツキがある. 各時期(period)の傾向を考慮したことで、傾きは0.13827から0.10386に減少しているが各グループの牛の数に応じて発症数は増加傾向にあることが理解できる.  時期(period)を考慮した場合、牛の数が1頭増えるごとに、発病する牛は0.10386頭増加する傾向にある.

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

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

# 例)時期1の場合の発症数の推定値
# 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(傾き=0、切片固定のモデル)

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))

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