生物医学統計を学ぶ！: 一般化線形モデル（GLM）でLD50を推定する

LD50

半数致死量（はんすうちしりょう、median lethal dose）とは、急性毒性の指標です。投与した動物の半数が死亡する用量を意味します。
"Lethal Dose, 50%"を略してLD50と書きます。

また、投与した動物の半数が最小限の効果を示す用量である半数効果用量（ED50）も同様の考え方で求めることができます。

一般化線形モデル（Generalized Linear Model, GML）

データが正規分布に従うとき、従属変数 Y と独立変数 X を、パラメーター β でモデル化できます。

\[g(y) = \beta_0x_0 + \beta_1x_1 + \cdots + \beta_nx_n \]

ただし、データが正規分布に従わないときは、各独立変数 X に対応する従属変数 Y の分散が一定でなくなり、データを正規分布と仮定した一般線形モデルが利用できなくなります。

このようなとき、分散が大きくなりすぎないようにする（いい具合に小さく見えるようにする）リンク関数（対数変換など）を導入して、予測性能を安定させたモデルがGMLです。

線形と名前が付いていますが、勾配の最適化も行う点で、最小二乗法と比べて、多少、融通の利くモデルになっています。

実践

サンプルデータとして、投与量と死亡数を記録したデータを利用します。

```
  # データをロード
  # https://drive.google.com/file/d/12_wZfUOmTLpqckCwHF7dk1OniaVjpWzo/view?usp=sharing
  !!curl -L -o "LD50.csv" "https://drive.google.com/uc?export=download&id=12_wZfUOmTLpqckCwHF7dk1OniaVjpWzo"
  # データをテーブルにします
  import pandas as pd
  import numpy as np
  # 各行が群（投与量別）ごとの記録です
  # 1つの群につき10個体ずつの結果です。
  # dose：投与量
  # death：死亡数
  data = pd.read_csv(filepath_or_buffer="LD50.csv",names=["dose","death"] ,encoding="utf-8")
  # statsmodel用に二項分布データを用意する
  alive = (np.zeros(6)+10 - data["death"].values).astype(np.int8)
  data = pd.DataFrame(data={
      "survived":alive,
      "death":data["death"],
      "dose":data["dose"]
      } )
  #生存した割合の変数を作成する
  data['survival-rate'] = data['survived'] / (data['survived'] + data['death'])
  data['mortality'] = 1 - data['survival-rate']
  data['TOTAL'] = data['survived'] + data['death']
  print(data)
  
  # 直線性の適合度検定
  from scipy import stats
  observed_values=data["dose"].values
  expected_values=np.ones(len(observed_values)) * np.mean(data["dose"])#今回は平均値を利用
  # print(observed_values)
  # print(expected_values)
  # p0.05（0.1でもよいくらい）にて直線性に問題ないことを確認
  stats.chisquare(observed_values, f_exp=expected_values)
  
  # 直線性の適合度検定(Chi-Square goodness of fit test)
  # https://www.statisticssolutions.com/chi-square-goodness-of-fit-test/
  # https://www.khanacademy.org/math/ap-statistics/chi-square-tests/chi-square-goodness-fit/v/goodness-of-fit-example
  '''
  死亡数に直線性があるを調べます。
  直線性があれば、容量が増えるごとに死亡数が増えていることが確認できます。
  直線性がない場合は、解析の信憑性が欠けると解釈されます。
  '''
  from scipy import stats
  '''
  確率が0、1となるデータを除外します。
  （直線になりやすいバイアスを除去）
  '''
  observed = data["mortality"].values[1:5] # 0.91
  print("0%, 100%を除外した残り",len(observed))
  print(observed)
  # pvalue>0.05にて直線性に差がないを採択
  # f_expは省略すると平均値で処理される
  # stats.chisquare(f_obs=observed, f_exp=expected, ddof=0)
  stats.chisquare(f_obs=observed, ddof=0)      
  
  # GLMのBinomial（二項分布+logit関数）で解くケース
  '''
  目的変数そのものである確率を予測する
  目的変数の値が仮に100としたら、GLMの予測値が0.5のとき、50%の確率で100であると解釈する。
  半致死線量の場合は、全体の50%の確率値となる説明変数を探せばよい。
  （p%LDは、全体のp％の確率値となる説明変数を探せばよい）

二項分布
  結果が成功か失敗のいずれかである試行（ベルヌーイ試行と呼ばれる）を
  独立に n 回行ったときの成功回数を確率変数とする離散確率分布
  ddof（自由度） = n-1として扱う。

logit関数
  0から1をとる確率に対して、
  逆S字を横向きにしたような値をとる関数
  '''
  # http://www.statsmodels.org/devel/examples/notebooks/generated/glm.html
  import statsmodels.api as sm
  import statsmodels.formula.api as smf

# データフレームからDatasetを作る
  dose_data = np.log10(data["dose"]) # 対数変換（戻すときは10**x）
  exog = sm.add_constant(dose_data) # X 説明変数
  endog = data[["survived","death"]] # Y 目的変数
   # GLMをBinominal + Logitで探索させる
  model = sm.GLM(endog,exog,family=sm.families.family.Binomial(link=sm.genmod.families.links.logit()))# logit is default
  res_binom = model.fit(method="newton") # newton法で勾配最適化
  # モデルの詳細
  print(res_binom.summary()) # constは切片（intercept）
  print('Parameters: ',res_binom.params)
  print('T-values: ', res_binom.tvalues)

# 可視化
  y = endog.values[:,0]/endog.values.sum(axis=1)
  print(y.shape)
  yhat = res_binom.mu
  print("y-hat",yhat)

from statsmodels.graphics.api import abline_plot
  import matplotlib.pyplot as plt
  fig, ax = plt.subplots()
  ax.scatter(yhat, y)
  line_fit = sm.OLS(y, sm.add_constant(yhat, prepend=True)).fit()
  abline_plot(model_results=line_fit, ax=ax)
  ax.set_title('Model Fit Plot')
  ax.set_ylabel('Observed values')
  ax.set_xlabel('Fitted values');

# for new test data
  # 50% 点を探す（もっといい方法ありそう）
  # おおよそ190程度がLD50
  testd = pd.DataFrame(
      {
          "dose":np.log10(np.linspace(np.min(data["dose"]),np.max(data["dose"]),200,endpoint=True))
      }
  )
  # https://stackoverflow.com/questions/51314284/predicting-confidence-interval-with-statsmodels
  pred = res_binom.predict(sm.add_constant(testd["dose"]))
  pred_res = res_binom.get_prediction(sm.add_constant(testd["dose"]))
  conf_intervals = pred_res.conf_int(obs=True,alpha=0.05)# PredictResultインスタンスから両側95%の信頼区間を得る
  # GLMは予測値が目的変数のパーセンテージになる。
  # 今回は、各群で10なので、10=99.99...%
  # for i in range(len(testd["dose"])):
  #   print("dose",10**testd["dose"][i],"pred_y",pred[i],"95ci",conf_intervals[i])

plt.figure()
  plt.plot(10**testd["dose"],pred)
  # print(pred[pred > 0.4999]) # 昇順と降順は要確認
  # 50%の値と位置を取得
  ld50_y = pred[pred > 0.4999][-1:].values
  ld50_pos = np.where(pred.values == ld50_y)
  ld50 = 10**testd["dose"].values[ld50_pos]
  plt.plot([10**testd["dose"].min(),10**testd["dose"].max()],np.repeat(ld50_y,2))
  plt.plot(np.repeat(ld50,2),[0,1])
  plt.title("LD50 = "+str(np.round(ld50,2)))
  plt.show()
                                      
```

Visionary Imaging Services, Inc.

Imaging CROサービスをご提供させていただきます。

お問い合わせお待ちしております。

https://www.vis-ionary.com/