生物医学統計を学ぶ！: 共分散分析

共分散分析の概要

分散分析では交絡要因を共変量として解析する共分散分析（ ANalysis of COVAriance, ANCOVA）という手法があります。
共分散分析は、言い換えれば、「相関の検定」と「対応のある1要因の多群の分散分析」とを組み合わせた手法です。
独立変数による従属変数値の変動の中から、従属変数と交絡のある共変量（剰余変数ともよばれます。背後に存在する原因と結果への影響因子、あるいは、結果と“共”に変わる原因と思われる“変量”です）の影響を極力排除して群を比較することを目的として用いられます。
例えば、血圧を下げる薬AとBを使って、薬剤Aの投与群と薬剤Bの投与群をつくり、年齢と、投与後の血圧を集めたとします。
年齢が高い人は血圧が高い傾向（相関）があったとします。
この場合、独立変数は薬の種類、従属変数は血圧、共変量（剰余変数）は年齢です。
このデータを使って共分散分析を行う場合、その帰無仮説は「年齢の影響を考慮しても群間の差はない」です。
N元配置分散分析との違いは、相関を加味するか否かです。

データの構造

ANCOVAでは共変量という役割をデータに設定します。
データは次のようになります。
※このデータは群を列方向に記載していますが、実際に解析に用いるデータは、行方向にまとめます。

薬A群			薬B群
ID	年齢（共変量）	血圧	ID	年齢（共変量）	血圧
1	45	190	101	60	210
2	58	200	202	57	180

フロー

共分散分析は基本的にパラメトリックな検定です（ノンパラメトリックな検定として扱う場合は、従属変数を順位尺度に変換してからANCOVAを実行するなどの手法なども検討されています）。

まず最初に、相関分析を行い、共変量として設定される変数とその他の変数とで相関があるかを確かめておきます。

他の分散分析と同様に、変数の正規性を各群で確認しておきます。

相関がない場合は共分散分析は行えません。

相関がある場合、線形回帰を確認します。

回帰直線が原点を通る場合は、共変量を除いた変数での分散分析に切り替えます。

回帰直線が原点を通る場合の分散分析は、従属変数を共変量で除した値を用いると、共変量の影響を無くすことができます。

回帰直線が原点を通らないことが確認できたら、回帰直線の平行性を確認します。

共変量ともう一方の変数とは相関が認められているため、直線回帰を仮定した上で、群ごとに回帰直線を計算します。

群間の回帰直線同士の勾配（傾き）が大きく異なる場合、平行性は否定され、群間比較は信用できないと考え、共分散分析は適切でないと判断します。

適切と判断された場合は、ようやくANCOVAを実行できます。

実践

2つの薬剤を対象に年齢と血圧の調査を実施したサンプルデータを用いて、2つの薬剤群の血圧に年齢（共変量）によって差があるかどうかを検討したいとします。

```
  # pingouinパッケージを利用
  !pip install --upgrade pingouin
  # データをロード
  # https://drive.google.com/file/d/1_Hezs310TWhuDDUvxKgDVKBkz__gjsLp/view?usp=sharing
  !curl -L -o "ANCOVA.csv" "https://drive.google.com/uc?export=download&id=1_Hezs310TWhuDDUvxKgDVKBkz__gjsLp"
  # データをテーブルにします
  import pandas as pd
  data = pd.read_csv(filepath_or_buffer="ANCOVA.csv",names=["Group","AGE","BP"] ,encoding="utf-8")
  print(data.head())
  # 正規性を確認します。
  from scipy import stats
  # データを群ごとに検定します。
  for name in ["A1","A2"]:
    one_group = data[data["Group"] == name]
    stat,p = stats.shapiro(one_group["BP"])
    print("norm dist test",name,stat,p) # 正規性が仮定できることが分かります
  # 次に等分散性を調べます。
  stat_l,p_l = stats.levene(data[data["Group"] == "A1"]["BP"],data[data["Group"] == "A2"]["BP"], center="mean")
  print("levene test",stat_l,p_l) # 等分散が仮定できることが分かります
  # 次に相関を調べます
  # 上記から正規性が確認できたので、ピアソンを使います。
  from scipy.stats import pearsonr
  # correlation coefficient, and 2-tailed p-value
  cc,p = pearsonr(data["AGE"],data["BP"])
  print(cc,p) # 相関がないという帰無仮説に対し、p0.05で相関があると仮定できます
  # 次に平行性(Parallelism)をt-testで確認します
  # F-testが良く用いられますが、t-testでもf-testと同じ傾向を示します。
  # Ref：Pharmaceutical Statistics: Practical and Clinical Applications, Fifth Edition
  # https://books.google.co.jp/books?id=vwHMBQAAQBAJ&pg=PA460&lpg=PA460&dq=parallelism+f+test+scipy+stat&source=bl&ots=eaZ_lalylU&sig=ACfU3U349FTAjZcqWVhf9FeVE2NaDOIEyg&hl=ja&sa=X&ved=2ahUKEwjopMOl-frnAhWBH3AKHedPD54Q6AEwAHoECAQQAQ#v=onepage&q=parallelism%20f%20test%20scipy%20stat&f=false
  # https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.stats.linregress.html
  from scipy import stats
  from scipy.stats import t
  import numpy as np
  # 自由度df degree of freedom
  df1 = len(data[data["Group"] == "A1"]["AGE"].values) -1
  df2 = len(data[data["Group"] == "A2"]["AGE"].values) -1
  NumOfSlope = 2
  df = df1 + df2 - NumOfSlope
  print("degree of freedom",df)
  # 共変量はX、従属変数はY
  slope1, intercept1, r_value1, p_value1, std_err1 = stats.linregress(data[data["Group"] == "A1"]["AGE"],data[data["Group"] == "A1"]["BP"])
  slope2, intercept2, r_value2, p_value2, std_err2 = stats.linregress(data[data["Group"] == "A2"]["AGE"],data[data["Group"] == "A2"]["BP"])
  # 傾きが同じといえる可能性
  t_val4S = (slope1 - slope2)/np.sqrt(std_err1**2 + std_err2**2)
  p_s = t.cdf(t_val4S,df,loc=0,scale=1)
  print("Slope:p",p_s)# 今回はp>=0.05で傾きが同じと判断できます

# この例では必要ないが、切片の検定も同様に行える。
  # （この例では原点を通るか否かが重要であるため、切片の有意差は不要）
  # t_val4I = (intercept1 - intercept2)/np.sqrt(std_err1**2 + std_err2**2)
  # p_i = t.cdf(t_val4I,df,loc=0,scale=1)
  # print("Intercept:p",p_i)

# 原点を通るかどうか確認します
  print("切片が0かどうか：",intercept1,intercept2)#37.72559906487433 37.0107338700268、原点は通らない。
  
  #ANCOVAを実行します
  from pingouin import ancova
  ancova(data=data, dv='BP', covar='AGE', between='Group') # AGEは有意に差がある0.000030
  ```

Visionary Imaging Services, Inc.