科研費シンポジウム「多変量データ解析法における理論と応用」
科学研究費・基盤研究（A) 「大規模複雑データの理論と方法論の総合的研究」
研究 代表者：青嶋 誠（筑波大学）, 課題番号：15H01678)■ 日時：2018年12月13日～12月15日■ 場所：広島大学理学部B707■ 開催責任者：栁原宏和 (広島大学)
12月13日13:00-13:10 Opening

セッション1, 座長：若木宏文 (広島大学)
13:10-13:50 矢田和善* (筑波大学)

高次元共分散構造に関する検定の一般化について

青嶋誠（筑波大学）
石井晶（東京理科大学）
13:50-14:30 新村秀一 (成蹊大学)

癌の高次元遺伝子解析の諸問題（2）

14:30-15:10 柳本武美 (統計数理研究所)

RCT と DNN が医療水準の向上を駆動する

企画セッション1「疫学における統計的挑戦」，座長：伊藤ゆり (大阪医科大学)(オーガナイザー：伊藤ゆり・福井敬祐 (大阪医科大学))
15:30-16:10 鈴木有佳 (大阪医科大学)

女性の経済・雇用に関する市区町村の指標と死亡率との関連

16:10-16:50 福井敬祐 (大阪医科大学)

マイクロシミュレーションによる大腸がん検診における死亡率減少効果

16:50-17:30 堀芽久美
(国立がん研究センター)

シミュレーションモデルによる乳がん予防・検診の効果

12月14日
セッション2, 座長：橋本真太郎 (広島大学)
9:30-10:10*

清水康希 (東京理科大学)
橋口博樹（東京理大学）

特異ウィシャート行列の最大・最小固有値分布について岩下登志也（東京理大学）
10:10-10:50 米永航志朗* (北海道大学)
鈴川晶夫 (北海道大学)
10:50-11:30 牧草夏実 (千葉大学)

ウィッシャート行列と正規ベクトルの積の密度関数とモーメント
ヒルベルト空間における正規性の検定の実際的側面

11:30-12:10 岩下登志也* (東京理科大学)

楕円対称性の必要性検定

Bernhard Klar
(カールスルーエ工科大学)
橋口博樹（東京理科大学）
企画セッション2 「社会疫学における統計解析の役割」，座長：福井敬祐 (大阪医科大学)(オーガナイザー：伊藤ゆり・福井敬祐 (大阪医科大学))
13:30-14:10 那波伸敏 (東京医科歯科大学) 地域のソーシャルキャピタル，親の心理的苦痛，
児童虐待の関連：足立区子どもの健康・生活実態
調査から14:10-14:50 伊角彩 (東京医科歯科大学)

子どもの貧困と虐待における親の心理的ストレスと個人レベルのソーシャル・キャピタルの媒介効
果特別セッション, 座長：山村麻理子(広島大学)(平成29年度ダイバーシティ研究環境実現イニシアティブ（牽引型）「国際型ダイバーシティ研究環境実現プログラム」協力)
15:10-15:50 伊藤ゆり(大阪医科大学)

がん登録資料を用いたがん患者の生存時間解析

15:50-16:30 ソルヴァン加藤 比呂子*

複雑な海洋生態システムにおける因果性推測(Inst. of Marine Research)
Sam Subbey(Inst. of Marine Research)
16:30-17:10 イリチュ美佳 (筑波大学)

データ融合に基づく要約手法

12月15日
セッション3, 座長：栁原宏和 (広島大学)
9:30-10:10

永井勇 (中京大学)

Plug-in optimization method for generalized
ridge regression for MLE in GMANOVA model

10:10-10:50 内藤貫太 (千葉大学)

歪曲度のノンパラメトリック推定

10:50-11:30 地道正行* (関西学院大学)

探索的財務ビッグデータ解析

阪智香 (関西学院大学)
宮本大輔 (東京大学)
永田修一 (関西学院大学)
11:30-12:10 大瀧慈* (広島大学 &広島大学ベンチャー企業(株)deCult)
大谷敬子 (広島大学 &広島大学ベンチャー企業(株)deCult)
12:10-12:20 まとめ

セミパラメトリック混合効果モデルの適用による繰り返し観測データの解析高次元共分散構造に関する検定の一般化について
矢田 和善
青嶋 誠
石井 晶

(筑波大数理物質)(筑波大数理物質)(東京理大理工)

1. はじめに共分散行列に p 次の半正定値行列 Σ をもつ母集団を考える．母集団から n (≥ 4)
個の p 次データベクトル x1 , ..., xn を無作為に抽出する．ただし，xj = (x1j , ..., xpj )T
とおく．ここで，σ = tr(Σ)/p，1p = (1, ..., 1)T ，Σ の j 番目の対角成分を σj とし，ΣD = diag(σ1 , ..., σp ) ，ΣIC = σ{(1 − ρ)I p + ρ1p 1Tp } とおく．ただし，ρ ∈ (0, 1) である．
そのとき，Schott (2005, Biometrika) や Srivastava et al. (2011, JMA) は，正規分布もしくはそれに類する仮定のもとで，以下の仮説における高次元検定方式を与えた．

H0a : Σ = ΣDvs. H1a : Σ ̸= ΣD(1.1)

一方で，Srivastava and Reid (2012, JMA) は，正規分布のもと以下の仮説における高次元検定方式を与えた．

H0b : Σ = ΣICvs. H1b : Σ ̸= ΣIC(1.2)

本講演では，まず，仮説 (1.1) と (1.2) を考えた．Yata and Aoshima (2013, JMA) で与えた「拡張クロスデータ行列法 (ECDM)」を用いることで，非正規分布においても頑健な高次元検定方式を与えた．さらに，ΣD や ΣIC を包含するような一般的な高次元共分
散構造の枠組みにおける検定方式も提供した．

2. 仮説 (1.1) の検定方式
いま，∆D = ||Σ − ΣD ||2F = tr(Σ2 ) − tr(Σ2D ) とおく．H0a のもと ∆D = 0 となることに注意する．Yata and Aoshima (2013) は，拡張クロスデータ行列法 (ECDM) という方法論を開発した．2 つの集合 Vn(1)(k) , Vn(2)(k) (k = 3, ..., 2n − 1) を次のように定義する．

Vn(1)(k) =Vn(2)(k) ={{⌊ k2 ⌋ − n(1) + 1, ..., ⌊ k2 ⌋},⌊ k2 ⌋ ≥ n(1) のとき,

{1, ..., ⌊ k2 ⌋} ∪ {⌊ k2 ⌋ + n(2) + 1, ..., n}, それ以外{
⌊ k2 ⌋ ≤ n(1) のとき,{⌊ k2 ⌋ + 1, ..., ⌊ k2 ⌋ + n(2) },{1, ..., ⌊ k2 ⌋ − n(1) } ∪ {⌊ k2 ⌋ + 1, ..., n},

それ以外

ここで，n(1) = ⌈n/2⌉, n(2) = n − n(1) , ⌊x⌋ は x 以下の最大の整数，⌈x⌉ は x 以上の最小の整数を表す．そのとき，拡張クロスデータ行列法を使えば，tr(Σ2 ) の不偏推定量は次のように構築できる．各 k (= 3, ..., 2n − 1) で 2 分割した集合について，標本平均をx̄n(1)(k) =∑1n(1) j∈Vn(1)(k)xj ,x̄n(2)(k) =∑1n(2) j∈Vn(2)(k)

xj

とし，
n
2cn ∑
{(xi − x̄n(1)(i+j) )T (xj − x̄n(2)(i+j) )}2
Wn =n(n − 1) i<jを定義する．ここで，cn = (n(1) − 1)−1 (n(2) − 1)−1 n(1) n(2) である．そのとき，E(Wn ) =tr(Σ2 ) となる．tr(Σ2D ) についても同様に，各 s (= 1, ..., p) と k (= 3, ..., 2n − 1) で∑∑−1xs(1)(k) = n−1
x
,
x
=
n
xsj
sj
s(2)(k)(1)(2)
j∈Vn(1)(k)j∈Vn(2)(k)

とし，
p
n
}2
2cn ∑ ∑ {(xsi − xs(1)(i+j) )(xsj − xs(2)(i+j) ) .
Un =n(n − 1) i<j s=1を定義する．そのとき，E(Un ) = tr(Σ2D ) となる．ここで，∆D の不偏推定量

b D = Wn − Un∆を考える．そのとき，次の定理が成り立つ．
定理 1.

(C-i)

m = min{p, n} とし，次を仮定する．
tr(Σ2 )
→ 0 as m → ∞.
n∆D

そのとき，いくつかの正則条件のもと次が成り立つ．

bD∆= 1 + oP (1) as m → ∞.
∆D
定理 2.

次を仮定する．
{ n∆ }
D(C-ii) lim sup< ∞.
tr(Σ2 )m→∞
そのとき，いくつかの正則条件のもと次が成り立つ．

b D − ∆D∆
⇒ N (0, 1) as m → ∞.
2tr(Σ2 )/n
ここで，H0a のもと tr(Σ2 ) = tr(Σ2D ) となることに注意し，仮説 (1.1) を次のように検
定する．
bD
n∆
rejecting H0a ⇐⇒> zα .
(2.1)
2Un
ただし，zα は N (0, 1) の上側 100α% を表す．そのとき，検定方式 (2.1) による第 1 種の過誤は α に収束する．
当日は，ECDM を用いて仮説 (1.2) の高次元検定方式も与えた．さらに，ΣD や ΣICを包含するような一般的な高次元共分散構造の枠組みにおける検定方式も提供した．

癌の高次元遺伝子解析の諸問題(2)
成蹊大学

名誉教授

新村秀一
現在、再び統計学が注目されている。そのキーワードの一つがビッグデータであり、量、多様性、速度が説明に用いられている。筆者は、既存の小標本を扱う推測統計学との関係で論じたい。小標本はｎと p が小さいものを指すが境界は個人にとって異なる。n が多くて p が小さい縦長データは、統計分析として従来の延長線上で扱える。n が小さく p が大きいものを高次元データあるいは横長データと呼ぶことにし、本研究のテーマである。n と p がともに大きいものがビッグデータであり、ユーザーが勝手に登録したデータを利用できる Apple の個人の顔認識や、Google はペイジランクを数理計画法(MP)のマルコフ過程の応用で特許を取り、大規模なデータをインターネットでコストをかけずに集めた統計手法の集計が主役である。横長データの代表は Microarray である。ハーバード大学の Golub 教授らが 1999 年にサイエンスに発表した論文
で 30 年以上 Microarray の遺伝子情報から 2 群の予測(判別)と新しい癌の Subclass の発見を目指してきたが芳しい成果を得ていない。すなわち、1970 年頃から行われてきたことになる。これらの研究で用いられた Microarray が無償で公開され、統計や機械学習の研究者が 2000 年以前に横長データの新しい研究テー
マとして研究をしてきた。2 群の判別分析が最も適した手法であるが、めぼしい成果が出ていない。
筆者は 1971 年に大学卒業後に大阪成人病センターで、心電図の自動診断論理を判別手法で開発を行ったが医師の開発した枝分かれ論理に全くかなわなかった。これが判別分析の研究の動機になり、多くの実証研
究を通して判別分析の 4 つの問題を見つけた。2015 年の科研費シンポジュームで石井さんの発表で上記の6 種のデータが公開されていることを知り、1970 年以降研究され成果の出ていない判別分析の問題 5 を再認識し、新しく開発した判別理論の応用問題とした。そして僅か 56 日間で簡単に解決し癌の遺伝子解析を完成させた。結論は最小誤分類数(MNM)基準による改定 IP-OLDF(RIP)で 6 種の Microarray の 2 群は完全に分かれ線形分離可能(LSD)である(Fact3)。さらにそれが 64 組から 239 組の遺伝子数が ni 個以下の部分空間(SM)と雑音空間に分割できた（Fact4）。本研究では、なぜ可能であるかを報告する。
1

癌の遺伝子解析6 種類の Microarray データを RIP で判別すると簡単に LSD すなわち MNM=0 である(Fact3)。そして、n 個以下の判別係数だけが 0 でなく、残りが全て 0 になる。さらに判別すると幾つかの係数が 0 になった。これを繰り返し判別し 0 になる判別係数が現れないモデルを SM1 と呼ぶことにした。そして SM1 を省いて判別すると別の SM2 が得られた。そこで LINGO を用いた新手法 2 を RIP で作成した。そして、多くの SM の信号部分空間と雑音空間（MNM>=1）に簡単に分離できた。
2癌の遺伝子診断
当初、癌の遺伝子診断までが統計の役割と考えていた。しかし SM は小標本であるので簡単に統計分析ができる。そこで 8 種類の判別分析と、標準的な統計手法で分析した。しかし RIP、改定 LP-OLDF、改定 IPLPOLDF、H-SVM とロジスティック回帰だけが全ての SM を MNM あるいは NM が 0 と正しく判別できた。それ以外の判別分析は全ての SM を正しく判別できなかった。Alon らで 56 個の SM を見つけ、RatioSV(SV 間の距離の2 を判別スコアの範囲で割った比(％))で評価したところ、RIP と H-SVM の範囲は[3%,25.3%]と[5.8%,32.6%]である。このような簡単な判別は研究対象にならないが、1970 年から行われ 2004 年以降の早い時期に NIH
が終息宣言を出し、医学研究が行われていない状況のようだ。6 種類の Microarray データで得られた全ての SM で求めた RIP の判別スコア(RipDSs)を遺伝子の代わりに変数とした新しい 56 次元のデータ(Signal
1

Data)を作成した。これを階層型クラスター分析で分析すると、2 群が完全に 2 つのクラスターに分かれた。
そして健常群と癌群を PCA で分析すると、図 1 のように Prin1 上の負の軸上に class1 の症例が布置し、正の軸上に class2 の症例が布置することが分かった。これは RIP に限らず改定 LP-OLDF あるいは H-SVM の判
別スコアで作った Signal Data でも、6 種類の Microarray でも、ほぼ同じであるので、2 群が LSD であることが信号の定義でないかと考えられる。しかし他の多くの研究では信号と雑音の定義が明確でない。
3

問題 5 が困難な理由
多くの研究で、癌の遺伝子解析が困難な３つの理由を取り上げているが適切でない。研究がひと段落し、あまりにも簡単に解決したので、恐らく多くの研究者が疑っているのではないかと考えた。そこで多次元遺伝子データの諸問題の事実や問題を以下のように整理することにした。
Fact3：Microarray データを正しく LSD と判別できるのは、RIP、改定 LP-OLDF と H-SVM だけである。統計的判別関数が対応できないのは LSD を正しく判別できないためである。分散共分散行列に基づき相関比
最大化基準では LSD 判別に正しく対応できない。青嶋と矢田は、多くの統計研究者と異なり高次元空間のMicroarray を PCA で分析すると第 1 固有値が極端に大きくなることから、2 群が高次元空間で Prin1 上にある 2 つの球上に布置することを示した。彼らの結果との整合性を探してきた。56 個の SM の和集合を RIPで判別すると、SV=-1 と SV=1 になる症例が多くの割合を示し、SV が-1 以下あるいは SV が 1 以上になる症
例が少ない。それを 56 個の SM に分割すると図 2 のよう SV 上の割合が少なく裾が長くなるようだ。
Fact4：さらに、RIP と改定 LP-OLDF は Microarray を多数の SM(MNM=0)とそれ以外の雑音空間(MNM>=1)に分割できる。しかし H-SVM は Golub らの Microarray で 2000 個程度の係数を 0 にするが、SM に分割できない。この理由は意外と簡単なことが分かった。MP の中で 2 次計画法(QP)だけが、目的関数を最大/最小にする最適解を定義域で 1 個見つける。
これは統計的判別関数でも同じである。このため、部分空間で Full Model
より良いモデル(部分空間すなわち癌の遺伝子)を見つけるため“Feature selection”が必要になる。P が 1
万個もあれば、NP-hard になる。しかし、線形計画法(LP)で定式化された改定 LP-OLDF は、簡単に n 個以下の最適な MNM=0 になる SM を見つけることができる。図 1 は n=3 で p=2(p<n)の一般的な事例である。そして
図中の三角形がこの判別モデルの実行可能領域になる。即ち三角形の全ての点が目的関数を最適化する点になる。この実行可能領域は通常の数理計画法の常識と異なり、どの点も全て最適解になる。LP は、3 つの三角形の頂点を選ぶ。この頂点は、一般的に n=3 の中の 2 個の症例の超平面の交点である。一方横長データ(n<<p)では、
最大 n 個の交点から選ばれる。
少し無理はあるが、3 番目のケースを省くと横長データの n=2<p=3の状況を作れる。この場合、ケース 1 とケース 2 と b1 軸で囲まれた三角形が実行可能領域になる。ケース
1 とケース 2 の交点は p=2 の点であるが、ケース 1 とケース 2 が b1 軸と作る交点は、b2=0 であり 1 次元の部分空間になる。1 万次元空間の実行可能領域のイメージは難しいが、各凸体の頂点は(p-n)個以上の遺伝子の判別係数を 0 にしたもので、容易に(p-n)個以上の判別係数が 0 になる。
以上の説明を聞いてもしっくり納得できない点が残る。なぜこれほど統計的判別関数や統計手法が役に立たないかである。報告書の図 3 は SM を主成分分析で分析し、Prin1 を横軸に、Prin2 と Prin3 の散布図である。この SM の遺伝子を用いた RIP の判別スコアは各軸に直交しないで、かつ長さも短い。主成分軸に投影すればさらに短くなる。即ち完全に大きなバラツキの中に埋もれている。癌と健常は明確に分離しているが、そのばらつきは全体の中で埋もれている。この 2 群が分離しているという信号は、RIP,改定 LPOLDF、H-SVM とロジスティック回帰でしかわからないためである。当初、SM で 2 群は完全に分かれているので、これに含まれる遺伝子を癌遺伝子と考えた。しかし、統計的判別関数や統計手法は LSD である兆候を示さない。適切な判別関数で作られる総合特性軸で初めて 2 群が分かれていることが分かる。しかし、2群のばらつきが全体に比べ小さく直交もしていないので、信号を検出することが困難であるが結論である。
2

 と  が医療水準の向上を駆動する
柳 本 武 美 統計数理研究所

 序       無作為化比較試験 と  は   
標題の内  は    深層神経ネット を表す。前者は       比較臨床試験 の略として使われることも多い。本稿では  を分類器に限って議論する。

 定着した  と途上の 
今日の日本では  は定着している。当初は患者に最善の治療を施す必要があるのに無作為に割り付けた実験 試験 を実施することに反感をもつ医療関係者が多くて定着しなかった。しかし、結核・リュウマチなどの分野から始まり広く受け入れられている。臨床研究の主流を占めるに至っている。一端  が導入されると真面目な医療関係者が多い我が国の文化の中で、質の高い試験が実施されている。この事情は以前の標本調査の導入と軌を一にしている。
近年の目覚ましい発展に伴い  への期待は大きい  。一方で懐疑的な意見をもつ研究者も多い。一時的なブームではないかとする研究者も多い  。人の感覚・知覚・認識を模して学習をさせようとする試みは、今日で言う機械学習の当初からあった。基本的な枠組みは今日でも変わっていない。特筆すべきブレークスルーがあったと言うより、計算機器と計算技能の向上とデー
タの保存・転送技術の進歩・普及によっている。基本的枠組みが変わっていないことは、その推測形式がごく常識的であることを示唆している。

 適用範囲の違いと広さ
は医療行為の効果を確認するために実施される。医療行為は一面では意図的な侵襲行為であるから、誤った印象を受けやすい。実際歴史的には無効な医療行為が歓迎されていた例が多い。一方、 では機械的処理で解決できそうなタスクが多いことから歓迎されている。そのタスクには、囲碁・将棋などゲームの戦術なお、多様な分野が含まれることが既に実証されている。


医療行為の効果を評価するタスクを  で実施することは困難である。医師の経験と印象による評価が危ういことは既に経験済みである。勿論  による効果の評価更には新しい治療法を探索する試みは一部で標榜されているが、著しい成功例は見当たらない。このことは  と の守備範囲が異なることを示している。治療の効果は医療行為の中核であるから、 の役割は強調するまでもない。治療法と効果は区々であって統一的な科学法則は見つけがたい。多くの各個の試験が必要とされる。
一方で、 が解決に大きく寄与すると期待されるタスクは多い。検査画像の分類を通じた診断支援はその実用化が視野に入っている。医療データの一つの特性は、検査時あるいは診察時の情報が、以後により精度の高い情報が得られる可能性が大きいことである。その結果多層的な推論の可能性が期待できる。
心理学では人が物事を認識・理解する段階を三つに分ける。  感覚     知覚     認識      である。従来の画像認識の用語は高い認識機能を目標にしたことを示している。 は感覚的データから知覚的反応を導くタスクに対して性能の良い結果を示した。他の方法で情報科学の分野で探ると、医療分野での成功例が乏しい。

 推論形式の類似性
と  の推論形式は見た目にも大きく異なる。しかし丁寧に見ると共通点が多い。これ

らの共通点は、常識的な科学観に基づく推論の限界を示唆する。

  統計的方法としての は将に統計的方法の典型である。単に解析の段階でのみ主要な役割を果たすのではない。
データ収集の骨格を解析手法に合わせてプロトコール記載する面でも、入門統計学記載されているままである。 も統計的方法と軌を一にする。実際分類器では以下に挙げる点で共通点が見いだされる。
 モデルの基本は線形回帰モデルある。
 活性化関数は、回帰モデルで用いられてきた関数である。
 分類問題が多項分布   の出現確率の推定を通じて分類法が構成される。
の推定は自然母数  !   を推定する一般化線形回帰モデルである 。
" 母数推定の基本は    # 最小化で最尤推定法である。
$ 分類器の性能は試験データで評価される。
 手法の単機能性
複雑なデータを解析する手法は大抵探索的方法である。既に分かっている知識だけではなくてデー
タから新しい知識を得ようとする。前節で言及した手法がそうである。 は将に用法・用量から対象患者・期待する効果まで細かく分かっていることを前提とする。試験の結果は、想定通りの効果が認められるか否かの判断基準が得られるだけである。入力と出力のみを扱う、単一の機能しかない。
 は期待するタスクを解決することを目標とする。訓練データから分類器を構築して、試験データで性能を示す。それ以外に何らかの知識を得ようとする配慮は成されない。 と同様に単一の機能しかない。関連分野の研究者からは、特徴量など関連知識が得られないことが重大な欠点として批判された。同じ批判は臨床薬理学に対する薬理学からの批判、ベイズ法に対する頻度論者の批判にも見られる。
解析法の性能の良さと説明の分かり易さのどちらを優先させるかは、科学観に依存する。しかし、良い性能も示せないのにその理由が分かりやすく説明できる場合は多くないと思われる。目的を一つに絞って高い性能を実現させるという括りでは  は  と同じである。

 統計的推論から著しい成功に拘わらず現行の  は統計学的視点から見て改善したい点が多く見受けられる。
既存の統計的手法とは大きく異なるから、新しい問題が多い。また、用語と概念には曖昧な点が散見される。中でも    % &' とか       は重要な用語であるが、その意味は曖昧である。
更に、解析法の重視とデータの収集の軽視は特に顕著である。訓練データを原料として機械学習によりタスクの解法という製品を作り出すとの視点がある。実際解析の流れを図示する、データが底辺にあって解析結果が頂上にある。データの収集法について論じることは殆どない。
アカデミアでの深層学習の研究は、データ集合の個別データを収集する費用が小さい場合が多い。
そのために、データ集合を効率よく育成する視点が弱い。実際、( %   ら )  でもデータについては   # とか    が紹介されているに過ぎない。
文献

%   *+  ,$+     *   
深層学習 近代科学社，東京  岡谷 貴之 ,"+ &. 東京 (,"+ *-  +

麻生英樹

女性の経済・雇用に関する市区町村の指標と死亡率との関連
鈴木有佳
大阪医科大学医学部医学科 社会・行動科学教室【背景】
個人の社会経済状況とその健康との関連については国内外に多くの研究が存在するが、女性の社会的地位に関する地域指標と健康の関連に関する研究は限られている。米国の研究では、女性の地位に関する地域指標が、女性の死亡率ならびに男性の死亡率にも影響すると示されている。(1) 日本は、経済参加等に関して、他国と比較して男女間の差が大きいことが知られている上、国内でも、地域間の社会格差が広がっていることが指摘されており、その健康影響を明らかにする必要がある。今回、日本人女性の経済・雇用に関する
市区町村の指標とそこに住む男女の死亡率との関連について検討したので、発表する。

【方法】
女性の経済・雇用に関する地域指標については、過去に米国で実施された先行研究(1)を参考に、日本の実情に合わせ、各自治体の以下の 7 指標を選定し、公的統計情報
からそれぞれ算出した。経済的自立度に関する指標として、女性の業主割合（対就業状況にある女性人口）(2)、女性の大学卒業以上の割合（対 25 歳以上女性）(2)、女性の生活保護率（対 20 歳以上女性人口（10 万人あたり））(3)の 3 指標を、雇用・所得に関する指標として、女性の労働力人口の割合（対 15 歳以上女性人口）(2)、就業者
中の女性割合（対全就業者）(2)、女性管理・専門職割合（対全就業者）(2)、賃金男女比（産業計）（女性/男性）(4)の 4 指標を算出した。ただし、生活保護被保護率ならびに賃金平均については、市区町村データが公開されていなかったため、都道府県データを代入した。
目的変数には、平成 22 年人口動態調査(5)ならびに平成 22 年国勢調査(2)から算出
した、各市区町村の男女別全死因、循環器疾患、がん、呼吸器疾患による死亡率を用いた。対象地域は人口 100 人以上の全 1,884 市区町村、対象年齢は全年齢とした。
それぞれの自治体について、各地域指標を平均値で除した値と男女の死亡率との関連について重回帰分析を行った。尚、調整変数には各自治体の年齢中央値(2)、人口密度(2)、年間収入のジニ係数(6)、平均等価可処分所得(7)および貧困率(7)を用いた。
【結果】
重回帰分析の結果、女性の業主割合は女性総死亡率と正に相関し（β=105.6, p＜
0.001）、女性の生活保護率（β=-29.2, p=0.01）、女性の大学卒業以上の割合（β=43.9, p=0.001）、就業者中の女性割合（β=-248.6, p=0.01）および賃金男女比（β=-395.1, p=0.03）は負に相関した。循環器疾患による死亡率も同様の結果であった。さ
らに、がん、呼吸器疾患による死亡率と就業者中の女性割合および女性管理・専門職
割合はそれぞれ正に相関し、女性の大学卒業以上の割合および女性の労働力人口の割合はそれぞれ負に相関した。
また、男性においては、女性の業主割合（β=38.3, p=0.02）、女性の大学卒業以上の割合（β=-64.9, p＜0.001）、女性管理・専門職割合（β=-55.9, p=0.03）と総死亡率に関連が認められた。循環器疾患による死亡率も同様の結果であった。さらに、がんによる死亡率は女性の生活保護率と正に相関し、女性の業主割合、女性の大学卒業以上の割合および女性の労働力人口の割合と負に相関した。呼吸器疾患による死亡率と女性の業主割合および女性管理・専門職割合は正に相関した。
【結論】
各市区町村における女性の経済的自立度ならびに雇用・所得に関する指標は、女性の死亡率に加え、その地域の男性の死亡率にも関連することが明らかになった。よって、女性の就労支援を進める等、市区町村の状況を改善することが、男女の死亡率の減少に寄与する可能性が示された。
【謝辞】
本研究は、JSPS 研究活動スタート支援 18H06403 ならびに平成 28 年度厚生労働科
学研究費補助金（女性の健康の包括的支援政策研究事業：H28-女性-一般-001）の助成を受けて行われた。
本シンポジウムでの発表にあたり、ご指導を賜った谷川武先生、池田愛先生（順天堂大学）、丸山広達先生（愛媛大学）、磯博康先生（大阪大学）、本庄かおり先生（大阪医科大学）に感謝の意を表する。
References
1.

Kawachi I, Kennedy BP, Gupta V, Prothrow-Stith D. Women's status and thehealth of women and men: a view from the States. Soc Sci Med. 1999;48:21-32.
2.

総務省統計局. 平成 22 年国勢調査. 2011.

3.

厚生労働省社会・援護局保護課. 平成 22 年度被保護者全国一斉調査

被保護人員数、年齢階級・都道府県－指定都市－中核市別.
4.

厚生労働省政策統括官付参事官付賃金福祉統計室. 平成 22 年賃金構造基本統計調査.

2011.
5.

厚生労働省. 平成 22 年人口動態調査. 2011.

6.

総務省統計局. 平成 26 年全国消費実態調査. 2016.

7.

厚生労働省. 平成 22 年国民生活基礎調査. 2011.

マイクロシミュレーションモデルを用いた大腸がん検診における受診年齢上限の検討
福井敬祐 1・加茂憲一 2・伊藤ゆり 1・片野田耕太 3・中山富雄 41 大阪医科大学

研究支援センター

2 札幌医科大学医療人育成センター
3 国立がん研究センターがん対策情報センター
4 国立がん研究センター社会と健康研究センター【背景】大腸がんは、わが国で二番目に多い

開発は, 図 1 にあげた自然史モデルにおいがんであり、早期発見の場合はほとんどがて, 各コンパートメント間の移動を数理的治癒可能である。1990 年代に導入された便にモデリングすることで行った. 使用した潜血検査による大腸がん検診の死亡率減少

ソフトは R 言語 ver 3.3.1 である.

効果は複数の無作為化比較試験で示されて作成した MS モデルにおいて, 2012 年時点

おり、科学的根拠に基づくがん検診としての性・年齢階級別検診受診率を基に現実を推奨されている. 一方で, 近年, 高齢のが反映したコホート（上限なしコホート）およん検診受診者の増加に伴い、検診による死

び検診年齢上限を 65, 70, 75, 80, 85 歳と設定

亡率減少効果等の利益のみならず検診によした際の仮想コホート（上限ありコホート）

る不利益とのバランスを鑑み, 検診受診年をシミュレートした.

齢の上限に関する議論が活発になっている.

と上限ありコホートにおける大腸がん死亡

実際に，厚生労働省の「がん検診のあり方にと有害事象発生を比較し, 回避死亡および関する検討会」においては, 検診対象年齢

有害事象発生をそれぞれ検診の利益およびの議論が行われている しかし, これまで

不利益と定義し，比較検討を行った.

上限なしコホートのように利益と不利益を改めて無作為化比
較試験により定量化し 具体的な年齢上限を検討することは資源や, 倫理的な面から現実的でなく、欧米諸国で多用されている
マイクロシミュレーション（MS）を用いるべきである.
【目的】日本のデータに基づき開発された大腸がんに関する MS モデルを用いて, 便潜血検査による大腸がん検診の年齢上限の検討に活用できる資料を作成する。
【方法】2012 年時点で 30 歳である男女 100
万人の仮想的コホートを対象に, 100 歳まで加齢する MS モデルを開発した.

図 1 大腸がん自然史図 2

マイクロシミュレーションモデルによる上限年齢の検討結果（左：男性, 右：女性）【結果】不利益および利益が年齢上限の設【考察】がん対策において重要な役割を担

定によりどう変化するかをみた. 男性におうがん検診は，対象者を適切に選択するこいては, 検診年齢上限を 85 歳から 80 歳にとにより有効性が高まることが期待され引き下げた場合、不利益 1 件減少あたりのる。特に年齢に関しては,体力の衰えが顕利益減少は 1.5 人となった。これが 80 歳

著となる高齢者における不利益を考慮するから 75 歳では 2.23 人、75 歳から 70 歳で必要があることが MS により示された. 今は 13.2 人、70 歳から 65 歳では 16.1 人、後, 年齢上限の設定のために, 本研究の成65 歳から 60 歳では 14.3 人と算出され, 75

果を一つの根拠資料として活用されること

歳前後で大きく傾向が変化することが示唆を期待したい.

された.

シミュレーションモデルによる乳がん予防・検診の評価
堀 芽久美 1, 齋藤 英子 1, 大久保 亮 2, 清水 陽一 2, 小手森 綾香 2, 街 勝憲 2, 高橋 宏和 2
1. 国立がん研究センター がん対策情報センター
2. 国立がん研究センター 社会と健康研究センター
1． がん対策とシミュレーション研究
シミュレーションの利点の一つは、現実に再現不可能な事象もモデルによって再現できることである。シミュレーション研究は質問に応じたシナリオを作成することで、今欲しい情報を生みだすことができる。
米国 NCI (National Cancer Institute)は大学・病院で組織されるシミュレーションモデリ
ングネットワーク CISNET (the Cancer Intervention and Surveillance Modeling Network)を支援している[1]。CISNET では個人のイベント発生を確率的に表現し、イベントの発生や選択行動の確率表現をシナリオ化することで、実際には観察できない状況を再現し、その状況変化の影響を推計している。CISNET は Centers for Disease Control (CDC)や Agency
for Health Care Research and Quality (AHRQ)等と連携しており、結果は最適ながん対策の策定に貢献している[2-3]。
CISNET では一つの部位に対して複数のモデリンググループがあり、それぞれが独立してモデルを構築しており、部位グループごとに NCI と外部の大学または病院にコーディネーターが選任されている。NCI におけるコーディネーターの役割には、政策決定に関わる機関との連携、がんの予防、治療、スクリーニングに関する臨床研究あるいは疫学研究に関わる研究者への情報提供体制の整備が含まれる。CISNET は研究で得た情報や結果を可能な限りオープンとし、情報を共有しており、同一分野の研究者やその結果の利用者が集まるプラットフォームを形成している。コラボレーションを望む研究者に対する参加承認手続きも明確で、国際的なコラボレーションがスムーズに行える体制が整っており、シミュレーション研究の発展こそが重要な課題とされている。NCI 外部のコーディネートセン
ターは、その部位グループを構成する複数のモデリンググループ間の調整を行っている。
スクリーニングによる過剰診断などそれぞれの研究テーマに対して、必要なデータやアウ
トカムの選択など、モデリンググループ間での調整が積極的に行われる。また、CISNET
では若手研究者育成は研究グループの成果として評価する仕組みで、研究計画の初期段階から育成プランが作成されている。規模の大きい計画を長期的に進めていくために必要不可欠であり、日本でも若手育成を評価する仕組みが重要だろう。
CISNET が構築した研究手法やネットワークは、日本におけるシミュレーション研究の発展に重要な手がかりである。シミュレーションですべての事象を説明するのは不可能であるが、実験研究や観察研究とともにシミュレーション研究で得られた知見を組み合わせることが、これまでより早く、効率的に、社会が求める質問に答える近道になるだろう。

2． Health behavior of breast cancer survivor simulation (HERBS) study
乳がんと診断された後の予後改善に関して日本人を対象とした研究は少なく、我が国ではいまだ乳がんサバイバーの予後改善につながる具体的なガイドラインは未整備である。
本研究は、サバイバーシップケア（運動・栄養介入、サーベイランス、患者ケア）を行った場合、乳がんサバイバーの予後がどの程度改善するかについて、現在まで介入研究では用いられていなかったマイクロシミュレーション手法を活用することで、乳がんサバイ
バーに効果的な介入を明らかにすることを目的とする。
本研究では以下のことを明らかにする。
1) 乳がんサバイバーにおいて、がんと診断された後の運動量がその後のアウトカム（死亡率、再発率、QOL）にどのような影響を与えているか。さらに、診断時の乳がん進行
度別・年齢別に、アウトカムを改善するための最適な運動プログラムを明らかにする。
2) 乳がんサバイバーにおいて、がんと診断された後の食物繊維やイソフラボンの摂取量がその後のアウトカムにどのように影響するか。また、日本人女性の年齢・BMI別の最適な栄養摂取基準を明らかにする。
3) 乳がんサバイバーにおける診断後サーベイランスについて、利益と害の双方から検証し、最適なサーベイランスのあり方（対象年齢、サーベイランスの間隔）を明らかにする。
本研究はマイクロシミュレーション手法を活用することで、乳がんサバイバーの生存率向上・再発防止・QOL 改善のための最適な健康行動介入プログラムを具体化し、今後のガイドライン立案の基盤となることが期待される。
3． 引用文献[1] National Cancer Institute (NCI). Cancer Intervention and Surveillance Modeling
Network. Available from: https://cisnet.cancer.gov/. [Accessed December 2018].
[2] Mandelblatt JS, Stout NK, Schechter CB, van den Broek JJ, Miglioretti DL, Krapcho M,
et al. Collaborative Modeling of the Benefits and Harms Associated With Different U.S.
Breast Cancer Screening Strategies. Ann Intern Med. 2016;16;164(4):215-25.
[3] Mandelblatt JS, Cronin KA, Bailey S, Berry DA, de Koning HJ, Draisma G, et al. Breast
Cancer Working Group of the Cancer Intervention and Surveillance Modeling Network.
Effects of mammography screening under different screening schedules: model
estimates of potential benefits and harms. Ann Intern Med. 2009;151:738-47.

特異ウィシャート行列の最大・最小固有値分布について

1

東京理科大学大学院

清水康希

東京理科大学・理学部

橋口博樹

東京理科大学・理工学部

岩下登志也

はじめに
ウィシャート行列の固有値分布は，主成分分析の結果の安定性・信頼性を検証する際に重要である. 近年 DNA マ

イクロアレイデータのような標本数 n より次元数 m が大きい場合のデータ解析の必要性が高まっている．このような観点から，n < m の場合での固有値分布論の理論の構築が必要である．
非特異ウィシャート行列の固有値の正確分布や近似分布の研究は数多くある．しかし特異の場合は，ウィシャー
ト行列の逆行列が存在しないことや固有値に 0 が含まれるため，近似分布の導出には，母集団固有値にある仮定をおくなどの方法が必要である．
本研究では，Srivastava(2003) が導出した特異ウィシャート行列の固有値の同時密度関数をはじめに取り上げる．
そこでスティーフェル多様体上の積分を直交群上の積分に変換する定理を与えることで，最大固有値の正確分布及び，近似分布が導出できることを示す．

2

特異ウィシャート行列の最大固有値の正確分布
多変量正規分布 Nm (0, Σ) に独立に従う m 次元ベクトルを x1 , x2 , . . . , xn とし X = (x1 , x2 , . . . , xn ) とする．このとき，ランダム行列 W = XX ′ は自由度 n，共分散行列 Σ の特異ウィシャート分布に従う．特異ウィシャート分布に従う W を特異ウィシャート行列という．Srivastava (2003) は特異ウィシャート行列 W のゼロでない n 個の固有値 ℓ1 , . . . , ℓn の同時密度関数を次のように求めた．

( n
)
(
)
∫
2
2−nm/2 π n /2 |Σ|−n/2 ∏ (m−n−1)/2 ∏
1 −1′
f (ℓ1 , . . . , ℓn ) =ℓi(ℓi − ℓj )etr − Σ H1 L1 H1 (dH1 ) ， (1)Γn ( n2 )Γn ( m
2
Vn,m
2)
i=1

i<j

ただし，Γn (a) は多変量ガンマ関数，Vn,m = {H1 : m × n | H1′ H1 = In } となるスティーフェル多様体，(dH1 ) は

Vn,m のハール測度，L1 = diag(ℓ1 , ℓ2 , . . . , ℓn ) を表す．もし n = m ならば，Vm,m は直交行列全体の集合 O(m) である．
特異ウィシャート行列の固有値分布を導出するために，異なる型の行列を引数にもつ超幾何関数を導入する．Xを m 次対称行列，Y を n 次対称行列として，m ≥ n とする．2 変量の超幾何関数 p Fq (m,n) を次のように定義する．
p Fq

(m,n)

(α, β; X, Y ) =∞ ∑∑(α1 )κ · · · (αp )κ Cκ (X)Cκ (Y )k=0 κ∈Pnk(β1 )κ · · · (βq )κ

k!Cκ (Im )，ただし，α = (α1 , . . . , αq )，β = (β1 , . . . , βq )，Pnk = { κ = (k1 , . . . , kn )|k1 + · · · + kn = k, k1 ≥ k2 · · · ≥ kn ≥ 0}，

Cκ (X) は分割 κ の複素ゾーナル多項式，(a)κ は分割 κ のポッホハマー記号である.
Theorem 1. A を m 次対称行列，B1 = diag(b1 . . . , bn )，m 次対角行列 B = diag(b1 . . . , bn , 0 . . . , 0)，H1 は m × n行列で，m 次直交行列 H の分割 H = (H1 , H2 ) から構成されるものとする．この時，以下が成立する．
∫∫
Vol(O(m))′′
p Fq (α; β; AH1 B1 H1 )(dH1 ) =p Fq (α; β; AHBH )(dH)
Vol(Vn,m )Vol(O(m − n)) H∈O(m)H1 ∈Vn,m=n

Vol(O(m))(m,n)(α; β; A, B1 ) ．
p Fq
Vol(Vn,m )Vol(O(m − n))nm/2

ただし，Vol(Vn,m ) = 2Γnπ(m/2) .

1

Corollary 1. L1 = diag(ℓ1 , . . . , ℓn )，m 次対角行列 L = diag(ℓ1 , . . . , ℓn , 0, . . . , 0) とする．この時，次式が成立する．
(
)
)
(
∫
1 −1
Vol(O(m))
1 −1
′
(m)
etr − Σ H1 L1 H1 (dH1 ) =− Σ ,L ．
0 F0
2
Vol(Vn,m )Vol(O(m − n))2Vn,mSugiyama(1967) は，次の Lemma 1 を導出し，非特異ウィシャート行列の同時密度関数に Lemma 1 を用いることで非特異フィシャート行列の最大固有値の分布関数を導いた．また Shinozaki et al.(2018) では，正規母集団の拡張である楕円母集団下で同様の手法を用いて，楕円ウィシャート行列の最大固有値の分布関数を導出した．特異ウィ
シャート行列の場合も同時密度関数に Lemma 1 を用いることで最大固有値の分布関数が導出できることを示す．

Lemma 1. 対角行列 X1 = diag(1, x1 , . . . xn )，X2 = diag(x2 , . . . xn )，x2 > · · · > xn > 0 とする.
∫
n
n
∏
∏
∏
|X2 |t−(n+1)/2 Cκ (X1 ) (1 − xi ) (xi − xj )dxi
1>x2 >···xn >0

i=2

i<j

i=1

Γn (t, κ)Γn ((n + 1)/2)Cκ (In )= (nt + k)(Γn (n/2)/π n /2 )，Γn (t + (n + 1)/2, κ)
2

ただし，Re(t) > (n − 1)/2，Γn (α, κ) = (α)κ Γ(α)．

Theorem 2. W ∼ Wm (n, Σ)， m > n とする．このとき，W の最大固有値 ℓ1 の分布関数は次のように表される．
∞

Γn (n + 1/2)( x2 )nm/2 ∑ ∑ (m/2)κ Cκ (− 21 xΣ−1 )Cκ (In )Vol(O(m))
Vol(Vn,m )Vol(O(m − n)) Γn (n + m + 1/2)|Σ|n/2(n + m + 1/2)κ k!Cκ (Im )k=0 κ∈Pnk()
Γn (n + 1/2)( x2 )nm/2
Vol(O(m))
m n + m + 1 1 −1
(m,n)
=
;
; − xΣ , In ．
1 F1
Vol(Vn,m )Vol(O(m − n)) Γn (n + m + 1/2)|Σ|n/2222

Pr(ℓ1 < x) =3

特異ウィシャート行列の固有値分布の χ2 近似
Corollary 1 の 0 F0 (m) に対して， Butler and Wood(2005) で与えられるラプラス近似を適用する次の近似式が成

立する.

(
0 F0

(m)

ただし，cij =1− Σ−1 , L2)(
) n ( )( )1/2mnn∏∑∏ 2π 1/2 ∏2π
1Vol(O(m
−
n))
ℓ
i
n
≈2
exp −
，
Vol(O(m))
2 i=1 λi i<j cijdij
i=1 j=n+1(2)(ℓi −ℓj )(λi −λj )ℓ (λ −λ )，dij = i λii λj j である．(1)，(2) より，f (ℓ1 . . . , ℓn ) のラプラス近似を得る．
λi λjTheorem 3. 特異ウィシャート行列 W ∼ Wm (n, Σ)，n < m，母共分散行列 Σ の固有値を λ1 ≥ λ2 ≥ · · · ≥ λm > 0，n
W の固有値を ℓ1 ≥ ℓ2 ≥ · · · ≥ ℓn > 0 この時隣接固有値の比 ρ が ρ = max{ λλ12 , λλ32 , · · · , λλn−1} → 0 の時，特異ウィ

シャート行列の第 i 固有値 ℓi に対して次の近似式が成立する.

Pr[ℓi < x] ∼ χ2n−i+1

(

x
λi

)

参考文献

[1] Butler, R. W., Wood, A. T. A. (2005). Laplace approximations to hypergeometric functions of two matrix
arguments, Journal of Multivariate Analysis, 94(1), 1–18.
[2] Shinozaki, A., Hashiguchi, H., Iwashita, T. (2018). Distribution of the largest eigenvalue of an elliptical
Wishart matrix and its simulation, Journal of the Japanese Society of Computational Statistics, 19(1)45–56.
[3] Srivastava, M. S. (2003). Singular Wishart and multivariate beta distributions, The Annals of Statistics,
31(5), 1537–1560.
[4] Sugiyama, T. (1967). On the distribution of the largest latent root of the covariance matrix, The Annals of
Mathematical Statistics, 38(4), 1148–1151.

2

ウィッシャート行列と正規ベクトルの積の密度関数とモーメント
米永 航志朗

鈴川晶夫

北海道大学経済学院

北海道大学公共政策大学院

本報告はウィッシャート行列と正規ベクトルの積に関する密度関数とモーメント公式についてである. 多変量正規分布とウィッシャート分布は多変量解析の様々な場面において現れる. 例えば多変量正規分布を母集団分布として想定した場合に, 代表的な多変量解析手法である主成分分
析や判別分析において, ウィッシャート分布, あるいはその関数が現れることは周知である. また,X1 , · · · , Xn を Nk (µ, Σ) からの独立同分布な標本としたときに, ただし, Nk (µ, Σ) は平均ベクトル∑nµ, 分散共分散行列 Σ の k 次元正規分布を表しているが, 母平均 µ の不偏推定量 X̄ = n1 i=1 Xi∑n
1′は Nk (µ, Σ/n) に従い, 母分散共分散行列の不偏推定量である S = n−1i=1 (Xi − X̄)(Xi − X̄)については (n − 1)S ∼ Wk (n − 1, Σ), すなわち自由度 n − 1, 分散共分散行列 Σ のウィッシャー
ト分布となる.
以 上 の よ う に, 多 変 量 正 規 分 布 と ウ ィ ッ シ ャ ー ト 分 布 は 多 変 量 解 析 に お い て 重 要 で ある の だ が, こ れ ま で 統 計 学 的 な 文 脈 に お い て ，ウ ィ ッ シ ャ ー ト 行 列 と 正 規 ベ ク ト ル の積 に 関 す る 結 果 は Bodnar&Okhrin(2011) [1], Bodnar&Mazur&Okhrin(2013) [2] そ し て

Kotsiuba&Mazur(2016) [3] を除いてあまり議論されてこなかった. ウィッシャート行列の逆行列は逆ウィッシャート行列と呼ばれるのであるが, Bodnar&Okhrin(2011) [1] は逆ウィッシャート
行列と正規ベクトルの積について密度関数の導出をした後，ファイナンスにおけるポートフォリオ
理論への応用や判別分析に対する応用を議論しており, Bodnar&Mazur&Okhrin(2013) [2] では
ウィッシャート行列と正規ベクトルの積の確率表現と密度関数を導出した後，その密度関数の近似について議論している. また, Kotsiuba&Mazur(2016) [3] においては逆ウィッシャート行列と正規ベクトルの積の確率表現を導出した後, それの極限分布と密度関数の近似について議論している.
ここで, 確率表現という言葉は, ある確率変数と同じ分布を持つような別の確率変数によってその確率変数を表すという意味合いであると考えられる.
ところで，判別分析においては通常, 逆ウィッシャート行列と正規ベクトルとの積が大切であるのだが, ベイズ統計学の観点からみた場合には, 分散共分散行列の事後分布として逆ウィッシャー
ト分布を得ることを考えると, ウィッシャート行列と正規ベクトルの積もまた重要である. 例として 2 群判別の例を考えてみたい. すなわち, 2 つの正規母集団 Nk (µ1 , Σ) と Nk (µ2 , Σ) があり, い
ま新しい観測ベクトル x = (x1 , · · · , xp ) が得られたとき, それがどちらの母集団に属するかを判別する数学的基準として線形判別関数 L(x) = (µ1 − µ2 )′ Σ−1 x を得る. そして線形判別関数の値によって新しい観測ベクトルがどちらの母集団に属するのかを決めることになる. 今, µ1 , µ2 は共に既知とすれば, ウィッシャート行列と確率ベクトルの積の 1 次結合の分布が重要となるのは明白である. さらに, 新たに得られた観測ベクトルはもとの正規母集団と同じ分散共分散行列を持っているとは限らない. 例えば, 先の 2 つの正規母集団が健常者と非健常者の母集団と考えてみたい. このとき何らかの理由で非健常者のコレステロール値の分散が変わるというケースも考えられる. このような場合には判別関数の分布がどうなるのかということに関心が出てくる. 本報告ではこのよ
うにウィッシャート行列と正規ベクトルの分散共分散行列が異なっている場合に対応した議論と平均ベクトルが既知の場合の判別関数の分布に対応した議論を行っていく.
本報告でははじめに, 異なる分散共分散行列を持つウィッシャート行列と正規ベクトルの積に関して確率表現を与えたのち, それを用いて密度関数を導出する. さらに, それらがいくつかの特殊な条件の下でどのように表されるかについて考察する. さらに分布の重要な情報を含むモーメントについて 4 次までのモーメント公式を与える. 本報告では 3 次モーメントと 4 次モーメントは確率表
現に基づいてモーメント公式が計算されているが, このうち 3 次モーメントについては条件付き期
待値の考え方を用いて別な表現を与えた. また, ウィッシャート行列と正規ベクトルの積に関する密度関数は積分を含むため, その形状を把握することができない. そこで, 数値積分によってグラフを描き, それとともにシミュレーションを行い, 数値積分によって描いたグラフと比較することで今回考察した理論的結果の正当性も確認する.

参考文献[1] Bodnar, T.& Okhrin, Y.(2011) On the Product of Inverse Wishart and Normal Distributions with Applications to Discriminant Analysis and Portfolio Theory. Scandi. J. Statist.
38, 311-331.
[2] Bodnar, T.& Mazur, S. & Okhrin, Y.(2013) On the exact and approximate distributionsof the product of a Wishart matrix with a normal vector. J. Multivariate Anal. 122, 70-81.
[3] Kotsiuba, I.& Mazur, S. (2016) On the asymptotic and approximate distributions of the
product of an inverse Wishart matrix and a Gaussian vector. Theor. Probability and Math.
Statist.93, 103-112.
[4] Muirhead, R. J(1982). Aspects of multivariate statistical theory. Wiley, New York.
[5] Gupta, A. K. & Nagar, D. K(2000). Matrix Variate Distributions. Chapman and Hall/CRC.
[6] Mathai, A. M. & Provost, S.B(1992) Quadratic forms in random variables. Marcel Dekker,
Inc, New York.

ヒルベルト空間における正規性の検定の実際的側面
千葉大・融合理工学府

牧草 夏実

1 はじめにi.i.d.

P を確率分布とするとき，データ X1 , . . . , Xn ∼ P に基づく
帰無仮説 H0 : P = N (m0 , Σ0 ) v.s 対立仮説 H1 : P ̸= N (m0 , Σ0 )の検定を “正規性の検定” という．ここで，m0 = E[X1 ]，Σ0 = V [X1 ] である．ユークリッド空間
では，1 変量，そして多変量でも様々な検定方法が議論されている．この正規性の検定の数理的な一般化を考え，ヒルベルト空間に値を取る確率変数に対する正規性の検定について考える．特に，Maximum Mean Discrepancy に基づく正規性の検定について，理論的側面と，実際的側面から議論を与える．

2 ヒルベルト空間における正規分布
ヒルベルト空間 H に値を取る確率変数 X が正規分布に従うとは，任意の h ∈ H に対して，H に

[

2

]

おける内積 ⟨X, h⟩H が 1 変量の正規分布の確率変数となることである．また，E ∥X∥H < ∞ とすると，任意の h, h′ ∈ H に対し，⟨m, h⟩H = EX [⟨X, h⟩H ]，⟨Σh, h′ ⟩H = cov (⟨X, h⟩H , ⟨X, h′ ⟩H )を満たす m ∈ H と作用素 Σ ∈ HS(H) が一意に存在し，m，Σ をそれぞれ，X の期待値，共分散
作用素といい，X の従う正規分布を N (m, Σ) と表す．ここで，HS(H) は H から H へのヒルベ
ルトシュミット作用素の空間である．

3 検定統計量の構築
正定値カーネル k ，k に対応する再生核ヒルベルト空間を H(k) とする．このとき，適当な条件の下で，H 上の 2 つの確率分布 P と N (m0 , Σ0 ) との違いは，MMD（Maximum Mean Discrepancy）∆(P, N (m0 , Σ0 )) = ||µ(P ) − µ(N (m0 , Σ0 ))||H(k)に よ っ て そ れ ら の 違 い を 見 る こ と が で き る ([1] 参 照)．こ こ で ，µ(P ) = EX∼P [k(·, X)]，µ(N (m0 , Σ0 )) = EX∼N (m0 ,Σ0 ) [k(·, X)] は確率分布 P ，N (m0 , Σ0 ) のカーネル k による埋め込み
i.i.d.

2である．H におけるデータ X1 , . . . , Xn ∼ P を用いて，∆2 = ||µ(P ) − µ(N (m0 , Σ0 ))||H(k) は

b2 =
∆

1∑
b
k(·, Xi ) − µ(N (m,
b Σ))
n i=1
n

2

H(k)

b は m，Σ の一致推定量である．さらに，µ(N (m, Σ)) が (m, Σ)によって推定される．ここで，mb ，Σ
b 2 は ∆2 の一致推定量となることが [3] によって知られている．
に関して連続であるとき，∆

本発表では，この MMD に基づく正規性の検定について理論的側面と実際的側面について議論を与えた．

b 2 の漸近挙動について報告をした．特に，帰無仮説 H0 : P =理論的結果として，検定統計量 ∆b 2 の漸近分布が，独立な自由度 1 の χ2 分布の重みつき無限和であるこN (m0 , Σ0 ) のもとで，n∆b 2 の漸近分布が，平均が 0 の正規分布にと，また，対立仮説 H1 : P ̸= N (m0 , Σ0 ) のもとでの n∆√√b2 のなっていること，局所対立仮説 P = Pn = (1 − 1/ n)N (m0 , Σ0 ) + (1/ n)Q のもとでの n∆漸近分布が，独立な非心 χ2 分布の重みつき無限和であることについて報告を行った．

b 2 による正規性の検定の実際的な挙動に関し，H = Rd ，k をガウスカーネルとした場合で調n∆b 2 の帰無分布の棄却点，および検出力について報告を行った．棄却
べた結果を報告した．特に，n∆点の推定として，従来のシミュレーションを用いた結果を報告するとともに，漸近帰無分布に基づ
く 2 つの方法との比較検討についても報告を行った．漸近帰無分布に基づく方法の一つとして，適当な定数 c > 0，r ∈ N の単一重み付きカイ 2 乗 cχ2r によって漸近帰無分布の裾を近似する方法を用いた．もう一つの方法は，あるグラム行列の固有値によって漸近帰無分布の重みを推定すること
で，漸近帰無分布を近似する方法であり，[2] でも議論されている．これらの方法を比較した結果，
従来のシミュレーションにより導出した棄却点と，漸近帰無分布に基づく方法による棄却点は，どちらもほぼ同じ値をとっており，高次元データに対しては，実行時間を考慮すると，漸近帰無分布に基づく方法によって棄却点を導出すれば実用上十分であることがわかった．
対立仮説のもとでの挙動としては，分布 P として一様分布，指数分布，Σ0 として，Σ0 = Id の場合と，適当な共分散をもつ場合の各ケースについての検出力をシミュレーションによって求めた
結果を報告した．結果としては，次元 d が大きくない場合は，どのケースでも検出力はほとんど 1に近い値をとったが，d を大きくした場合，特に一様分布に対しては，検出力はあまり大きな値を取らない結果となった．

参考文献[1] A. Gretton, K. M. Borgwardt, M. Rasch, B. Schlkopf, and A. Smola (2007). A kernel
method for the two-sample-problem. Advances in Neural Information Procesing Systems,
volume 19 of MIT Press, Cambridge.
[2] A. Gretton, K. Fukumizu, Z. Harchaoui, and B. Sriperumbudur (2009). A Fast, Consistent
Kernel Two-Sample Test. Advances in Neural Information Processing Systems volume 22.
[3] J. Kellner and A. Celisse (2015). A One-Sample Test for Normality with Kernel Methods,
arXiv preprint arXiv:1507.02904v1.

楕円対称性の必要性検定
東京理科大学理工
岩下 登志也

Karlsruher Institut für Technologie
Bernhard Klar
東京理科大学理橋口 博樹

1. はじめに楕円分布 (族) は, 多変量正規分布 の一般化であり, 重要な多変量分布のクラスである. それゆえ, 母集団として楕円分布を仮定することは, 多くの多変量解析にとって極めて重要であり, このような観点から, 採られた標本が楕円母集団から採られたものであるか否かを検定することは, 欠くことのできないものである.

X 1 , . . . , X N を p-次元 (列) 確率ベクトルとする. X = [X 1 , . . . , X N ] により対応する p × N 観測行列 とすると, 標本平均ベクトル 及び標本共分散行列 それぞれX̄ = N −1 Xj N ,(1)

S = n−1 XQN X ′ ,n = N − 1 ≥ p,(2)と表される. ここに,

QN = IN − N −1 j N j ′N .

j N = (1, . . . , 1)′ ∈ RN ,(3)

本報告では, スケール化された残差の同時分布, 即ち p × N の確率行列

W = [W 1 , . . . , W N ] = S −1/2 [X 1 , . . . , X N ] QN = S −1/2 XQN ,(4)に基づく新たな統計量を考案し, 楕円対称性の必要条件検定の新たな手法の提案をした. ここに,W i = S −1/2 (X i − X̄),i = 1, . . . , N,(5)

S −1/2 は S の対称平方根の逆行列を表す。
2. 主たる理論的結果
X 1 , . . . , X N を X ∼ ECp (0, Λ) のランダムコピー, X = [X 1 , . . . , X N ] を p × N の観測行列 (observation
matrix) とする. 次のような左球形分布 LSp×N (ϕ) のサブクラス Fp×N
Fp×N = {X(p × N ) ∼ LSp×N (ϕX ); the distribution of Xa depends on a ∈ RN only through a′ a}と Iwashita and Klar (2014) の結果により次のような結果を得た.

1(6)

Theorem 1.

X 1 , . . . , X N は独立に ECp (0, Λ) に従うとし X = [X 1 , . . . , X N ] とおく. S を (2) により定義される標本共分散行列として p × N 確率行列 Y = S −1/2 X とする. このとき

Y ∼ SSp×N (ϕY ).

(7)

さらに, (3) の定数行列 QN は rank(QN ) = n の直交射影行列であるから, KK ′ = QN , K ′ K = In を満足する N × n (実) 行列 K が存在する. これを利用して Theorem 1 を発展されると, 次の結果を得る.

Theorem 2. Let X 1 , X 2 , . . . , X N は独立に ECp (µ, Λ) 従う確率ベクトルとし, X = [X 1 , X 2 , . . . , X N ] とおく. さらに, S を (2) により定義される標本共分散行列とする. このとき, n × p 確率行列

U = K ′ QN X ′ (nS)−1/2 = K ′ X ′ (nS)−1/2(8)は Stiefel 多様体 O(n, p) 上の一様分布に従う. ここに, QN は (3) で定義された N × N 行列である.

3. 提案する検定手法と数値実験結果(k){X i }Ni=1 (k = 1, 2, . . . , m) を p-次元確率ベクトル X のランダムコピー,
[
]
(k)
(k)
′
X(k) = X 1 , . . . , X N ,S(k) = n−1 X(k) QN X(k)
,
として
′
′
Uk = K(k)X(k)(nS(k) )−1/2 ,(k)とすると, Theorem 2 の結果から, X j

n = N − 1 ≥ p,

(9)

∼ ECp (µ, Λ) (j = 1, . . . , N ; k = 1, . . . , m) ならば Uk は独立に Stiefel

多様体 O(n, p) 上の一様分布に従うことになるので Pycke (2010) が提案した円周上のの一様分布に関する検定法を応用して, 検定手順を構成した. 提案する検定法に基づいて数値実験を実行した結果, 楕円母集団の下で,提案した検定手法が有効であるが, 歪楕円母集団の下で, 検出力を有しないことが判明した. また，非正規性の指摘がされている Iris Setosa deta についても，変数選択をすることで帰無仮説を棄却できたが，すべてで帰無仮説を棄却することはできなかった．検定手順を含め, 今後の研究課題と考えている.
謝辞
本研究は JSPS 科研費 18K11198(岩下), 18K03428(橋口，岩下) の助成をうけたものです。

参考文献

Iwashita, T. and Klar, B.(2013). The joint distribution of Studentized residuals under elliptical distributions.
J. Multivariate Anal.,128, 203–209.
J.-R. Pycke. (2010). Some tests for uniformity of circular distributions powerfull against multimodal alternatives, Can. J. Statist., 38, 80–96.

2Community-level social capital, parental psychological distress, and child physicalabuse: A multilevel mediation analysis
Nobutoshi Nawa, Aya Isumi, Takeo Fujiwara
Department of Global Health Promotion, Tokyo Medical and Dental University, Tokyo,
Japan.
Abstract
Clarifying modifiable risk factors and mediators is crucial for developing a strategy toprevent child maltreatment. The purpose of this study was to investigate the association
between community-level social capital and physical abuse towards children, and the
mediating effect of parental psychological distress by multilevel mediation analyses.
We analyzed data from a population-based study of first-grade elementary school
children (6–7 years old) in Adachi Ward, Tokyo, Japan, which was conducted in 2015.
The caregivers of first-grade students from all elementary schools in Adachi Ward (N =5,355) were asked to respond to a questionnaire assessing parents’ self-reported
physical abuse (hitting and beating) and neighborhood social capital. Among them,
4,291 parents returned valid responses (response rate: 80.1%). We performed multilevel
analyses to determine the relationships between community-level parental social capitaland physical abuse, and further multilevel mediation analyses were performed todetermine whether parental psychological distress mediated the association. Low
community-level social capital was positively associated with physical abuse (bothbeating and hitting) after adjustment for other individual covariates. Multilevel
mediation analyses revealed that community-level parental psychological distress did
not mediate the association. Fostering community-level social capital might beimportant for developing a strategy to prevent child maltreatment, which may have a
direct impact on abusive behavior towards children.
ACKNOWLEDGEMENTS
We are particularly grateful to the staff members, central office of Adachi City Hall for
conducting the survey and everyone who participated in the surveys. In particular, we
would also like to thank Mayor Yayoi Kondo, Mr. Syuichiro Akiu, Mr. Hideaki Otaka,and Ms. Yuko Baba of Adachi City Hall, all of whom contributed significantly tocompletion of this study. This study was supported by a Health Labour Sciences
Research Grant, Comprehensive Research on Lifestyle Disease from the Japanese
Ministry of Health, Labour and Welfare (H27-Jyunkankito-ippan-002), andGrants-in-Aid for Scientific Research from the Japan Society for the Promotion of
Science (JSPS KAKENHI Grant Number 16H03276, 16K21669 and 16J11423), St.
Luke’s Life Science Institute Grants, and the Japan Health Foundation Grants.

子どもの貧困と虐待における親の心理的ストレスと個人レベルのソーシャル・キャピタルの媒介効果

Mediating effects of parental psychological distress and individual-level social capitalon the association between child poverty and maltreatment in Japan伊角
東京医科歯科大学

彩

国際健康推進医学分野

1. は じ め に
貧 困 は 子 ど も 虐 待 の リ ス ク と 考 え ら れ て い る が 、親 の 心 理 的 ス ト レ ス や 個 人
レ ベ ル の ソ ー シ ャ ル・キ ャ ピ タ ル（ SC）が 貧 困 と 子 ど も 虐 待 と の 関 連 を 媒 介 して い る か は 検 討 さ れ て い な い 。そ こ で 本 研 究 で は 、足 立 区 の 小 学 1 年 生 を 対 象と し た 悉 皆 調 査 を 用 い て 、 親 の 心 理 的 ス ト レ ス や SC の 媒 介 効 果 を 明 ら か に することを目的とした。

2. 方 法平 成 27 年 7 月 お よ び 11 月 に 、足 立 区 立 小 学 校 に 在 籍 す る 小 学 1 年 生 の 保 護 者
5,355 名（ 有 効 回 答 率 80.1%）を 対 象 に 実 施 さ れ た 質 問 紙 調 査「 子 ど も の 健 康 ・
生 活 実 態 調 査 」の デ ー タ を 用 い て 、ロ ジ ス テ ィ ッ ク 回 帰 分 析 お よ び ポ ワ ソ ン 回
帰 分 析 を 行 っ た 。本 研 究 で は 、①世 帯 年 収 300 万 円 未 満 、②生 活 必 需 品 の 非 所
有 、③支 払 い 困 難 経 験 の い ず れ か ひ と つ で も 該 当 す る 場 合 、生 活 困 難 に あ る 状
態 と 定 義 し た 。虐 待 傾 向 に つ い て は 、身 体 的 虐 待 、ネ グ レ ク ト 、心 理 的 虐 待 を
含 む 9 種 類 の 虐 待 の 頻 度 を 尋 ね た 。分 析 に は 、生 活 困 難 と 虐 待 傾 向 に 欠 損 の な
い 3,944 名 の デ ー タ を 用 い た 。

3. 結 果生 活 困 難 と 虐 待 は ど の 虐 待 に お い て も 関 連 が あ る こ と が 確 認 さ れ た （ 表 1） 。
さ ら に 媒 介 分 析 に よ り 、親 の 心 理 的 ス ト レ ス は 生 活 困 難 と 身 体 的 虐 待 、ま た 生
活 困 難 と 心 理 的 虐 待 と の 関 連 の 60%以 上 を 説 明 す る こ と が 明 ら か に な っ た（ 表2）。一 方 、SC は 生 活 困 難 と い ず れ か の 虐 待 の 10%し か 説 明 し て い な か っ た（ 表2）。さ ら に 構 造 方 程 式 モ デ リ ン グ に よ っ て 、親 の 心 理 的 ス ト レ ス と SC が 同 時に 媒 介 し て い る こ と が わ か っ た （ 図 1） 。

4. 考 察本 研 究 よ り 、親 の 心 理 的 ス ト レ ス が 軽 減 す る よ う に サ ポ ー ト す る こ と が 、貧 困が子ども虐待に与える影響を緩和させるのに効果的かもしれないことが示唆された。

表 1 ポワソン回帰分析およびロジスティック回帰分析の結果
Crude
Prevalence of
Outcome: Any Maltreatment
95% CI
maltreatment (%) IRR
Ref
37.1
No
Poverty
1.29 1.15–1.44
47.8
Yes
0.98 0.97–0.99
N/A
Age
0.91 0.85–0.97
N/A
Education
1.05 1.04–1.06
N/A
K6 total score1.14 1.07–1.21
N/A
SC mean score
Outcome: Physical Abuse
95% CI
OR
Ref
10.8
No
Poverty
1.67 1.36–2.05
16.8
Yes
0.95 0.93–0.97
N/A
Age
0.74 0.64–0.84
N/A
Education
1.10 1.09–1.13
N/A
K6 total score1.39 1.24–1.57
N/A
SC mean score
Outcome: Neglect
95% CI
OR
Ref
12.2
No
Poverty
1.67 1.38–2.04
18.8
Yes
0.99 0.97–1.01
N/A
Age
0.97 0.86–1.09
N/A
Education
1.17 1.05–1.09
N/A
K6 total score1.30 1.16–1.45
N/A
SC mean score
Outcome: Psychological Abuse
95% CI
IRR
Ref
28.8
No
Poverty
1.34 1.18–1.51
38.5
Yes
0.98 0.97–0.99
N/A
Age
0.90 0.84–0.97
N/A
Education
1.06 1.05–1.07
N/A
K6 total score1.15 1.07–1.24
N/A
SC mean score
IRR: Incidence Rate Ratios, OR: Odds Ratios, Bold: p<0.05

Model 1
95% CI
IRR
Ref
1.23 1.10–1.38
0.99 0.98–1.00
0.94 0.88–1.01

OR
Ref
1.43
0.96
0.80

OR
Ref
1.69
1.00
1.04

IRR
Ref
1.27
0.98
0.95

95% CI
1.16–1.78
0.94–0.98
0.70–0.91

95% CI
1.38–2.08
0.98–1.01
0.92–1.18

95% CI
1.12–1.44
0.97–0.99
0.88–1.02

Model 2
95% CI
IRR
Ref
1.08 0.96–1.22
0.99 0.98–1.00
0.96 0.90–1.03
1.15 1.04–1.06
95% CI

OR
Ref
1.10
0.96
0.83
1.10

0.88–1.38
0.94–0.98
0.72–0.95
1.08–1.12
95% CI

OR
Ref
1.44
1.00
1.07
1.07

1.27–1.79
0.98–1.02
0.95–1.22
1.05–1.09
95% CI

IRR
Ref
1.09
0.98
0.97
1.06

0.95–1.24
0.97–0.99
0.90–1.05
1.04–1.07

Model 3
95% CI
IRR
Ref
1.21 1.08–1.35
0.99 0.98–1.00
0.95 0.89–1.02
1.11
OR
Ref
1.35
0.96
0.82

1.04–1.18
95% CI
1.09–1.68
0.94–0.98
0.72–0.94

1.30
OR
Ref
1.62
1.00
1.07

1.15–1.47
95% CI

1.25
IRR
Ref
1.24
0.98
0.96

1.12–1.41
95% CI

1.12

1.04–1.20

1.32–1.99
0.98–1.01
0.94–1.21

1.09–1.41
0.97–0.99
0.89–1.04

Model 4
95% CI
IRR
Ref
1.07 0.95–1.21
0.99 0.98–1.00
0.97 0.90–1.04
1.05 1.04–1.06
1.05 0.98–1.12
95% CI
OR
Ref
1.07 0.85–1.35
0.96 0.94–0.98
0.84 0.73–0.96
1.19 1.07–1.12
1.18 1.04–1.33
95% CIOR
Ref
1.41 1.14–1.75
1.00 0.98–1.02
1.09 0.96–1.24
1.06 1.04–1.08
1.18 1.05–1.32
95% CI
IRR
Ref
1.08 0.94–1.23
0.98 0.97–0.99
0.97 0.90–1.05
1.05 1.04–1.07
1.05 0.98–1.13

表 2 媒介分析結果
CDE (IRR)

95% CI

NIE (IRR)

95% CITE (IRR)

95% CI

Percent
mediated

1.08
1.21

0.99 ―1.18
1.11 ―1.32

1.12
1.02

1.09 ―1.15
1.01 ―1.03

1.21
1.23

1.11 ―1.32
1.13 ―1.34

59.1
9.6

CDE (OR)

95% CI

NIE (OR)

95% CITE (OR)

95% CI

Percent
mediated

Psychological distress
SC

1.10
1.35

0.89 ―1.39
1.10 ―1.68

1.25
1.05

1.18 ―1.32
1.03 ―1.09

1.38
1.43

1.11 ―1.73
1.16 ―1.77

69.9
14.8

Outcome: Neglect

CDE (OR)

95% CI

NIE (OR)

95% CITE (OR)

95% CI

Percent
mediated

Psychological distress
SC

1.44
1.62

1.17 ―1.79
1.31 ―1.98

1.16
1.05

1.11 ―1.22
1.02 ―1.08

1.67
1.69

1.35 ―2.05
1.37 ―2.07

28.7
8.6

95% CI

NIE (IRR)

95% CITE (IRR)

95% CI

Percent
mediated

1.24
1.27

1.11 ―1.37
1.14 ―1.41

60.2
9.2

Outcome: Any Maltreatment
Psychological distressSC
Outcome: Physical Abuse

Outcome: Psychological Abuse CDE (IRR)

Psychological distress
1.09
0.98 ―1.21
1.14
1.11 ―1.17
SC
1.24
1.11 ―1.38
1.02
1.01 ―1.04
CDE: Controlled Direct Effects, NIE: Natural Indirect Effects, TE: Total Effects
IRR: Incidence Rate Ratios, OR: Odds Ratios, Bold: p<0.05

Social capital

.20

.15

Householdincome.73
Material
deprivationChildmaltreatment.09.64
Experience ofbeing unable
to pay

.30

.07

.30

Poverty.40

Psychologicaldistress.60
Physicalabuse.34
Neglect

図 1 構造方程式モデリングによる結果.62

PsychologicalAbuseがん登録資料を用いたがん患者の生存時間解析
大阪医科大学

研究支援センター

医療統計室

伊藤

ゆり

はじめに
がん登録はがんの発生状況とその予後を把握し、がん対策に活用するために、わが国では 1950 年代から地域がん登録として宮城県から開始した。地域において発生したがんを全て登録する悉皆調査である。これにより、がんの罹患率（人口 10 万人あたりの年間がん発生率）が把握される。登録されたがん患者情報は人口動態統計および住民基本台帳（住民票）と照合され予後が把握される。これにより、がん患者の生存率が把握可能となる。本報告では、がん登録資料を活用したがん患者の生存時間解
析の各手法について、適用事例を交えて紹介する。

図 1. がん過剰死亡の考え方
長期生存率の計測：Period analysis
がん患者の長期生存率を計測する際、通常の方法（Cohort 法）では、同一期間に診断された患者コホートを追跡し、累積生存率を算出する。その場合、10 年、20 年
生存率と長期予後を知りたい場合、それだ

Competing risk の扱い：Net survival, excess

け古くさかのぼったがん患者の情報に基づ

hazard model

くこととなり、得られた生存率は outdate

がん患者の生存時間解析におけるイベン

なものとなる。簡易生命表により平均余命

ト発生は当該がん死亡であることが理想とを推定する方法と同じ方法を用いて、長期されるが、その正確な把握は困難である。

生存率を推計するのが period analysis でがん登録資料を用いた生存時間解析においある。最新の医療の状況を反映した長期生ては、全死因の死亡を用いることとなる存率が提供可能となるとして、わが国でもが、当該がん死亡以外の死亡の影響を極力適用が始まった。

除去したいと考える。そこで、対象集団の期待死亡を生命表に基づいて算出し、それよりも過剰な死亡（excess mortality,
excess hazard）を当該がん死のイベントととらえる考え方を適用している（図1）。それにより推定される生存率を net
survival と呼ぶ。また多変量解析では、
Poisson regression model を適用した
excess hazard model を使用する。

新しい予後指標：Conditional survival
通常、がんの診療の場面で参照される生存率は 5 年生存率であり、患者は診断から5 年間診断時に提示されたがん生存率を参照し続けることとなる。5 年が治癒の目安と為れてきたためである。しかし、がんの長期生存が可能となった今、患者の社会復帰において、診断からの経過年数に応じた
その後の予後に関する指標が求められている。上述の period analysis を適用した 10
年生存率を用いて、診断から 1 年以上生存した患者のその後の 5 年生存率（1 年以上生存者の 6 年生存率に該当）、2 年以
上・・・と、X 年以上生存した者の X+5
年生存率を Conditional 5-year survival として計測し、患者や医療現場で活用してもらう試みをはじめた（図 2）。

した。
予後予測モデル：樹木構造接近法（CART）
生存時間解析における樹木構造接近法はSegal や LeBlanc and Crowley らにより紹介された。これをがん登録資料に適用するために、Competing risk を考慮した生存解析に拡張した（杉本. 2013）。がん専門
病院における院内がん登録資料と診療科データベースをリンケージし、非小細胞肺がん患者の予後予測モデルを構築した。診断から 10 年間予後追跡された 1,612 名の症例データを用いて、相対生存に基づく樹木
構造接近法を適用したところ、12 の予後の異なるグループに分類された（図 3）。

図 2．Conditional 5-year survival
治癒モデル：Cure fraction model
過剰死亡という考え方を用いると、過剰死亡が起こらなくなったときがすなわち一般集団と死亡リスクが同じになった時点と考えられ、統計的な治癒（PopulationCure）として推定可能になる。非治癒患者の生存関数にワイブル関数や経験的な分布を当てはめて、治癒割合および非治癒患者の median survival time を推定する。
この二つの推定値のトレンドを用いて、がん医療の評価を行った。例えば卵巣がんでは、非治癒患者の生存時間が延長するだけ
でなく、標準的化学療法がガイドラインにより適用されてから治癒割合が大きく向上

図 3. 非小細胞肺がんの回帰樹木

結びにかえて
がん登録資料を活用した生存解析の結果
は、がん患者、医療現場で役立てるために強く求められているが、わが国におけるその適用事例の報告はまだ少ない。また、統計的な側面での発展性も見込まれるため、多くの研究者の参入が期待される。臨床・
統計・疫学の研究者が協働してさまざまな課題解決に取り組みながら社会のニーズに合致した統計資料を発信していきたい。

」㞧࡞ᾏὒ⏕ែࢩࢫࢸ࣒࡟࠾ࡅࡿᅉᯝᛶ᥎
ࢯࣦࣝ࢓ࣥຍ⸨ẚ࿅Ꮚ6DP 6XEEH\
0DULQH PDPPDOV UHVHDUFK JURXS ,QVWLWXWH RI0DULQH 5HVHDUFK %HUJHQ 1RUZD\㏆ᖺࡢᆅ⌫ ᬮ໬࡟క࠺ᾏὒ⏕ែࢩࢫࢸ࣒࡜ᾏὒ⏕≀✀ࡢ㈨※㔞ࡢኚ໬㸪⏕≀✀࡜⎔ቃせᅉ࡜ࡢᅉᯝ㛵ಀࢆண ࡍࡿࡇ࡜ࡣ㸪Ꮫ⾡ⓗ࣭ၟᴗⓗ࡟㔜せ࡛࠶ࡿ㸬ࣂ࣮ࣞࣥࢶᾏࡢ⁺ᴗ⟶⌮࡟㛵ࡍࡿᾏὒㄪᰝ࡛ࡣࠊᾏὒ⏕≀✀ࡸ⎔ቃせᅉࢆྵࡵࡓ㧗ḟඖࡢ᫬⣔ิࢹ࣮ࢱࡀほ ࡉࢀࡿ㸬ࡑࢀࡽࡢࢹ࣮ࢱࢆඖ࡟ࠊ,QWHUQDWLRQDO&RXQFLO IRUWKH([SORUDWLRQ RIWKH6HD ,&(6 ࡢ࣮࣡࢟ࣥࢢࢢ࣮ࣝࣉࡀ
⤫ྜⓗ࡞ᾏὒ⏕ែࢩࢫࢸ࣒ᰝᐃࢆ࠾ࡇ࡞ࡗ࡚࠸ࡿ㸬ẖᖺࡢ࣏࣮ࣞࢺ࡟ࡣ ,&(6ࡀ᥎ዡࡍࡿከኚ㔞ゎᯒ࡟ᇶ࡙࠸ࡓศᯒ⤖ᯝࡀࡲ࡜ࡵࡽࢀ࡚࠸ࡿ㸬ࡑࢀࡽࡢᡭἲࡣ᫬㛫㍈ୖࡢࢲ࢖ࢼ࣑ࢡࢫࢆ⪃៖ࡋ࡞࠸ࡓࡵ࡟ࠊᾏὒ࢚ࢥࢩࢫࢸ࣒࡟㛵ࡍࡿ⏕≀✀࡜⎔ቃせᅉࡢᅉᯝ㛵ಀࡣ᥈ࡿࡇ࡜ࡣ࡛ࡁ࡞࠸㸬ࡑ ࡇ ࡛ ࠊ2]DNL㸦  ࡀ ᥦ ᱌ ࡋࡓ㸪*UDQJHU  ࡜ *HZHNH  ࡢ ࣌ ࢔࣡࢖ࢬᅉᯝᛶ᥎
3DUWLDOSDLUZLVH FDXVDOLW\ ࡜㉥ụ  ࡢ┦ᑐࣃ࣮࣡ᐤ୚⋡㸦7RWDOFDXVDOLW\㸧ࢆ⤫ྜࡋࡓᅉᯝᛶ᥎
ࢆኚ㔞㛫┦஫㛵ಀࢆ᥎ ࡍࡿ᪉ἲ࡜ࡋ࡚㐺⏝ࡋࡓ㸬ほࡉࢀࡓ N ḟඖ᫬⣔ิࢹ࣮ࢱ [W[ W  [ W   [N W  W   1 ࡜ࡍࡿ㸦ࡇࡇ࡛㸪  ࡣ㌿⨨ࡢグྕ㸧㸬ࡇࢀࡽࡢࢹ࣮ࢱࡣ㸪௨ୗ࡟♧ࡍࡼ࠺࡞ከኚ㔞⮬ᕫᅇᖐ㸦0XOWLYDULDWHDXWRUHJUHVVLYH
0$5㸧㐣⛬ࡢᐇ⌧್࡜௬ᐃࡍࡿ㸸[W¦ P  $P[W P  İW 0

ࡇࡇ࡛ 0 ࡣ⮬ᕫᅇᖐࡢḟᩘ㸪 $ P ࡣ⮬ᕫᅇᖐࣔࢹࣝࡢಀᩘ㸪ࡑࡋ࡚ İ W ࡣᖹᆒࢮࣟ࣋ࢡࢺࣝ㸪ศ ᩓඹศᩓ⾜ิ 6 ࡟ᚑ࠺ከኚ㔞ṇつศᕸ࡟ᚑ࠺࡜ࡍࡿ㸬⮬ᕫᅇᖐಀᩘࡣ᭱ᑠ஧஌ἲࡸ <XOH
:DONHU ἲ࡞࡝࡛᥎ᐃࡉࢀࠊண ㄗᕪ⣔ิ࡟ࡼࡾඹศᩓ⾜ิࡀィ⟬࡛ࡁࡿ㸬⮬ᕫᅇᖐಀᩘ࡜ࣇ࣮
࢚ࣜኚ᥮࠿ࡽ㸪࿘Ἴᩘ I ࡟ᑐࡍࡿ࿘Ἴᩘᛂ⟅㛵ᩘ ) I ࡀồࡲࡾ㸪௨ୗࡢࡼ࠺࡞ࣃ࣮࣡ࢫ࣌ࢡࢺ
ࣝࡀồࡲࡿ㸸

3I) I 6) I§ S I S I  SNI ·¨¸¨ S I S I  S NI ¸¨¸   d I d '¨¸
¨ S S  S ¸
NNI ¹© N I N  Iࡇࡇ࡛ )I ࡣ ) I ࡢඹᙺ」⣲ᩘ㸪 ' ࡣほ್ࡢࢧࣥࣉࣜࣥࢢ㛫㝸࡜ࡍࡿ㸬 3 I ࡢ㠀ᑐゅᡂศࡣࢡࣟࢫࣃ࣮࣡ࢫ࣌ࢡࢺ࡛ࣝ࠶ࡿ㸬ࡶࡋ㸪 6 ࡢ㠀ᑐゅᡂศࡀ㠀ᖖ࡟ᑠࡉ࠸㸪ࡍ࡞ࢃࡕྛኚ㔞ࡢணㄗᕪࡣ⊂❧࡛࠶ࡿࡇ࡜ࡀ௬ᐃ࡛ࡁࡿ࡜㸪 L ␒┠ࡢኚ㔞 [L ࡢࣃ࣮࣡ࢫ࣌ࢡࢺࣝࡣ௚ࡢኚ㔞 [ M ࠿ࡽࡢ࿘Ἴᩘᛂ⟅㛵ᩘ )LM I ࡜ண ㄗᕪࡢศᩓ V MM ࡢᙳ㡪ࢆྵࡴ㡯ࡢ࿴

SLL I







)L I V     )LL I V LL    )LN I V NNࡢᙧ࡛♧ࡍࡇ࡜ࡀ࡛ࡁ㸪 M ␒┠ࡢኚ㔞 [ M ࠿ࡽࡢࣀ࢖ࢬࡢᙳ㡪 ULM I ࢆ௨ୗࡢࡼ࠺࡟♧ࡍࡇ࡜ࡀ࡛
ࡁࡿ㸸


ULM I
ண

)LM I V MM
SLL I

 >  @ㄗᕪ V  ࡜ࣃ࣮࣡ࢫ࣌ࢡࢺࣝ S I ࡜ࡢ㛵ಀࢆ⏝࠸ࡿ࡜ࠊࡇࡢࣀ࢖ࢬᐤ୚⋡ࡣ௨ୗࡢࡼ࠺

࡟ࠊ඲࡚ࡢኚᩘࡀ㛵ಀࡋ࡚࠸ࡿࣔࢹࣝ࡜ࠊ඲࡚ࡢኚᩘ࠿ࡽ  ࡘࡔࡅኚᩘࢆእࡋࡓࣔࢹࣝࡢண
ㄗᕪࡢᕪࡢᙧ࡜ࡋ࡚ᣑᙇࡍࡿࡇ࡜ࡀ࡛ࡁࡿ㸸

ORJ V LA M  ORJ V LL

³



 

ORJ SLL M I GI  ³



 



ORJ SLL I GI

SLL M I
ORJ
³  SLL I GI


³



 

ORJ  ³

D LM I V MM
SLL I



 

ORJ

SLL I  D LM I V MM
SLL I

GI

GI

ᡃࠎࡣࡉࡽ࡟㉥ụ᝟ሗ㔞つ‽࡟ᇶ࡙ࡃᯟ⤌ࡳ࡛㸪ࡇࡢᅉᯝ㛵ಀࡢ᭷ពᛶ࡟㛵ࡍࡿつ‽ࢆຍ࠼ࡓ㸬 ኚ㔞࡜  ኚ㔞ࡢ᫬⣔ิࢹ࣮ࢱࢆ⏕ᡂࡋࠊ*UDQJHU ࡜ *HZHNH ࡢᅉᯝᛶศᯒ࡜ࡢẚ㍑ࠊࢧࣥࣉࣝ
ᩘࡢ㐪࠸࡟ࡼࡿ 6HQVLWLYLW\ࠊᅉᯝᛶࡢᙉᗘ࡟㛵ࡋ࡚ࠊᥦ᱌ᡭἲࡢࣃࣇ࢛࣮࣐ࣥࢫࢆ᳨ドࡋࡓ㸬ࣂ࣮ࣞࣥࢶᾏࡢ⏕ែࢩࢫࢸ࣒࡟࠾ࡅࡿ㣗≀㐃㙐࡛㔜せ࡞  ✀㢮ࡢᾏὒ⏕≀㸦ࢩࢩࣕࣔ㸪ࢱࣛ㸪࢜࢟࢔࣑㸪ࢽࢩࣥ㸧ࡢࣂ࢖࣐࢜ࢫ᫬⣔ิࢹ࣮ࢱࢆ⏝࠸㸪ࢩࢩࣕࣔࡢᖺ㱋ẖ࡟⏕≀✀㛫ࡢࣇ࢕࣮ࢻࣂࢵࢡ㛵ಀࢆ⪃ᐹࡋࡓ㸬㸫 ᖺ㱋ࡢࢩࢩࣕࣔࢆྵࡴࣇ࢕࣮ࢻࣂࢵࢡࢩࢫࢸ࣒ࡣ㸪 ᖺ㱋ࡶࡋࡃࡣ  ᖺ㱋ࡢࢩࢩࣕࣔࢆྵࡴࣇ࢕࣮ࢻࣂࢵࢡࢩࢫࢸ࣒ࡼࡾࡶࡼࡾከࡃࡢ⏕≀✀㛫ࡢ┦஫㛵ಀࢆ♧ࡋࡓ㸬ࡑࢀࡽࡣࣂ࣮ࣞࣥࢶᾏࡢ㣗≀㐃㙐࡟㛵ࡍࡿࡇࢀࡲ࡛ࡢඛ⾜◊✲ࢆ⿬௜ࡅ㸪ࢩࢩࣕࣔࡀ
ࣂ࣮ࣞࣥࢶᾏᇦ⏕ែࢩࢫࢸ࣒ࡢ㣗≀㐃㙐࡟㛵㐃ࡍࡿ⏕≀✀㛫ࡢ㔜せ࡞㥑ື※࡟࡞ࡿࡇ࡜ࢆ᫂ࡽ
࠿࡟ࡋࡓ㸬ᮏ✏ࡀᥦ᱌ࡍࡿᡭἲࡣ㸪ᾏὒᏛ◊✲࡛ᑐ㇟࡜࡞ࡿ」㞧࡞⏕ែ⣔࡟࠾࠸࡚㸪⏕≀✀㛫㸪⎔ቃᅉᩘ㛫ࡢᅉᯝ㛵ಀࢆ᥎ ࡍࡿ୍ᡭẁ࡜ࡋ࡚᭷⏝࡛࠶ࡿ࡜⪃࠼ࡽࢀࡿ㸬ཧ⪃ᩥ⊩$NDLNH +  2QWKH XVHRIDOLQHDU PRGHO IRUWKH LGHQWLILFDWLRQ RIIHHGEDFNV\VWHPV $QQDOVRIWKH,QVWLWXWHRI6WDWWLVWLFDO0DWKHPDWLFV*HZHNH-  0HDVXUHPHQW RIOLQHDU GHSHQGHQFH DQGIHHGEDFNEHWZHHQ PXOWLSOH WLPH VHULHV-RXUQDORIWKH$PHULFDQ6WDWLVWLFDO$VVRFLDWLRQ*UDQJHU &:-  ,QYHVWLJDWLQJ FDXVDO UHODWLRQV E\HFRQRPLHWULF PRGHOV DQGFURVVVSHFWUDOPHWKRGV (FRQRPHWULND
2]DNL 7  7LPH VHULHV PRGHOLQJ RIQHXURVFLHQFH GDWD&KDSPDQ +DOO&5& %RFD5DWRQSummarization Method Based On Data Fusion
Mika Sato-Ilic
Faculty of Engineering, Information and Systems, University of Tsukuba,
Tsukuba, Japanmika@risk.tsukuba.ac.jp
1. Introduction
The recent growth of data requires a proportional increase in speed of analysis, aswith the problem of how to predict a future value when multiple data sets areobserved in real time. At this time, how to summarize the multiple data sets withdata fusion is important. Therefore, this presentation discusses a data fusion
method as a summarization method based on our recent researches [1], [2]. The
research focused on how to obtain comparable estimates of predicting values ofdependent variables in regression analysis when we observe multiple different data
sets simultaneously. These data sets have various forms. For example, if we
observe high-dimension and low-sample size (HDLSS) data and the data isclassified into several subclasses by using an adaptable classification of the highdimension, then the obtained subclasses are one of the multiple data sets.
The main role of ordinary linear regression analysis is to obtain the estimate ofpredicted values as projected values in a linear subspace spanned by vectors of
independent variables. However, if the data set has been observed simultaneously
from multiple different data sources, then we must create different linear
subspaces to estimate the different predicted values corresponding to the different
datasets. So, we cannot compare the different predicted values, since the linear
subspaces are different.
In order to solve this problem, we have proposed a method to obtain comparable
predicted values obtained from different datasets by utilizing a fuzzy clusteringresult and an orthogonal projector which projects two different vectors
corresponded with the two different dependent variables to the same intersectionof the two different linear subspaces. [1], [2] From this, since the different predictedvalues from different data sources can be obtained in the common space, we cancompare the different predicted values.
2. Comparable Predicted Values in Regression Analysis
Two n × p matrices 𝐴 and 𝐴 which consist of independent p different variables
( )
( )
( )
( )
( )
are shown as: 𝐴 = 𝒂 , ⋯ , 𝒂 , 𝒂 = 𝑎 , ⋯ , 𝑎 ′, 𝑟 = 1, ⋯ , 𝑝, 𝑡 = 1,2. 𝒃 and𝒃 are vectors of dependent variables, 𝜷 and 𝜷 are vectors of regression
coefficients, and 𝒆 and 𝒆 are error vectors whose elements follow the same
normal distribution. We utilize a result of fuzzy clustering for 3-way data [3] as acommon scale of 𝐴 and 𝐴 , since in this method, the same K clusters areobtained through 𝐴 and 𝐴 when T = 2. Then, we obtain the following result of
( )
( )
( )
( )
( )
the fuzzy clustering: 𝑈 = 𝒖 , ⋯ , 𝒖 , 𝒖 = 𝑢 , ⋯ , 𝑢 ′, 𝑘 = 1, ⋯ , 𝐾, 𝑡 = 1,2.
( )( )We define the following criterion:𝑐 ≡ ∑𝑎𝑏𝑠(𝑐𝑜𝑟(𝒖 , 𝒂 )), 𝑙 = 1, ⋯ ,2𝑝, 𝑡 = 1,2,
where cor(𝐚, 𝐛) means the correlation between 𝐚 and 𝐛, abs(∗) means an absolute( )( ) ( )
( )( )value of ∗ , and {𝒂 , ⋯ , 𝒂 , 𝒂 , ⋯ , 𝒂 } ≡ {𝒂 , ⋯ , 𝒂 , 𝒂 , ⋯ , 𝒂 }. Then 𝑐shows a degree of contribution of 𝒂 to t-th data 𝐴 , 𝑡 = 1,2. Therefore, we can( )( )recreate the independent vectors for each t as follows: 𝐴 = (𝒂 , ⋯ , 𝒂 ),( ) ( )( )( )( ) ( )
{𝒂 |𝑐 > 𝜀} , 𝑗 = 1, ⋯ , 𝑠, 𝐴 = (𝒂 , ⋯ , 𝒂 ), {𝒂 |𝑐 > 𝜀} , 𝑗 = 1, ⋯ , 𝑤, 𝑠 +𝑤 = 2𝑝, where ε is a given positive threshold value which is selected as at least
one same vector is included for both 𝐴 and 𝐴 . Then we can calculate orthogonal
projectors 𝑃 and 𝑃 using 𝐴 and 𝐴 as: 𝑃 ≡ 𝐴 (𝐴 ′𝐴 ) 𝐴 ′, 𝑡 = 1,2, which projecta vector to two different linear subspaces, 𝑉 and 𝑉 . In order to obtain a common
subspace of 𝑉 and 𝑉 , we utilize an orthogonal projector, 𝑃 ∩ , which projects a
vector to the intersection of subspaces 𝑉 and 𝑉 , denoted as 𝑉 ∩ . Then 𝑃 ∩ is
obtained as follows by using 𝑃 and 𝑃 [4], [5], [6].
𝑃 ∩ = 2𝑃 (𝑃 + 𝑃 ) 𝑃 = 2𝑃 (𝑃 + 𝑃 ) 𝑃 ,where (∗) shows the Moore-Penrose generalized inverse matrix of ∗ . We obtainthe following equations:𝒃 = 𝑃 ∩ 𝒃 , 𝑡 = 1,2.
Since the estimates 𝒃 and 𝒃 are obtained in the same common linear subspace,𝑉 ∩ , we can compare 𝒃 and 𝒃 mathematically.
References[1] Sato-Ilic, M., Knowledge-based Comparable Predicted Values in Regression
Analysis, In C. H. Dagli (ed) Procedia Computer Science, Elsevier, 114, 216223 (2017)[2] Sato-Ilic, M., Cluster-Scaled Regression Analysis for High-Dimension and LowSample Size Data, Advances in Smart Systems Research, 7, 1, 1-10 (2018)[3] Sato-Ilic, M., Individual Compositional Cluster Analysis, In C. H. Dagli (ed)Procedia Computer Science, Elsevier, 95, 254-263 (2016)[4] Anderson, W. N., and R. J. Duffin., Series and Parallel Addition of Matrices,
Journal of Mathematical Analysis and Applications, 26, 576-594 (1969)[5] Ben-Israel, A., and T. N. E. Greville, Generalized Inverses, Theory andApplications, 2nd ed., Springer (2003)[6] Krafft, O., An Arithmetic-Harmonic-Mean Inequality for Nonnegative Matrices,
Linear Algebra and Its Applications, 268, 243-246 (1998)

Plug-in optimization method for generalized ridge regression
for MLE in GMANOVA model
中京大学 国際教養学部 永井 勇
各個体に対して経時的に測定することで得られるデータは経時測定データと呼ばれ, 様々な分野で収集・分析がされている. この経時測定データの分析においては, データに潜む経時的な変動 (経時変動) を上手く捉えることが一つの目的である. 本研究では, 全ての個体で測定時点が揃っている, つまり, n 個全ての個体に対して時点 t1 , . . . , tp で測定した経時測定データを考えた. このとき, Potthoﬀ and Roy (1964) により提案された次の GMANOVA(一般化多変量分散分析) モデルを用いることで経時変動が推定できることが知られている;
Y = 1nµ′X ′ + AΞX ′ + E,
ここで, Y は各行が各個体の経時測定データからなる n × p 行列, 1n は全ての成分が 1 の n次元ベクトル, µ は q 次元未知ベクトル, X は p × q 既知行列 (詳しくは後述), A は各行が各個体の性別などからなる n × k 既知 (説明変数) 行列, Ξ は p × q 未知行列, E は n × p 誤
差行列である. 本研究では, rank(X) = q, A は rank(A) = k で中心化されている (つまりA′ 1n = 0k , 0k は k 次元ゼロベクトル) とし, E の各行が独立に Np (0p , Σ) に従う確率変数ベ
クトル, Σ は p × p 未知正定値行列, n − k − p + q − 2 > 0 とした. ここで, (1n µ′ X ′ + AΞX ′ )の部分が Y の経時変動を表しており, 例えば X の j 行目を (1, tj , . . . , tjq−1 ) とすることは,
経時変動を測定時点 t1 , . . . , tp の (q − 1) 次の多項式で推定することに対応している.
このモデルにおいて, µ や Ξ の最尤推定量はそれぞれ以下となる;{µ̂ = (X ′ S −1 X)−1 X ′ S −1 Y ′ 1n /n,Ξ̂ = (A′ A)−1 A′ Y S −1 X(X ′ S −1 X)−1 ,ここで S = Y ′ {In − 1n 1′n /n − A(A′ A)−1 A′ }Y /(n − k − 1) である. このとき, A′ 1n = 0k とE の各行への正規分布などの仮定より, (n−k −1)S ∼ Wp (n−k −1, Σ) (自由度が (n−k −1)で尺度行列が Σ の p 次元ウィッシャート分布, Schott (2017) Sec. 11.7 など参照) である.
この推定量において, A′ A の固有値 d1 , . . . , dk の中に小さな値が現れると Ξ̂ が不安定になってしまう. そこで, Yanagihara, Nagai and Satoh (2009), Nagai (2011) などで用いられている多変量一般化リッジ回帰による推定法を用いた次の推定量を考えた;

1Ξ̂(θ)= Mθ−1A′Y S −1X(X ′S −1X)−1,n
ここで, Mθ = A′ A + Qdiag(θ)Q′ , Q は Q′ A′ AQ = diag(d1 , . . . , dk ) となる直交行列,

θ = (θ1 , . . . , θk )′ (θi ≥ 0) はリッジパラメータと呼ばれるパラメータである.
次にこの θ の最適化を考えた. そのために, Ξ̂(θ) を用いた予測値 Ŷ (θ) = 1n µ̂′ X ′ +AΞ̂(θ)X ′ を評価する関数として, 次で定義される予測平均二乗誤差 (PMSE) を用いた;[ [ {()()′ }]]U − Ŷ (θ) Σ−1 U − Ŷ (θ)
,
1
PMSE[
Ŷ (θ)] = EU EY trn

ここで U は Y と独立で同一の分布に従う確率変数, EB [·] は確率変数 B に関する期待値を表している. この PMSE[Ŷ (θ)] は, 今のデータ Y から構築した予測値 Ŷ (θ) と新たに Y と同じモデルから得られたデータ U との平均残差を評価している関数と見ることもできる.

本研究では, この PMSE[Ŷ (θ)] を最小にする最適なリッジパラメータ θ を求めることを考えた. このとき, U が Y と独立に同一分布に従う確率変数より PMSE[Ŷ (θ)] = EY [tr{(Y −Ŷ (θ))Σ−1 (Y − Ŷ (θ))′ }] + 2q{1 + tr(Mθ−1 M0k )} なので, 第一項の期待値を求めて θ に関する項を最小にすることで最適な θ が得られる. しかしながら, そのまま最適な θ を求めても未知の行列である Ξ と Σ が残るため, 以下の二種類の手法が考えられる;
1 既知の行列からなる EY [tr{(Y − Ŷ (θ))Σ−1 (Y − Ŷ (θ))′ }] の推定値を作り, その最小化⃝を行う (Yanagihara, Nagai & Satoh (2009), Yanagihara & Satoh (2010) など参照)
2 PMSE[Ŷ (θ)] を最小にする θ を求め, 未知の行列にそれぞれの推定量を代入する (Nagai,⃝
Yanagihara & Satoh (2012) など参照)
2 の手法を提案した. このとき, Ŷ (θ) に S −1 が必要な点や A′ 1n =本研究では Ξ̂(θ) において⃝
0k などに注意して, 分割行列の逆行列などの行列自体の性質 (Lütkepohl (1996), Koll & vonRosen (2005) など参照), Wishart 分布や二次形式の独立性の性質など (Muirhead (1984),
Siotani, Hayakawa & Fujikoshi (1985), Fujikoshi, Ulyanov & Shimizu (2010), Schott (2017)など参照) を用いて, Satoh, Kobayashi and Fujikoshi (1997) と同様に期待値を計算した.
その結果, リッジパラメータ θ に関する項に着目すると, 適当な関数 fi (θi |Ξ, Σ) を用い,k∑
PMSE[Ŷ (θ)] = [θ に依存しない項] +fi (θi |Ξ, Σ),i=1という形で展開できた. よって, PMSE[Ŷ (θ)] を最小にする最適な θ (各成分は非負) は, 各
成分ごとに非負の範囲で fi (θi |Ξ, Σ) を最小にする θi を求めればよい. つまり θ̃i (Ξ, Σ) =arg minθi ≥0 fi (θi |Ξ, Σ) を各 i で求めて, それを並べれば最適な θ が得られる. このとき,
θ̃i (Ξ, Σ) に未知の行列 Ξ と Σ が残ってしまうため, 推定量 Ξ̂ と S をそれぞれに代入した.
fi (θi |Ξ, Σ) や θ̃i (Ξ, Σ) の具体的な形などや数値実験による比較などは当日報告した.
引用文献:[1] Fujikoshi, Y., Ulyanov, V. V. & Shimizu, R. (2010). Multivariate Statistics, John Willey& Sons.
[2] Kollo, T. & von Rosen, D. (2005). Advanced Multivariate Statistics with Matrices,
Springer.
[3] Lütkepohl, H. (1996). Handbook of Matrices, John Wiley & Sons.
[4] Muirhead, R. J. (1984). Aspects of Multivariate Statistical Theory, John Willey & Sons.
[5] Nagai, I. (2011). Modiﬁed Cp criterion for optimizing ridge and smooth parameters in
the MGR estimator for the nonparametric GMANOVA model. Open J. Stat., 1, 1–14.
[6] Nagai, I., Yanagihara, H. & Satoh, K. (2012). Optimization of Ridge Parameters in
Multivariate Generalized Ridge Regression by Plug-in Methods. Hiroshima Math. J.,
42, 301–324.
[7] Potthoﬀ, R. F. & Roy, S. N. (1964). A generalized multivariate analysis of variancemodel useful especially for growth curve problems. Biometrika, 51, 313–326.
[8] Satoh, K., Kobayashi, M. & Fujikoshi, Y. (1997). Variable selection for the growth
curve model. J. Multivariate Anal., 60, 277–292.
[9] Schott, J. R. (2017). Matrix Analysis for Statistics (Third ed.), John willer & Sons.
[10] Siotani, M., Hayakawa, T. & Fujikoshi, Y. (1985). Modern Multivariate Statistical
Analysis: A Graduate Course and Handbook, American Sciences Press.
[11] Yanagihara, H., Nagai, I. & Satoh, K. (2009). A bias-corrected Cp criterion for optimizing ridge parameters in multivariate generalized ridge regression. Jpn. J. Appl.
Stat., 38, 151–172 (in Japanese).
[12] Yanagihara, H. & Satoh, K. (2010). An unbiased Cp criterion for multivariate ridge
regression. J. Multivariate Anal., 101, 1226–1238.

歪曲度のノンパラメトリック推定
千葉大学大学院理学研究院

内藤貫太

はじめに: Rd におけるサイズ N の対応のある 2 つの点群 X = [X1 · · · XN ]T と Y =[Y1 · · · YN ]T を考える．ここで，Xi と Yi は d 次元確率ベクトルである．特にその組(Xi , Yi ) は i 番目のオブジェクトに関する (X, Y ) の観測値である．
ここで議論する問題は，この Rd の 2 つの変数 X と Y の類似性，“調和度”，をどのように捉え，そしてデータである点群 X と Y に基づきどのように推定するかということである．
具体例としては例えば d = 2 として，ヒトの (右大腿長，右上腕長)(= X) と同一人の (左大腿長，
左上腕長)(= Y ) が挙げられる．この場合，腕と足の長さの左右での調
和をどのように捉えるのかを考えることになる．
本講演では，この “調和度”として，擬等角写像論で現れる歪曲度を用いることを考える．そして歪曲度のノンパラメトリック核型推定について議論する．
歪曲度: この 2 つの変数を対応づける変換 f の存在を仮定する. 変数 X と Y の調和
度はこの変換 f の性質を通して議論しよう．特に，擬等角写像論において等角写像
からの乖離を測る尺度である歪曲度と呼ばれる指標をもって調和度とすることを提
案する．
√∑
√
d
2
x = [x1 · · · xd ]T ∈ Rd に対しそのノルムを |x| = xT x =i=1 xi とする．d × d
行列 A の作用素ノルムは||A|| = sup{|Ay| : y ∈ Rd , |y| = 1}.
で定義される．A の作用素ノルムは AT A の最大固有値の平方根に他ならない. Rd
から Rd の関数

g1 (x)

g(x) =  ...  , x ∈ Rd ,gd (x)
を考える．ここで，gi は Rd から R への滑らかな関数である．g の x におけるヤコビ
行列は()∂gi (x)Dg(x) =
∂xj 1≤i,j≤dで定義され，この行列式の絶対値はヤコビアンと呼ばれ Jg (x) = |det(Dg(x))| で表す.
いま，∆ を Rd の領域とし， f を ∆ から Rd への同相写像とする. また，

κ(x) = κ(x|f ) =||Df (x)||d,
Jf (x)(1)

とする．f が K-擬等角とは， ∆ 上で κ(x|f ) ≤ K となるような K > 0 が存在することである（[1]）. 変換 f の歪曲度 d = d(f ) は

d = d(f ) = max κ(x|f ).
x∈∆

(2)

と定義される．∆ 上で常に κ(x|f ) ≥ 1 であり，κ(x|f ) ≡ 1，すなわち 1-擬等角は等角写像であることと同値である．歪曲度で調和を測るためには，κ(x) の推定が必要である．この κ(x) が各 x においての f の等角写像からの乖離を測っていることになる．
素朴な考えでは，適当な平滑化で得られる f のノンパラメトリック推定量 fˆ を用
いて b
d = d(fˆ) とすることで，歪曲度のノンパラメトリック推定量が構成される．しかしながらこれはあまり筋がよろしくなく，特に fˆ の微分を要するので，挙動の不安定性が知られている．
したがって，d に現れる κ(x)，これを構成している f のヤコビ行列 Df (x) それ自体をノンパラメトリックに推定することを考える．
ノンパラメトリック推定: 多変量回帰モデル

Yi = f (Xi ) + εi ,i = 1, ..., N,(3)を想定する．ここで，εi (i = 1, ..., N ) は i.i.d. で E[εi ] = 0，V [εi ] = Σ > O を満たす．
このモデルにおいて，局所 2 次回帰（[2]）を用いることにより，Df (x) のノンパ
ラメトリック核型推定量が構成される．κ(x|f ) はその定義から，√
d∏λ1 (x)
κ(x) =
,
(4)
λ
(x)
i
i=1
と表現される．ここで，λ1 (x) ≥ λ2 (x) ≥ · · · ≥ λd (x) > 0 は Df (x)Df (x)T の固有値である.
したがって，Df (x) のノンパラメトリック推定量から Df (x)Df (x)T の推定量が構成されるから，その行列の固有値を用いることで κ(x|f ) のノンパラメトリック推定量が構成される．
研究集会においては，このようにして得られる歪曲度のノンパラメトリック推定量の漸近挙動と数値的挙動について報告を行った．

参考文献[1] Väisälä J. Lectures on n-Dimensional Quasiconformal Mappings. Lecture Notein Mathematics 229: Springer,1971.
[2] Ruppert D and Wand MP. Multivariate locally weighted least squares regression.
Ann. Statist., 1994; 22: 1346-1370.

探索的財務ビッグデータ解析
地道 正行 ∗ , 宮本 大輔 ∗∗ , 阪 智香 ∗ , 永田 修一 ∗∗ 関西学院大学 商学部∗∗ 東京大学 大学院情報理工学系研究科
2018 年 12 月 15 日
科学研究費シンポジウム「多変量解析法における理論と応用」
広島大学理学部 B707

概要
本研究では, Bureau van Dijk (BvD) 社 1) から提供されるデータベース Osiris から抽出された世界 157 カ国の全上場企業 (一般事業会社, 上場廃止企業含む) の主要財務情報 (売上高, 営業利益, 総資産など 84 項目, 33 年分) の財
務データ (粗データ) を扱った.
このデータが収められたファイルを, コンピュータ (ソフトウェア) で利用できる形式に変換する工程を「前処理」 (preprocessing) と呼び, その整形されたファイルを R に読み込み, 実際にデータ解析が行えるオブジェクト形
式に変換する工程を, データ解析環境 R の統合開発環境である RStudio 上で行うため, Grolemund and Wickham(2016) にならって「データラングリング」(data wrangling) と呼んでいる. 前処理には, Unix コマンドとインタプ
リタを用いて, データ解析環境 R に読み込める形式に変換し, R を用いてデータの不備などの修正を行いデータを分析・解析できるファイル形式 (CSV, RDS ファイル) として出力することによって行った. また, データラングリ
ングは, 近年注目されているクラスター・コンピューティング・システム Spark と SparkR パッケージを用いて分析できるオブジェクト形式に変換した. さらに, Tukey (1977) によって提唱された探索的データ解析 (Exploratory

Data Analysis: EDA) に基づき, 可視化 (data visualization) によって得られた知見を利用して, 統計モデリング(statistical modeling) を行った.
次に, このデータに対して探索的データ解析を行った. 具体的には, 2015 年に時点を固定し, データ可視化 (datavisualization) を行った結果として, 売上高, 従業員数, 総資産のそれぞの (1 変量) 分布は極度に右に歪んだものであることがわかり, 各ペアの 2 変量分布も同様の結果となった. この結果に対して, 対数をとることによって歪みを修正したところ, 若干左に歪んでいることが分かった. この知見を利用して, Azzalini (1985), Azzalini and Capitanio(2014) によって提案された非対称分布 (非対称正規分, 非対称ティー分布) を売上高の対数に当てはめたところ, 非対称ティー分布の当てはまりが良いことがわかった. この結果から, 売上高を従業員数と総資産によって説明するためにコブ・ダグラス型生産関数 (Cobb and Douglas (1928)) を応用し, 両辺の対数をとった, いわゆる両対数モデ
ル (double-log model) を利用して統計モデリング (statistical modeling) を行ったところ, 誤差項に非対称ティー分布を仮定したものが回帰診断の結果として最も良いことが分かった. さらに, 赤池情報量規準 (Akaike Information

Criterion: AIC) を用いてモデル選択 (model selection) を行い, 交差確認法によって評価 (model evaluation) を行った結果も, 上記の結果が肯定されるものとなった.
本研究のデータ (ファイル) を処理する工程は Unix のシェルスクリプトと R スクリプトで一元管理し, 探索的
データ解析の過程も意味のある結果を文書化する工程は, Sweave を利用して LATEX ファイルに R コードを埋め込む
形式で動的に生成した. さらに, これらの全工程を Makefile にスクリプトを記述し, Unix の make コマンドにより
1) https://www.bvdinfo.com/en-gb/

自動実行することによって, 再現可能研究 (reproducible research) を行うことを試みた.

参考文献[1] Akaike, H. (1973) Information theory and an extension of the maximum likelihood principle, Proceedingsof the 2nd International Symposium on Information Theory, Petrov, B. N., and Caski, F. (eds.), Akadimiai
Kiado, Budapest: pp. 267–281.
[2] Azzalini, A. (1985) A class of distributions which includes the normal ones, Scandinavian Journal of
Statistics, Vol. 12, No. 2, pp. 171–178.
[3] Azzalini, A. with the collaboration of A. Capitanio (2014) The Skew-Normal and Related Families, Cambridge University Press, Institute of Mathematical Statistics Monographs.
[4] Cobb, C. W. and P. H. Douglas (1928) A theory of production, American Economic Review, Vol. 18, pp.
139–165.
[5] Efron, B. and T. Hastie (2016) Computer Age Statistical Inference: Algorithms, Evidence, and Data Science,
Cambridge University Press.
[6] Gandrud, C. (2015) Reproducible Research with R and RStudio, Second Edition, CRC Press.
[7] James, G., D. Witten, T. Hastie, and R. Tibshirani (2013) An Introduction to Statistical Learning with
Applications in R, Springer.
[8] 地道正行 (2017-a)『R による対数非対称正規線形モデルによる財務データの統計モデリング』, 商学論究, 第 64
巻, 第 5 号, pp. 159–185, 2017 年 3 月, 関西学院大学商学研究会.
[9] 地道正行 (2017-b)『R を利用した対数非対称分布族にもとづく財務データの統計モデリング』, 経済学論究, 第71 巻, 第 2 号, pp. 141–174, 2017 年 9 月, 関西学院大学経済学部研究会.
[10] 地道正行 (2018-a)『探索的財務ビッグデータ解析 –前処理, データラングリング, 再現可能性–』, 商学論究, 第66 巻, 第 1 号, pp. 1–32, 2018 年 9 月, 関西学院大学商学研究会.
[11] 地道正行 (2018-b)『探索的財務ビッグデータ解析 –データ可視化, 統計モデリング, モデル選択, モデル評価, 動的文書生成, 再現可能研究–』, 商学論究, 第 66 巻, 第 2 号, pp. 1–41, 2018 年 12 月, 関西学院大学商学研究会.
[12] Jimichi, M., Miyamoto, D., Saka, C. and Nagata, S. (2018) Visualization and Statistical Modeling of
Financial Big Data: Log-Linear Modeling with Skew Error, Japanese Journal of Statistics and Data Science,Vol. 1, No. 2, pp. 347–371, https://doi.org/10.1007/s42081-018-0019-1[13] 地道正行, 豊原法彦 (2018) 『景気先行指数の動的文書生成にもとづく再現可能研究』, 豊原法彦編著『関西経済の構造分析』, 第 5 章, pp. 77-111, 中央経済社.
[14] Konishi, S. and G. Kitagawa (2008) Information Criteria and Statistical Modeling, Springer.
[15] Leisch, F. (2002) Sweave: Dynamic generation of statistical reports using literate data analysis, In Wolfgang
Härdle and Bernd Rönz, editors, Compstat 2002 - Proceedings in Computational Statistics, pp. 575-580.
Physica Verlag, Heidelberg. ISBN 3-7908-1517-9.
[16] Ryza, S., U. Laserson, S. Owen, and J. Wills (2016) Advanced Analytics with Spark, O’Reilly.
(玉川 竜司訳 (2016)『Spark による実践データ解析』, オライリー・ジャパン.)[17] Saka, C. and M. Jimichi (2017) Evidence of inequality from accounting data visualisation, Taiwan Accounting Review, Vol. 13, No. 2, pp. 193–234.
[18] 下田倫大, 師岡一成, 今井雄太, 石川 有, 田中裕一, 小宮篤史, 加嵜長門 (2016)『詳解 Apache Spark』, 技術評
論社.
[19] Tukey, J. W. (1977) Exploratory Data Analysis, Addison-Wesley Publishing Co.
[20] Unwin, A. (2015). Graphical Data Analysis with R, Chapman and Hall/CRC.
[21] Xie, Y. (2015) Dynamic Documents with R and knitr, Second Edition, CRC Press.

謝辞
本研究の一部は以下の研究費より助成を得ていることに感謝の意を述べたい:• 科学研究費 基盤研究 C:「グラフィカル・データ・アナリシスによる格差研究と社会環境会計による解決方法の提案」(2016 年 ∼2018 年), 課題番号: 16K04022, 研究代表者: 阪 智香
• 平成 30 年度 学際大規模情報基盤共同利用・共同研究拠点 (JHPCN) 課題: 「財務ビッグデータの可視化と統計モデリング」, 課題番号: jh181001-NWJ, 研究代表者: 地道 正行• 関西学院大学 図書館 図書費 B, 個人研究費
また, BvD 社の増田 歩氏にはデータの抽出に関して多大なるご協力いただいた. ここに感謝の意を述べる.

セミパラメトリック混合効果モデルの適用による繰り返し観測データの解析
大瀧 慈、大谷敬子（広島大学、(株) deCult

１．はじめに必ずしも直線的でない経時的トレンドが想定される個体差を伴う繰り返し観測データに対して、汎用性に富み扱い易く解析結果の解釈が容易と思われるセミパラメトリック混合効果モデルおよびその適用による解析法について、実臨床データを対象にした R の関数による分析を例示しながら紹介する。
[脳梗塞患者の回復状況に関する ADC データ]
脳梗塞の患者さんの「病変」と「正常部」の ADC という細胞密度を表す（MRI により計測される）数値の比（rADC)に関心を持っています。(表１) これまでの経験的知見として、rADC は、脳梗塞急性期（発症 1-2 週）は 1.0 未満となり、亜急性期から慢性期（発症 3-6 週）にかけて 1.0 以上に上昇する傾向を把握しています。(図１) これは、細胞が虚血により膨化して、あたかも密度が上昇したような状態となり、その後細胞が死んで脱落することにより、密度が低下することに起因します。ここで、脳梗塞の原因別（心臓が原因の心原性、血管が原因の動脈硬化性など）で、この密度の変化の速度に差があるか否か、検討してみたい。
（これは、かつて広島大学大学病院脳神経外科からの実際に受けた検討依頼状です。
）
表１．脳梗塞の患者さんの ADC データ (出典:広島大学大学病院脳神経外科)

２．モデル
ADC データについて、患者 (i ) の病型 (typei ) を 2 個のダミー変数（ z1i : typei = 2 の場合
1、その他の場合 0、と z2i : type i = 3 の場合１、その他の場合 0）で表現し、性別 (sex i ) と発症時年齢 ( agei ) を、それぞれ、 =si 1.5 − sex i=ai(男性で 0.5、女性で-0.5) および、

( agei − 70 ) /10 と変換し、その第 j 番目の計測値の測定日 (dij ) を時間変数とし、rADCの常用対数値 ( yij ) を目的変数とし、下記のような繰り返し測定値に関するパラレルプロフ
ァイルモデルを想定した。
3'yij (1, z1i , z2i ) Θ0 + ∑ ( dij − cl ) ⋅ H ( dij − cl ) Θl  (1, si , ai , si ⋅ ai ) + ηi + ε ij ,
=
l =1



(1)

(0)
 β λλλλγ (0) δ (0) λ (0)  (1)（固定効果）未知母数行列であり、
ただし、Θλλλλλγ (1) δ (1) λ (1)  , λ 0,1, 2, は、
=
 β=
(2)
(2)
(2)
(2)
 β λλλλ
γ
δ
λ 


( c0 , c1 , c2 ) = ( 0, 7,14 ) である。 また、 H (t ) = 1 ( t ≥ 0 ) , = 0 ( t < 0 ) は Heaviside 関数であ2==mi , ηi  N (0,τ 2 ), i 1,..., n, 互いに独立な測定誤差および個り、 ε ij  N (0, σ ), j 1,...,

人差を表すランダム変動項である。
モデル（１）は、病型（群）を規定する長さｇの 2 値ダミー変数と、性別や年齢などの共変数が、それぞれ、 z i = (1, z i1 ,..., z g −1 ) ' および xi = (1, x1( i ) ,..., xu( i ) ) 'で表され、目的変数のトレンドが折れ線で表現できることを想定すると、下記のように一般化できる。

 (hi' 1 ⊗ I u )xi' 1 ⊗ z i 



yi 
=
 vec ( Θ0 , Θ1 ,..., Θ k ) + ηi 1mi + ε i ,'' (him ⊗ I u )xim ⊗ z i 
i
i



(2)

ただし、 h ij =(1, ( dij − c1 ) H ( dij − c1 ),..., ( dij − ck ) H ( dij − ck )) ', j =1,..., mi , i =1,..., n.

モデル(2)は、目的変数の経時的トレンドが折れ線で記述されているため、モデル中で使用されている母数と時間変数や共変数との関連性の把握が容易である。また、各モデル間に自然な階層関係を導入することができるので、変数選択による最適モデル探索を含んだ
データ解析が R などの既存のソフトウェアパッケージを利用して容易に実施できる。
なお、
本モデルは、通常の線形重回帰モデルの自然な拡張となっている。

３．データ解析と ADC データの解析結果
既述の ADC データに対して、モデル(1)を用いて R の関数 lmer による解析を行った。その方法と結果については、研究会の当日に発表する。