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

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

since2016 ときどきTEXのメモ

分散分析(繰り返し 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!!!

sample

解析の練習で使うデータサンプル集 随時更新

sample1

症例 <- c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20)
性別 <- c("男性","男性","女性","男性","男性","女性","男性","男性","男性","男性","女性","女性","女性","女性","男性","女性","女性","女性","女性","男性")
治療 <- 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変化量 <- c(12,2,17,3,3,10,18,3,7,10,3,11,13,8,4,10,10,20,6,4)
sample1 <- data.frame(症例,性別,治療,治療前BP,治療後BP,BP変化量 )

sample2
対応のある2群

ID <- c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10)
pre <- c(18, 8, 14, 15, 16, 13, 23, 20, 14, 20)
post <- c(23, 27, 15, 22, 15, 20, 20, 20, 18, 25)
sample2 <- data.frame(ID, pre, post)

sample3
独立した2群

グループ <- c("A","A","B","B","A","B","A","B","B","A","A","A","B","A","B","B","B","B","A","A")
ROMT <- c(18, 8, 23, 27, 14, 15, 15, 22, 15, 16, 13, 23, 20, 20, 20, 20, 18, 25, 14, 20)
sample3 <- data.frame(グループ, ROMT)

sample4
クロス表

sample4
判定 <- c(rep("非該当",30), rep("要支援", 25))
歩行 <- c(rep("可能", 20), rep("不可能", 10), rep("可能", 10), rep("不可能", 15))
sample4 <- data.frame(判定, 歩行)