Home

多変量重回帰モデルのコントラスト検定

 

多変量重回帰モデルのコントラスト検定をgeneral fundamental least-squares theorem（Timm, 1975, p. 191）に基づいて行う方法（岡本、2014、第3.4節）のためのスクリプトを、Pythonで作成した。スクリプトをリスト１に示す。スクリプトファイルなどは、圧縮ファイルmvmregfiles.zipにまとめた。ダウンロード・展開して自由に利用されたい。

リスト１のスクリプトをコマンド

 

PS ・・・\mvmregfiles> python MVmreg.py

 

で実行すると、入力データファイル名を聞いてくるので、次のように入力する。

 

入力データファイル名(*.txt) = Data.txt

 

上のデータファイルData.txtは、拙著「心理学データ分析と測定」の表3.2.1のデータのファイルである。リスト２にファイルData.txtのリストを示す。

 

表3.2.1

A校の生徒の点数

B校の生徒の点数

(41, 77)     (35, 77)

(74, 29)     (51, 67)

(58, 55)     (50, 69)

(61, 49)     (67, 53)

(60, 52)     (51, 67)

(44, 61)     (53, 64)

(62, 54)     (84, 30)

(61, 45)     (77, 41)

(67, 46)     (68, 41)

(53, 60)     (39, 76)

(55, 64)     (54, 77)

(43, 87)     (83, 47)

(61, 63)     (84, 42)

(62, 62)     (86, 27)

(60, 63)     (35, 89)

(73, 50)     (64, 48)

(79, 47)     (48, 79)

(58, 74)     (54, 66)

(68, 54)     (90, 29)

(74, 43)     (49, 69)

 

入力データファイルData.txtは、以下の形式で用意されている。

 

/       Y1  Y2  X1  X2

        2       2

        41        77         1         0

        35        77         1         0

        74        29         1         0

　　　　　・

　　　　　・

　　　　　・

        90        29         0         1

        74        43         0         1

        49        69         0         1

/     データの終わり

 

 

/　　コントラストC

        1

        1               -1

       

/　　A

        1

        1

        0

       

/　　コントラストC

        1

        1               -1

       

/　　A

        1

        0

        1

       

/　　コントラストC

        1

        1               -1

       

/　　A

        1

        1

        1

       

/　　　終わり

        -1

       

 

データファイルは、スラッシュ / で始まる行によって区切られたテキストファイルである。スラッシュ / の右側は、無視されるのでコメントを書くことができるが、日本語を用いるときには注意が必要である。入力データファイルは、次のコードで読み込んでいる。

 

with open(infilenm, 'r', encoding = 'utf-8') as fin:  #   日本語はUTF-8で設定

    raw_data = fin.readlines()

 

 

データ作成時に、日本語コードが「encoding = 'utf-8'」の設定に合わないときは、読み込み時にエラーとなる。使用した日本語コードにencodingの設定を合わせるか、日本語を用いないでデータファイルを作成する必要がある。

最初のスラッシュ / で始まる行の次の行に、従属変数の数pと独立変数の数qを空白文字で区切って並べる。

次の行から、１行に１ケースのデータを、従属変数、独立変数と並べる。i番目のケースのデータが、従属変数Y(i,1)、・・・、Y(i,p)、独立変数X(i,1)、・・・、X(i,q)であれば、1行に

Y(i,1)    Y(i,2)  ・・・  Y(i,p)   X(i,1)  X(i,2)  ・・・  X(i,q)

と書く。すべてのデータを書き終わったら、スラッシュ / で始まる行を置いてデータの終わりであることを示す。

次のスラッシュ / で始まる行の次の行には、コントラスト行列Cの行数を書き、続く行から１行に１つずつコントラストを書いていく。

Cの値を書き終わったら、また、スラッシュ / で始まる行を置いて、その次の行にはコントラスト行列Aの列数を書き、続く行からAの値を行単位で書いていく。行列Aの値を読み込むと、そのCとAの組み合わせに対して検定が行われる。

このCとAの組を書くことを必要な回数繰り返す。CとAによる制約条件については、拙著「心理学データ分析と測定」における次式

の説明を参照されたい。プログラムでは

の場合を扱っている。

 

Cの行数として負の値を書くと、入力の終了となり、スクリプトの実行終了である。

端末画面には、次のメッセージが表示されている。

 

Results.txt was saved.

 

出力ファイルResults.txtを開くと以下のようになっている。

 

入力データファイル = Data.txt

 

 

        Y1        Y2        X1        X2

      41.0      77.0       1.0       0.0

      35.0      77.0       1.0       0.0

      74.0      29.0       1.0       0.0

　　　　　・

　　　　　・

　　　　　・

      90.0      29.0       0.0       1.0

      74.0      43.0       0.0       1.0

      49.0      69.0       0.0       1.0

 

 

 

C =

[[ 1. -1.]]

 

A =

[[1.]

 [0.]]

 

Lambda = 0.95204,   df1 = 1,  df2 = 1,  df3 = 38

 

F = 1.91414,   df1 = 1,  df2 = 38,  p = 0.17458

 

Chi^2 = 1.84292,   df = 1,   p = 0.17461

 

 

C =

[[ 1. -1.]]

 

A =

[[0.]

 [1.]]

 

Lambda = 0.98866,   df1 = 1,  df2 = 1,  df3 = 38

 

F = 0.43593,   df1 = 1,  df2 = 38,  p = 0.51307

 

Chi^2 = 0.42775,   df = 1,   p = 0.51310

 

 

C =

[[ 1. -1.]]

 

A =

[[1.]

 [1.]]

 

Lambda = 0.53633,   df1 = 1,  df2 = 1,  df3 = 38

 

F = 32.85175,   df1 = 1,  df2 = 38,  p = 0.00000

 

Chi^2 = 23.36263,   df = 1,   p = 0.00000

 

 

入力データの出力後、各コントラスト行列CとAの組に対する検定統計量が出力されている。

Lambda = に続く値は、に従う統計量の値である。

F の値は、コントラスト行列Aの列数u が１のときに出力される。

Chi-square の値は、バートレット近似（Bartlett’s approximation）の値である。

 

上の出力結果の意味は、拙著「心理学データ分析と測定」の「3.2　分析例２」を参照されたい。

 

 

 

引用文献

岡本安晴（2014）心理学データ分析と測定．勁草書房

Tim, N. H. (1975). Multivariate analysis with applications in education and psychology. Brooks/Cole publishing company.

 

 

リスト１　多変量重回帰モデルのコントラスト検定（MVmreg.py）

 

"""

        N. H. Timm (1975).

          Multivariate analysis with applications in

                education and psychology.

             Brooks/Cole publishing company

 

        岡本安晴（2014）心理学データ分析と測定．

        　　第３章　偏回帰係数を比べる：コントラスト

        　勁草書房　　　　　

"""

import numpy as np

import scipy.stats as ss

import sys

 

infilenm = input('入力データファイル名(*.txt) = ') 

with open(infilenm, 'r', encoding = 'utf-8') as fin:  #   日本語はUTF-8で設定

    raw_data = fin.readlines()

 

fout = open('Results.txt', 'w')

fout.write(f'入力データファイル = {infilenm}\n')

 

pre_data = ''

for v in raw_data:

    pre_data += v.replace('\t', ' ').replace('\n', ',')

pre1_data = pre_data.split(',')

ck = 0

pos = 0

while True:

    if len(pre1_data[pos]) > 0:

        if pre1_data[pos][0] == '/':

            pos += 1

            break

    pos += 1

p_q = np.array(pre1_data[pos].split(), dtype='int')

p = p_q[0]                 #   従属変数の個数

q = p_q[1]                 #   独立変数の個数

print('p =',p, '  q =',q)

#

#           データY,Xの読み込み

#

pos += 1

preX = []

while True:

    if pre1_data[pos][0] == '/':

        pos += 1

        break

    v = np.array(pre1_data[pos].split(), dtype='float')

    print(v)

    preX.append(v)

    pos += 1

preX = np.array(preX)

#print('preX =\n', preX)

 

Y = preX[:,0:p]

X = preX[:,p:p+q]

r = np.linalg.matrix_rank(X)

N = len(X)

if p > N-r:

    print(f'\nthat p <= N-r should hold.\np = {p},  N = {N},  r = {r}\n')

    sys.exit()

 

fout.write('\n\n')

for i in range(p):

    vs = f'Y{i+1}'

    fout.write(f'{vs:>10s}')

for i in range(q):

    vs = f'X{i+1}'

    fout.write(f'{vs:>10s}')

fout.write('\n')

for i in range(N):

    for j in range(p):

        vs = f'{Y[i][j]}'

        fout.write(f'{vs:>10s}')

    for j in range(q):

        vs = f'{X[i][j]}'

        fout.write(f'{vs:>10s}')

    fout.write('\n')

   

ginvXX = np.linalg.pinv(X.T @ X)

H = ginvXX @ X.T @ X

Bhat = ginvXX @ X.T @ Y

print('Bhat =\n', Bhat)

#

#            コントラストの検定

#

while True:

    while True:

        if len(pre1_data[pos]) > 0:

               if pre1_data[pos][0] == '/':

                   pos += 1

                   break

        pos += 1

    g = int(pre1_data[pos])

    if g < 1:                 #   検定の終了

        break

    if g > r:

        print(f'g(={g}) should be smaller or equal to f(={r})!\n')

        continue

    pos += 1

    #

    #               Cの読み込み

    #

    C = []

    for i in range(g):

        C.append(pre1_data[pos].split())

        pos += 1

    C = np.array(C, dtype='float')

    print('\n\nC =\n', C)

    fout.write('\n\n')

    fout.write('\nC =\n')

    fout.write(C.__str__())

    if np.trace((C - C@H) @ (C - C@H).T) > 1.0e-9:

        print('C is not estimable\nC =\n', C)

        continue

    while True:

        if len(pre1_data[pos]) > 0:

               if pre1_data[pos][0] == '/':

                   pos += 1

                   break

        pos += 1

       

    u = int(pre1_data[pos])

    pos += 1

    if u > p:

        print(f'\nu(={u}) should be smaller or eual to p(={p})\n')

    #

    #        Aの読み込み

    #

    A = []

    for i in range(p):

        A.append(pre1_data[pos].split())

        pos += 1

    A = np.array(A, dtype='float')

    print('A =\n', A)

    fout.write('\n\nA =\n')

    fout.write(A.__str__())

    #

    #      岡本安晴（2014）、pp.45-47.

    #

    preQea = (Y - X @ Bhat) @ A

    Qea = preQea.T @ preQea

    Qta = Qea + (C @ Bhat @ A).T @ np.linalg.inv(C @ ginvXX @ C.T) @\

          (C @ Bhat @ A)

    Lmbd = np.linalg.det(Qea) / np.linalg.det(Qta)

    print('\nLambda =', Lmbd, ' df1 =', u, ' df2 =', g, ' df3 =', N-r)

    fout.write('\n\n')

    fout.write(f'Lambda = {Lmbd:.5f},   df1 = {u},  df2 = {g},  df3 = {N-r}')

    if u == 1:

        F = ((1.0 - Lmbd)/g) / ((Lmbd/(N-r)))

        print(f'\nF = {F:.5f},   df1 = {g}, df2 = {N-r},', end ='')

        cdf = ss.f.cdf(F, g, N-r)

        print(f'     p = {1-cdf:.5f}')

        fout.write('\n\n')

        fout.write(f'F = {F:.5f},   df1 = {g},  df2 = {N-r},')

        fout.write(f'  p = {1-cdf:.5f}')

 

    chi2 = -((N-r) -0.5*(u-g+1)) * np.log(Lmbd)

    print(f'\nChi~2 = {chi2:.5f},   df = {u*g},', end = '')

    cdf = ss.chi2.cdf(chi2, u*g)

    print(f'    p = {1-cdf:.5f}')

    fout.write('\n\n')

    fout.write(f'Chi^2 = {chi2:.5f},   df = {u*g},   p = {1-cdf:.5f}')

 

fout.close()

print('\nResults.txt was saved.')

 

 

リスト２．データファイル（Data.txt）

 

/       Y1  Y2  X1  X2

        2       2

        41        77         1         0

        35        77         1         0

        74        29         1         0

        51        67         1         0

        58        55         1         0

        50        69         1         0

        61        49         1         0

        67        53         1         0

        60        52         1         0

        51        67         1         0

        44        61         1         0

        53        64         1         0

        62        54         1         0

        84        30         1         0

        61        45         1         0

        77        41         1         0

        67        46         1         0

        68        41         1         0

        53        60         1         0

        39        76         1         0

        55        64         0         1

        54        77         0         1

        43        87         0         1

        83        47         0         1

        61        63         0         1

        84        42         0         1

        62        62         0         1

        86        27         0         1

        60        63         0         1

        35        89         0         1

        73        50         0         1

        64        48         0         1

        79        47         0         1

        48        79         0         1

        58        74         0         1

        54        66         0         1

        68        54         0         1

        90        29         0         1

        74        43         0         1

        49        69         0         1

/     データの終わり

 

 

/　　コントラストC

        1

        1               -1

       

/　　A

        1

        1

        0

       

/　　コントラストC

        1

        1               -1

       

/　　A

        1

        0

        1

       

/　　コントラストC

        1

        1               -1

       

/　　A

        1

        1

        1

       

/　　　終わり

        -1

       

 

 

 

Home