尤度比検定
対応のあるデータの差の検定
モデル
がモデル
のパラメータに制約を加えたものであるとき、制約モデルが正しいかどうかの検定に尤度比が用いられる。それぞれのモデルにおけるパラメータの最尤推定値に対する尤度関数値（最大尤度）を
と
で表すとき、尤度比
を考える。制約モデルが正しければ、尤度比
はあまり小さくならないと予想できる。
制約モデルが正しいとき、
の分布は、データ数が多くなると自由度
のカイ二乗分布に近付く。ここで、
と
は、それぞれモデルとのパラメータの張る空間の次元（個数）である（岡本、2009、p．153参照）。
フルモデルと比べて制約モデルのデータに対する当てはまりが悪いときはの値は小さくなり、
の値は大きくなる。
データ数が多くなると、帰無仮説「制約モデルが正しい」の下で対数尤度比の−2倍の分布が自由度
のカイ二乗分布に近付くという性質を用いて制約モデルの検定を行うスクリプトを、本ウェブサイトの後半に示すリスト１のように用意してみた。スクリプトファイルとサンプルデータファイルは、まとめて圧縮ファイルfilesMLLRatio.zipとした。ダウンロード解凍すれば、自由に利用できる。
リスト１では、最尤推定値を、最小値を求める関数scipy.optimize.minimizeを以下のように
Results =
scipy.optimize.minimize(self.OptFunc, x)
用いて求めるため、尤度関数の符号を反転した関数OptFuncの最小値を求めるという形にしている。関数OptFuncは、例えば、以下のように与えられている。
def
OptFunc(self, x): # x[0] = mu, x[1]**2 = sd > 0.0
return -ss.norm.logpdf(self.difs, loc = x[0], scale = x[1]**2).sum()
リスト１のスクリプトは、図１のようなデータを分析するものである。
図１ Excelデータファイル（DiffData.xlsx）
図１のデータは、２０人の人の健康法実行前と実行後の血圧を測ったものである（仮想データ）。データファイルの第１行目には、変数名が書かれている。データは２行目以降に書く。１列目は、データの識別値で図１の場合は通し番号が用いられているが、文字列でもよい。また、同じ識別値を繰り返し用いてもよい。２列目と３列目にデータ値が書かれている。
第２列目のデータを
、第３列目のデータを
としたとき、差
の分布が正規分布に従う
とするのがフルモデルであり、平均０の正規分布に従う
とするのが制約モデルである。
したがって、
である。
リスト１では、フルモデルがクラス型FullModelで、制約モデルがクラス型NoDifModelで表されている。
リスト１を実行すると、まず入力データファイル名の設定が求められる（図２）。
図２
入力データファイルは、図１の形式で用意されたものを指定する。図２では、図１のファイルDiffData.xlsxが設定されている。
入力データファイル名を設定すると、次に計算結果の出力ファイル名の設定が求められる（図２）。ファイル拡張子が「.txt」である任意のテキストファイル名を設定する。
出力テキストファイル名を設定すると、次はグラフの画像保存ファイルの出力ファイル名の設定が求められる（図２）。
画像保存ファイル名を設定すると、データが読み込まれて計算が始まる。計算が終了すると、図３のグラフが表示される。
図３
図３のウィンドウを閉じると、スクリプトの実行終了となる（図４）。
図４
保存されたファイル名が表示されている。これらは、スクリプトの実行終了後、開くことができる。
出力テキストファイルをPython IDLEのエディタで開くとき、図５のメッセージボックスが表示されることがある。
図５
「s-jis」と設定して「OK」ボタンをクリックすると、開くことができる。
図６
図１のデータに対する出力ファイルを開くと、以下のような内容である。
入力データファイル: DiffData.xlsx
前 後
133.00 123.00
133.00 127.00
136.00 129.00
136.00 132.00
134.00 127.00
136.00 127.00
131.00 118.00
139.00 134.00
135.00 129.00
140.00 130.00
132.00 124.00
130.00 122.00
138.00 127.00
138.00 124.00
130.00 118.00
137.00 119.00
131.00 123.00
135.00 128.00
136.00 120.00
134.00 129.00
前：
平均値＝134.70 標準偏差＝2.88
後：
平均値＝125.50 標準偏差＝4.48
差（前 - 後）：
[10, 6, 7, 4, 7, 9, 13, 5, 6, 10, 8,
8, 11, 14, 12, 18, 8, 7, 16, 5]
平均値＝9.20 標準偏差＝3.71
フルモデル：
mu = 9.20 sd = 3.71
制約モデル：
mu = 0.0 sd = 9.92
―２対数尤度比 = 39.34550
p値 = 0.00000 df = 1
ｐ値が0.00000と非常に小さいが、これは図３のグラフにおいて制約モデルのグラフとフルモデルのグラフが著しく異なっていることとよく対応している。
参考文献
岡本安晴（2009）データ分析のための統計学入門．おうふう．https://ci.nii.ac.jp/ncid/BA9003407X
import numpy as np
import scipy.stats as ss
import scipy.optimize as so
import pandas as pd
import matplotlib.pyplot as plt
finnm = input('入力データファイル（*.xlsx）= ')
txtoutnm = input('出力ファイル（*.txt）= ')
foutnm = input('画像保存ファイル（*.png）= ')
txtout = open(txtoutnm, 'w')
#
# Excelファイルの読み込み
# モジュール xlrd がインストール済みであること
#
xlsx_fl = pd.ExcelFile(finnm)
# ExcelファイルからExcelFile型オブジェクトを生成
data_xlsx =
pd.read_excel(xlsx_fl, header =
None) # ExcelFile型オブジェクトからDataFrame型オブジェクトを生成
data = data_xlsx.values
# Dataframe型オブジェクトからnumpyの配列オブジェクトを生成
print('\n', data)
txtout.write('\n入力データファイル: {}\n'.format(finnm))
VarNameX = data[0, 1]
VarNameY = data[0, 2]
print(VarNameX, VarNameY)
txtout.write('\n{0:>10s}{1:>10s}\n'.format(VarNameX,
VarNameY))
X = data[1:, 1]
Y = data[1:, 2]
for v in zip(X, Y):
txtout.write('{0:>10.2f}{1:>10.2f}\n'.format(v[0], v[1]))
meanX = np.mean(X)
meanY = np.mean(Y)
sdX = np.std(X)
sdY = np.std(Y)
print(meanX)
print(meanY)
print(sdX)
print(sdY)
txtout.write('\n\n{}：\n'.format(VarNameX))
txtout.write(' 平均値＝{0:.2f} 標準偏差＝{1:.2f}\n'.format(meanX, sdX))
txtout.write('\n{}：\n'.format(VarNameY))
txtout.write(' 平均値＝{0:.2f} 標準偏差＝{1:.2f}\n'.format(meanY, sdY))
DiffXY = list(X - Y)
# 差をリスト型に変換
print(DiffXY)
meanXY = np.mean(DiffXY)
sdXY = np.std(DiffXY)
print('\n\n差（{0} - {1}）：\n'.format(VarNameX, VarNameY))
print(DiffXY)
print(' 平均値＝{0:.2f} 標準偏差＝{1:.2f}\n'.format(meanXY, sdXY))
txtout.write('\n\n差（{0} - {1}）：\n'.format(VarNameX, VarNameY))
txtout.write('{}\n'.format(DiffXY))
txtout.write('\n 平均値＝{0:.2f} 標準偏差＝{1:.2f}\n'.format(meanXY, sdXY))
class FullModel:
"""
平均値と標準偏差をパラメータとするモデル
"""
def __init__(self, difs):
self.difs = difs
def OptFunc(self, x): # x[0] = mu, x[1]**2 = sd > 0.0
"""
―対数尤度： 最小化を図る目的関数
"""
return -ss.norm.logpdf(self.difs, loc = x[0], scale = x[1]**2).sum()
def EstParams(self):
x =
[0.0, 1.0]
Results = so.minimize(self.OptFunc, x)
return
Results.x[0], Results.x[1]**2
class NoDifModel:
"""
平均値＝０ とするモデル
"""
def __init__(self, difs):
self.difs = difs
def OptFunc(self, x): # x[0]**2 = sd > 0.0
"""
―対数尤度： 最小化を図る目的関数
"""
return -ss.norm.logpdf(self.difs, loc = 0.0, scale = x[0]**2).sum()
def EstParams(self):
x =
[1.0]
Results = so.minimize(self.OptFunc, x)
return Results.x[0]**2
flmdl = FullModel(DiffXY) # フルモデル用オブジェクト
muf, sdf = flmdl.EstParams()
print('\nフルモデル：\nmu = {0:.2f} sd = {1:.2f}'.format(muf, sdf))
txtout.write('\n\nフルモデル：\nmu = {0:.2f} sd = {1:.2f}\n'.format(muf, sdf))
nodifmdl = NoDifModel(DiffXY) # 制約モデル用オブジェクト
sd0 = nodifmdl.EstParams()
print('\n制約モデル：\nmu = 0.0 sd = {0:.2f}'.format(sd0))
txtout.write('\n\n制約モデル：\nmu = 0.0 sd = {0:.2f}\n'.format(sd0))
#
# ―２対数尤度比
#
Neg2LogLmbd = -2 *
(ss.norm.logpdf(DiffXY, loc = 0.0, scale = sd0).sum() - \
ss.norm.logpdf(DiffXY, loc = muf, scale = sdf).sum())
print('\n―２対数尤度比 =', Neg2LogLmbd)
txtout.write('\n―２対数尤度比 =
{0:.5f}\n'.format(Neg2LogLmbd))
df = 2 - 1
# 自由度
p_v = 1.0 - ss.chi2.cdf(Neg2LogLmbd,
df) # ｐ値
print('p値 = {0:.5f} df = {1}'.format(p_v, df))
txtout.write('\np値 = {0:.5f} df = {1}'.format(p_v, df))
# 日本語表記の設定
plt.rcParams['font.family'] = ['Yu
Gothic', 'Noto Sans CJK JP', 'Hiragino Maru Gothic Pro']
x_values = np.arange(min(DiffXY),
max(DiffXY), (max(DiffXY) - min(DiffXY)) / 100)
m1_values = ss.norm.pdf(x_values, muf,
sdf)
m0_values = ss.norm.pdf(x_values, 0.0,
sd0)
plt.title('-2log$\Lambda$ =
{0:.2f}, df =
{1}'.format(Neg2LogLmbd, df), fontsize = 20)
plt.hist(DiffXY, density = True)
plt.plot(x_values, m1_values, lw = 3,
label = 'フルモデル')
plt.plot(x_values, m0_values, lw = 3,
label = '制約モデル')
plt.legend(fontsize = 12)
plt.yticks([])
plt.xlabel('{0} -
{1}'.format(VarNameX, VarNameY), fontsize = 16)
plt.savefig(foutnm)
plt.show()
txtout.close()
print('\n保存ファイル：')
print(' ', txtoutnm)
print(' ', foutnm)