スパース経時測定データに対する 関数データ解析 ～方法論とRによる実装方法の紹介～ 滋賀大学データサイエンス学部 松井 秀俊

2 • 観測個体それぞれが経時的に測定された データを分析する方法を紹介 • 観測個体それぞれがまばらな時点でのみ 観測値を得たデータに対して 主成分分析を行う方法と これを実行するRパッケージを説明 （こちらの内容の続編） • 観測個体それぞれを 低次元のベクトル値に変換することで まばらに観測されたデータの経時変化の 特徴を定量化できる 本資料の内容 主成分得点plot （横：第１，縦：第２） 子どもの身長のデータ （点が観測値） 一人のデータを 同じ色で表している

3 • はじめに • 関数データ解析 • スパース経時測定データ • スパース経時測定データのモデル化 • 平均関数の推定 • 共分散関数の推定 • スパース経時測定データに対する関数主成分分析 • 適用例とRによる実装 目次

4 • 右図のように，観測個体（子供）それぞれが 時間等の経過に伴い繰り返し計測された 形式のデータ（経時測定データ）を考える • このようなデータに対して各個体を時間の関数 として表し，観測データ（ベクトルデータ）の 代わりに関数をデータとして扱い分析する 方法は関数データ解析とよばれている • 関数データ解析を適用することで，観測時点が 不均一であったり，観測時点・時点数が 個体ごとに異なっていたとしても 比較的容易に分析ができる 関数データ解析 (Functional Data Analysis; FDA) 子供の身長の推移データ Ramsay & Silverman (2005) Kokoszka & Reimherr (2017)

5 • 関数データ解析では，観測個体それぞれに対して 基底関数展開などを用いて独立して関数化する ことで関数データを得る方法が多い （関数（曲線）1本1本を独立に推定） • この方法は各個体において観測時点数が十分あり， かつ観測時点が定義域を十分カバーしている という前提であれば適切に関数化できる • 下図のデータ（上図のデータを，各個体で 時点数が10になるよう間引いたもの）のように， 個体によっては観測時点や時点数が まばらな場合もある （例：患者の診察データ） スパース経時測定データ 子供の身長データ（上）と 観測点を間引いたデータ（下）

6 • 下図のようなデータはスパース(経時測定)データ (sparse data)とよばれる • それと対比して，上図のようなデータは デンスデータ (dense data)とよばれる • スパースデータに対しては，各個体に対して 独立して関数化を行おうとしても うまく関数を推定できない場合が多い • 全ての観測データの情報を共有（プール）して 分析することで，観測個体それぞれの特徴を 捉えたい スパース経時測定データ 子供の身長データ（上）と 観測点を間引いたデータ（下）

7 • スパース経時測定データに対する分析手法の1つである PACE (Principal Analysis by Conditional Estimation) について紹介 (Yao et al., 2005a; Kokoszka & Reimherr, 2017) ✓スパース経時測定データから主成分得点を推定 ✓主成分得点をデータとして扱うことで 一般的な多変量解析と同様の方法を適用可能 • PACEを実装したRパッケージ”fdapace”の使い方を紹介 • どういう形式のデータを入力しどういう出力が得られるか ※「方法論は置いといてとりあえず実装方法を知りたい」という方は スライド33へ 本資料の目的

8 • はじめに • 関数データ解析 • スパース経時測定データ • スパース経時測定データのモデル化 • 平均関数の推定 • 共分散関数の推定 • スパース経時測定データに対する関数主成分分析 • 適用例とRによる実装 目次

9 • スパースデータに対しては，観測時点数がまばらなことから デンスデータとは異なり個々の関数データを直接得ることは難しい • スパースデータに対する分析アプローチとして，平均関数や共分散関数を 推定する方法をはじめに紹介 • これらの推定には，基底関数展開または局所多項式回帰が用いられる • スパースデータのモデリング方法には，主に次の2種類がある （どちらかというと前者の方が最近は多く用いられているように感じる（私見）） • 共分散関数平滑化に基づく方法 ←PACEはこちら 推定された平均関数と共分散関数から，関数主成分得点を推定する • 混合効果モデルに基づく方法 (James et al., 2000ほか) スパースデータの発生メカニズムを混合効果モデルと捉え 平均関数（固定効果）と個体変動（変量効果）を推定する問題に帰着させる スパースデータに対するアプローチ

10 • 第 𝑖 番目の観測における 𝛼 番目の時点 𝑡𝑖𝛼 ∈ 𝒯 (𝒯 ⊂ ℝ：定義域) における観測値を 𝑥𝑖𝛼 とおく • データ𝑥𝑖𝛼 は次のように，「平均関数」と「個体変動」と「ノイズ」の ３つの要素により生成されていると仮定 𝑥𝑖𝛼 = 𝑥𝑖 𝑡𝑖𝛼 = 𝜇 𝑡𝑖𝛼 + 𝜀𝑖 𝑡𝑖𝛼 + 𝛿𝑖𝛼 𝑖 = 1, … , 𝑛：個体番号, 𝛼 = 1, … , 𝑁𝑖 ：時点番号，𝑖によって異なる 𝜇 𝑡 ：平均関数, 𝛿𝑖𝛼 ：誤差 (𝑖, 𝛼両方について独立) 𝜀𝑖 𝑡 ：第𝑖個体の変動，異なる時点間で相関を仮定 • 全ての時点 𝑡𝑖𝛼 ; 𝑖 = 1, … , 𝑛, 𝛼 = 1, , … , 𝑁𝑖 は，（各𝑖では不十分でも） 定義域 𝒯上に十分に分布しているものとする スパース経時測定データの定式化

11 • 基本的な考え方は，デンスデータに対する関数化と同じ • ただし各個体ではなく，全てのデータを用いて平滑化を行う • 平均関数𝜇 𝑡 は次のように基底関数展開で表されると仮定 𝜇 𝑡 = ෍ 𝑚=1 𝑀 𝑤𝑚 𝜙𝑚 𝑡 = 𝒘𝑇𝝓 𝑡 𝒘 = 𝑤1 , … , 𝑤𝑀 𝑇：係数, 𝝓 𝑡 = 𝜙1 𝑡 , … , 𝜙𝑀 𝑡 𝑇：基底関数 • 次の（全観測 𝑥11 , … , 𝑥𝑛𝑁𝑛 に基づく）２乗誤差最小化により 係数ベクトル𝒘を推定 𝐿 𝒘 = ෍ 𝑖=1 𝑛 ෍ 𝛼=1 𝑁𝑖 𝑥𝑖𝛼 − 𝒘𝑇𝝓 𝑡𝑖𝛼 2 平均関数の推定法１：基底関数展開（１） (Green & Silverman, 1993)

12 • 全ての観測値からなるσ 𝑖 𝑛 𝑁𝑖 次元ベクトルとσ 𝑖 𝑛 𝑁𝑖 × 𝑀基底関数行列 𝒙 = 𝑥11 , … , 𝑥𝑛𝑁𝑛 𝑇 , Φ = 𝝓 𝑡11 , … , 𝝓 𝑡𝑛𝑁𝑛 𝑇 を用いることで，２乗誤差𝐿 𝒘 は次で表される 𝐿 𝒖 = 𝒙 − Φ𝒘 𝑇 𝒙 − Φ𝒘 • これより，𝒘の最小２乗推定量として次を得る（σ 𝑖 𝑛 𝑁𝑖 > 𝑀のとき） ෝ 𝒘 = Φ𝑇Φ −1 Φ𝑇𝒙 • 最小２乗法の代わりに𝐿2 正則化法を用いる場合は，推定量は次で与えられる ෝ 𝒘 = Φ𝑇Φ + 𝜆Ω −1 Φ𝑇𝒙 𝜆：正則化パラメータ，Ω：非負値定符号行列 平均関数の推定法１：基底関数展開（２）

13 • 基底関数展開と同様，すべての観測値 𝑥11 , … , 𝑥𝑛𝑁𝑛 に対して平滑化を行う • 次の基準を最小化する𝜷 = 𝛽0 , … , 𝛽𝑝 𝑇を推定する (Loader, 1999) 𝐿 𝜷 = ෍ 𝑖=1 𝑛 ෍ 𝛼=1 𝑁𝑖 𝐾 𝑡 − 𝑡𝑖𝛼 ℎ 𝑥𝑖𝛼 − 𝛽0 ෍ 𝑗=1 𝑚 𝛽𝑗 𝑡 − 𝑡𝑖𝛼 𝑗 2 𝐾 ⋅ ：カーネル関数 • 任意の 𝑡 について，𝐿 𝜷 最小化により得られる𝛽0 の推定値 መ 𝛽0 を用いて Ƹ 𝜇 𝑡 = መ 𝛽0 𝑡 とする • 次数は 𝑗 = 1 （局所線形回帰）が用いられることが多い 平均関数の推定法２：局所多項式回帰（１）

14 • 次のようにベクトル・行列を定義 𝒙 = 𝑥11 , … , 𝑥𝑛𝑁𝑛 𝑇 , 𝑍 = 𝑡11 − 𝑡 0 ⋯ 𝑡11 − 𝑡 𝑚 ⋮ ⋱ ⋮ 𝑡𝑛𝑁𝑛 − 𝑡 0 ⋯ 𝑡𝑛𝑁𝑛 − 𝑡 𝑚 𝐾ℎ = diag 𝐾 𝑡 − 𝑡11 ℎ , … , 𝐾 𝑡 − 𝑡𝑛𝑁𝑛 ℎ • このとき， 𝐿 𝜷 は次で表される 𝐿 𝜷 = 𝒙 − 𝑍𝜷 𝑇𝐾ℎ 𝒙 − 𝑍𝜷 • これを𝜷について解くことで，次の推定値を得る ෡ 𝜷 = 𝑍𝑇𝐾ℎ 𝑍 −1 𝑍𝑇𝐾ℎ 𝒙 • バンド幅ℎはCVなどを用いて決定する 平均関数の推定法２：局所多項式回帰（２）

15 • 基底関数展開・局所平滑化それぞれを使って スパースな成長データの平均曲線を推定（調整パラメータは適当） • いずれの方法でも平均的な振る舞いを捉えられている 平均関数の推定・適用例 基底関数展開による平均曲線(赤) 局所平滑化による平均曲線(赤)

16 • 共分散関数の推定は，基本的には平均関数の推定である １次元データ（曲線）の平滑化を２次元（曲面）に置き換えたもの • しかし，単なる曲面平滑化とは次の２点が異なる • 推定された関数に正定値性と対称性が必要 • ノイズ𝛿𝑖𝛼 のために，曲面の対角成分に特別な対処が必要 • 個体変動𝜀𝑖 𝑡 の共分散関数を𝑐 𝑠, 𝑡 ，ノイズ𝛿𝑖𝛼 の分散を𝜎2 𝑡𝑖𝛼 とおく (𝜎2 𝑡 は連続関数とする) • このとき，ランダム関数 𝑋 ⋅ の共分散関数は次で与えられる Cov 𝑋 𝑠 , 𝑋 𝑡 = 𝑐 𝑠, 𝑡 + 𝜎2 𝑡 𝐼 𝑡 = 𝑠 共分散関数の推定

17 Cov 𝑋 𝑠 , 𝑋 𝑡 = 𝑐 𝑠, 𝑡 + 𝜎2 𝑡 𝐼 𝑡 = 𝑠 • 共分散関数の対角線上で𝜎2 𝑡 だけ「かさ上げ」 され不連続に（右図） • デンスデータであればこのかさ上げは 無視できるか平滑化により除去できるが スパースデータであれば無視できなくなる • そこで， 𝑐 𝑠, 𝑡 と𝜎2 𝑡 を分けて推定する 共分散関数の推定：成分の分解 Wang et al. (2016)より

18 • 曲面 𝑐 𝑠, 𝑡 について，（不連続な）対角線上の点の集合𝑐 𝑡, 𝑡 を除いた点で 曲面（2次元）平滑化を行う • まず，෤ 𝑥𝑖𝛼1 ෤ 𝑥𝑖𝛼2 = 𝑥𝑖𝛼1 − Ƹ 𝜇 𝑡𝑖𝛼1 𝑥𝑖𝛼2 − Ƹ 𝜇 𝑡𝑖𝛼2 とおき， ෥ 𝝌 = ෤ 𝑥𝑖𝛼1 ෤ 𝑥𝑖𝛼2 |𝑖 = 1, … , 𝑛, 𝛼1 = 1, … , 𝑁𝑖 , 𝛼2 = 1, … , 𝑁𝑖 , 𝛼1 ≠ 𝛼2 を目的変数， ෤ 𝒕 = 𝑡𝑖𝛼1 , 𝑡𝑖𝛼2 |𝑖 = 1, … , 𝑛, 𝛼1 = 1, … , 𝑁𝑖 , 𝛼2 = 1, … , 𝑁𝑖 , 𝛼1 ≠ 𝛼2 を2変数説明変数 とした回帰モデルを構築し，その推定により得られる曲面を ҧ 𝑐 𝑠, 𝑡 とおく • ҧ 𝑐 𝑠, 𝑡 を得る方法としては，平均関数と同様に基底関数展開または 局所平滑化による曲面平滑化が用いられる（次頁・次々頁） 共分散関数の推定：曲面平滑化

19 • 説明変数の集合 ෤ 𝒕 = ǁ 𝑠, ǁ 𝑡 の ǁ 𝑠, ǁ 𝑡 ，目的変数の集合 ෥ 𝝌 をここではベクトル （𝑁′ = σ 𝑖=1 𝑛 𝑁𝑖 𝑁𝑖 − 1 次元）とみなし，その要素番号を 𝑖′ とおく • 2変数基底関数展開を用いて，共分散関数𝑐 𝑠, 𝑡 𝑠 ≠ 𝑡 は次で表されると仮定 （基底関数としてはthin-plate splineや動径基底関数などが用いられる） 𝑐 𝑠, 𝑡 = ෍ 𝑚=1 𝑀 𝑤𝑚 𝜙𝑚 𝑠, 𝑡 = 𝒘𝑇𝝓 𝑠, 𝑡 • 1変数基底関数を用いて次のように表してもよい（加法モデル） 𝑐 𝑠, 𝑡 = 𝒘𝑠 𝑇𝝓 𝑠 + 𝒘𝑡 𝑇𝝓 𝑡 • 次の２乗誤差最小化により係数ベクトル𝒘を推定することで， ҧ 𝑐 𝑠, 𝑡 を得る 𝐿 𝒘 = ෍ 𝑖′=1 𝑁′ ෤ 𝜒𝑖′ − 𝒘𝑇𝝓 ǁ 𝑠𝑖′ , ǁ 𝑡𝑖′ 2 曲面平滑化の方法１：基底関数展開 𝒘 = 𝑤1 , … , 𝑤𝑀 𝑇：係数 𝝓 𝑠, 𝑡 = 𝜙1 𝑠, 𝑡 , … , 𝜙𝑀 𝑠, 𝑡 𝑇：基底関数

20 • 前頁と同様，ベクトルとみなした説明変数 ǁ 𝑠, ǁ 𝑡 ，目的変数 ෥ 𝝌 を扱う • 次の基準を最小化する𝜷 = 𝛽0 , … , 𝛽𝑝 𝑇を推定する 𝐿 𝜷 = ෍ 𝑖′=1 𝑁′ 𝐾 𝑠 − ǁ 𝑠𝑖′ ℎ , 𝑡 − ǁ 𝑡𝑖′ ℎ 𝑥𝑖𝛼 − 𝛽0 − ෍ 𝑗=1 𝑚 𝛽𝑠𝑗 𝑠 − ǁ 𝑠𝑖′ 𝑗 − ෍ 𝑗=1 𝑚 𝛽𝑡𝑗 𝑡 − ǁ 𝑡𝑖′ 𝑗 2 𝐾 ⋅,⋅ ：カーネル関数 • 任意の 𝑠, 𝑡 について，𝐿 𝜷 最小化により得られる𝛽0 の推定値 መ 𝛽0 を用いて ҧ 𝑐 𝑠, 𝑡 = መ 𝛽0 𝑠, 𝑡 とする 曲面平滑化の方法２：局所平滑化

21 • 平滑化により得られる曲面は一般的に対称でないうえ 正定値性も保証されない • そこで，推定曲面を, ׬ 𝑐𝑋 𝑠, 𝑡 𝑢 𝑡 𝑑𝑡 = 𝜆𝑢 𝑠 に基づき 固有関数により展開し，固有値が負のパートを削る つまり， ҧ 𝑐 𝑠, 𝑡 の固有値（昇順）を𝜆1 ≥ 𝜆2 ≥ ⋯としたとき， 初めて負になる固有値の番号を𝐾として ҧ 𝑐 𝑠, 𝑡 ≈ ෍ 𝑘=1 𝐾−1 𝜆𝑘 𝑢𝑘 𝑠 𝑢𝑘 𝑡 = Ƹ 𝑐 𝑠, 𝑡 のように近似する（𝑢𝑘 𝑡 ：固有値 𝜆𝑘 に対応する固有関数） • これにより対称性と正定値性の２つを満たす • 線形代数でいうところの，任意の行列を正の固有値のみを使って 低ランク近似するようなもの 共分散関数の推定：固有関数展開

22 • 続いて，分散関数 Var 𝑋 𝑡 (= Cov 𝑋 𝑡 , 𝑋 𝑡 )を ෥ 𝒙 と ෤ 𝒕 の対角線上の曲線（1次元）平滑化により推定 つまり， ෥ 𝒙 = ෤ 𝑥𝑖𝛼 2 |𝑖 = 1, … , 𝑛, 𝛼 = 1, … , 𝑁𝑖 , 𝛼 = 𝛼1 = 𝛼2 を目的変数， ෤ 𝒕 = 𝑡𝑖𝛼 |𝑖 = 1, … , 𝑛, 𝛼 = 1, … , 𝑁𝑖 , 𝛼 = 𝛼1 = 𝛼2 を1変数説明変数 として回帰分析を行う．得られた推定値を ǁ 𝑐 𝑡, 𝑡 とおく • これを用いて， 𝜎2 𝑡 （かさ上げ分）を ො 𝜎2 𝑡 = ǁ 𝑐 𝑡, 𝑡 − Ƹ 𝑐 𝑡, 𝑡 で推定 つまり， Cov 𝑋 𝑠 , 𝑋 𝑡 の推定値は ǁ 𝑐 𝑡, 𝑡 = Ƹ 𝑐 𝑡, 𝑡 + ො 𝜎2 𝑡 • あるいは， ෤ 𝑥𝑖𝛼 2 − Ƹ 𝑐 𝑡𝑖𝛼 , 𝑡𝑖𝛼 を平滑化することでも ො 𝜎2 𝑡 を推定できる 共分散関数の推定：対角成分の推定

23 • はじめに • 関数データ解析 • スパース経時測定データ • スパース経時測定データのモデル化 • 平均関数の推定 • 共分散関数の推定 • スパース経時測定データに対する関数主成分分析 • 適用例とRによる実装 目次

24 • 前回のスライドで紹介した関数主成分分析により，ランダム関数𝑋 𝑡 を 主成分得点𝑍1 , 𝑍2 , … と固有関数𝑢1 𝑡 , 𝑢2 𝑡 , …を用いて次のように表した 𝑋 𝑡 = 𝜇 𝑡 + ෍ 𝑘=1 ∞ 𝑍𝑘 𝑢𝑘 𝑡 • スパースデータでは，観測データ 𝑥𝑖𝛼 ; 𝛼 = 1, … , 𝑁𝑖 に代えて主成分得点 𝑧𝑖1 , 𝑧𝑖2 , …を扱うことで，次のようなメリットがある • 不均一なデータを均一データ（普通の多変量データ）として扱える • 主成分分析を用いているので、主成分得点は一般的に低次元になる • 全定義域 𝒯 上での関数データ𝑥𝑖 𝑡 を復元できる（けど注意が必要） 関数主成分得点の利用

25 • 主成分得点を求めるには，内積 𝑧𝑖𝑘 = 𝑥𝑖 − 𝜇, 𝑢𝑘 = ׬ 𝑥𝑖 𝑡 − 𝜇 𝑡 𝑢𝑘 𝑡 𝑑𝑡 𝑖 = 1, … , 𝑛 の分散の最大化問題を解けばよい • ・・・と言いたいところだが，スパースデータでは（デンスデータの場合と は異なり）直接関数データ𝑥𝑖 𝑡 が得られないため，この最適化問題を 解くことができない • そこで，前述の方法で推定された共分散関数 Ƹ 𝑐 𝑠, 𝑡 から求められる 固有値 መ 𝜆1 , መ 𝜆2 , …，固有関数ො 𝑢1 𝑡 , ො 𝑢2 𝑡 , …から主成分得点𝑧𝑖1 , 𝑧𝑖2 , … を推定する • より具体的には，条件付き期待値 E 𝑍1 𝑋 𝑡1 , … , 𝑋 𝑡𝑁 を 第１主成分得点𝑍1 の推定量とする (Yao et al., 2005a) 関数主成分得点の推定

26 • 簡単のために𝜇 𝑡 = 0とし，また𝑍1 , 𝑋 𝑡1 , … , 𝑋 𝑡𝑁 は全て正規分布に従うと 仮定 (𝑡1 , … , 𝑡𝑁 ：観測時点) • このとき， 𝑍1 と𝑋 𝑡𝛼 の共分散は Cov 𝑍1 , 𝑋 𝑡𝛼 = E න𝑋 𝑡 𝑢1 𝑡 𝑑𝑡 𝑋 𝑡𝛼 = න𝑐(𝑡, 𝑡𝛼 )𝑢1 𝑡 𝑑𝑡 = 𝜆1 𝑢1 𝑡𝛼 で与えられる • これより， 𝑍1 , 𝑋 𝑡1 , … , 𝑋 𝑡𝑀 𝑇の分散共分散行列 Ξ は次で与えられる Ξ = 𝜆1 𝜆1 𝑢1 𝑡1 ⋯ 𝜆1 𝑢1 𝑡𝑀 𝜆1 𝑢1 𝑡1 𝑐 𝑡1 , 𝑡1 + 𝜎2 𝑡1 ⋯ 𝑐 𝑡1 , 𝑡𝑁 ⋮ ⋮ ⋱ ⋮ 𝜆1 𝑢1 𝑡𝑀 𝑐 𝑡𝑁 , 𝑡1 ⋯ 𝑐 𝑡𝑁 , 𝑡𝑁 + 𝜎2 𝑡𝑁 関数主成分得点の計算手順（1/3）

27 • 下記の命題を利用すると，𝑍1 |𝑋 𝑡1 , … , 𝑋 𝑡𝑁 は期待値・分散共分散行列が それぞれ次の正規分布に従う 期待値： 𝜆1 𝑢1 𝑡1 ⋮ 𝜆1 𝑢1 𝑡𝑁 𝑇 𝑐 𝑡1 , 𝑡1 + 𝜎2 𝑡1 ⋯ 𝑐 𝑡1 , 𝑡𝑁 ⋮ ⋱ ⋮ 𝑐 𝑡𝑁 , 𝑡1 ⋯ 𝑐 𝑡𝑁 , 𝑡𝑁 + 𝜎2 𝑡𝑁 −1 𝑋 𝑡1 ⋮ 𝑋 𝑡𝑁 分散共分散行列： 𝜆1 − 𝜆1 𝑢1 𝑡1 ⋮ 𝜆1 𝑢1 𝑡𝑁 𝑇 𝑐 𝑡1 , 𝑡1 + 𝜎2 𝑡1 ⋯ 𝑐 𝑡1 , 𝑡𝑁 ⋮ ⋱ ⋮ 𝑐 𝑡𝑁 , 𝑡1 ⋯ 𝑐 𝑡𝑁 , 𝑡𝑁 + 𝜎2 𝑡𝑁 −1 𝜆1 𝑢1 𝑡1 ⋮ 𝜆1 𝑢1 𝑡𝑁 関数主成分得点の計算手順（2/3） 命題 2つの確率変数ベクトル𝒀1 , 𝒀2 を並べたベクトルが次の正規分布に従うとする 𝒀1 𝒀2 ~𝑁 𝝁1 𝝁2 , Σ11 Σ12 Σ21 Σ22 このとき， 𝒀1 |𝒀2 の条件付き期待値，分散共分散行列はそれぞれ次で与えられる E 𝒀1 |𝒀2 = 𝜇1 + Σ12 Σ22 −1 𝒀2 − 𝝁2 , Var 𝒀1 |𝒀2 = Σ11 − Σ12 Σ22 −1Σ21

28 • この条件付き期待値を，第 𝑖 観測の第𝑘主成分得点の推定値 Ƹ 𝑧𝑖𝑘 = E 𝑧𝑖𝑘 |𝑥𝑖 𝑡𝑖1 , … , 𝑥𝑖 𝑡𝑁𝑖 とする • さらに，固有関数の推定値ො 𝑢1 𝑡 , ො 𝑢2 𝑡 , …を 共分散関数の推定値 Ƹ 𝑐 𝑠, 𝑡 の固有関数展開で求めることで，関数データの予測値 𝑥𝑖 𝑡 = ෍ 𝑘=1 𝐾 Ƹ 𝑧𝑖𝑘 ො 𝑢𝑘 𝑡 を得る（𝐾 は主成分の個数） これは全定義域 𝒯 での関数データを復元したものとみなせる • 得られる（推定される）順序をまとめると次の通り： • デンスデータ：𝑥𝑖𝑡 → 𝑥 𝑡 → ො 𝑢𝑘 𝑡 → Ƹ 𝑧𝑘 • スパースデータ： 𝑥𝑖𝑡 → Ƹ 𝑐 𝑠, 𝑡 → ො 𝑢𝑘 𝑡 → Ƹ 𝑧𝑘 → 𝑥 𝑡 関数主成分得点の計算手順（3/3）

29 • 前述の方法で推定された関数主成分得点を用いて、スパースデータに対する 関数回帰モデルの推定ができる (Yao et al., 2005b) • 関数回帰モデルにはさまざまな種類のものがあるが、どの種類に対しても 大体適用可能 • Scalar-on-function モデル（説明変数が関数データ・目的変数がスカラーデータ） • Function-on-scalar モデル（説明変数がスカラーデータ・目的変数が関数データ） • Function-on-function モデル（説明変数・目的変数ともに関数データ） • 上記の3モデルそれぞれについて，関数データがスパースデータである場合の 推定方法を簡単に紹介 （簡単のため，説明変数，目的変数共に中心化され切片=0とする） 関数回帰分析

30 • 関数説明変数𝑥𝑖 𝑡 と，スカラー目的変数𝑦𝑖 との関係を次でモデル化 𝑦𝑖 = න𝛽 𝑡 𝑥𝑖 𝑡 𝑑𝑡 + 𝜀𝑖 𝑖 = 1, … , 𝑛, 𝛽 𝑡 : 係数, 𝜀𝑖 : 誤差 • 𝑥𝑖 𝑡 に対応するスパースデータに対して，前述の方法により 主成分得点 ො 𝒛𝑖 = Ƹ 𝑧𝑖1 , … , Ƹ 𝑧𝑖𝐾 𝑇と固有関数 ෝ 𝒖 𝑡 = 𝑢1 𝑡 , … , 𝑢𝐾 𝑡 𝑇を推定 • 係数関数𝛽 𝑡 は，固有関数の線形結合により表されると仮定 𝛽 𝑡 = ෍ 𝑘=1 𝐾 𝛽𝑘 ො 𝑢𝑘 𝑡 = 𝜷𝑇ෝ 𝒖 𝑡 𝜷 = 𝛽1 , … , 𝛽𝐾 𝑇 • このとき， ො 𝑢1 𝑡 , … , ො 𝑢𝐾 𝑡 の正規直交性より関数回帰モデルは次で表される 𝑦𝑖 = ො 𝒛𝑖 𝑇𝜷 + 𝜀𝑖 • 通常の線形回帰モデルと同様，最小２乗法などを用いて𝜷を推定すれば 係数関数は መ 𝛽 𝑡 = ෡ 𝜷𝑇ෝ 𝒖 𝑡 と推定できる 関数回帰モデル：Scalar-on-function

31 • スカラー説明変数𝑥𝑖 と，関数目的変数𝑦𝑖 𝑡 との関係を次でモデル化 𝑦𝑖 𝑡 = 𝑥𝑖 𝛽 𝑡 + 𝜀𝑖 𝑡 𝑖 = 1, … , 𝑛, 𝛽 𝑡 : 係数, 𝜀𝑖 𝑡 : 誤差 • 前頁と同様， 𝑦𝑖 𝑡 に対する関数主成分得点ො 𝒛𝑖 と固有関数ෝ 𝒖 𝑡 を推定し， さらに回帰係数が𝛽 𝑡 = 𝜷𝑇ෝ 𝒖 𝑡 と表されると仮定 • このとき，関数回帰モデルは次で表される ො 𝒛𝑖 𝑇ෝ 𝒖 𝑡 = 𝑥𝑖 𝜷𝑇ෝ 𝒖 𝑡 + 𝜀𝑖 𝑡 • 積分２乗誤差 න ෍ 𝑖=1 𝑛 ො 𝒛𝑖 𝑇ෝ 𝒖 𝑡 − 𝑥𝑖 𝜷𝑇ෝ 𝒖 𝑡 2 𝑑𝑡 の最小化を考えることで，𝛽 𝑡 の推定量として次を得る መ 𝛽 𝑡 = ෡ 𝜷𝑇ෝ 𝒖 𝑡 , ෡ 𝜷 = 1 σ 𝑖 𝑥𝑖 2 Σ𝑖 𝑥𝑖 Ƹ 𝑧𝑖1 , … , Σ𝑖 𝑥𝑖 Ƹ 𝑧𝑖𝐾 𝑇 関数回帰モデル：Function-on-scalar

32 • 関数説明変数𝑥𝑖 𝑠 と，関数目的変数𝑦𝑖 𝑡 との関係を次でモデル化 （前頁のモデルと比べて，異なる時点𝑠, 𝑡での関連も見られる） 𝑦𝑖 𝑡 = න𝛽 𝑠, 𝑡 𝑥𝑖 𝑠 𝑑𝑠 + 𝜀𝑖 𝑡 𝑖 = 1, … , 𝑛, 𝛽 𝑠, 𝑡 : 係数, 𝜀𝑖 𝑡 : 誤差 • 𝑥𝑖 𝑠 , 𝑦𝑖 𝑡 それぞれに対して関数主成分得点ො 𝒛𝑖 , ෝ 𝒘𝑖 と固有関数ෝ 𝒖 𝑡 , ෝ 𝒗 𝑡 を 推定し， 𝑥𝑖 𝑠 = ො 𝒛𝑖 𝑇ෝ 𝒖 𝑠 ，𝑦𝑖 𝑡 = ෝ 𝒘𝑖 𝑇ෝ 𝒗 𝑡 で近似 • 回帰係数は𝛽 𝑠, 𝑡 = ෝ 𝒖 𝑠 𝑇𝐵ෝ 𝒗 𝑡 と表されると仮定 (𝐵 = 𝑏𝑘𝑙 𝑘𝑙 ) • このとき，関数回帰モデルは次で表される ෝ 𝒘𝑖 𝑇ෝ 𝒗 𝑡 = ො 𝒛𝑖 𝑇𝐵ෝ 𝒗 𝑡 + 𝜀𝑖 𝑡 • 積分２乗誤差の最小化により，𝛽 𝑠, 𝑡 の推定量として次を得る መ 𝛽 𝑠, 𝑡 = ෝ 𝒖 𝑠 𝑇 ෠ 𝐵ෝ 𝒗 𝑡 , ෠ 𝑏𝑘𝑙 = σ𝑖 Ƹ 𝑧𝑖𝑘 ෝ 𝑤𝑖𝑙 σ 𝑖 Ƹ 𝑧𝑖𝑘 2 関数回帰モデル：Function-on-function

33 • はじめに • 関数データ解析 • スパース経時測定データ • スパース経時測定データのモデル化 • 平均関数の推定 • 共分散関数の推定 • スパース経時測定データに対する関数主成分分析 • 適用例とRによる実装 目次

34 • 39人の子供の身長の推移データ（下図左）に対して観測値をいくつか間引き 疑似的にスパースデータを生成（下図右） • このスパースデータに対してPACEを適用し 主成分得点の推定と元のデータの復元を試みる 適用例 子供の身長の推移データ 左図から観測値を間引いたもの ※ 視認性のため 10人のみ図示

35 # ライブラリ読み込み library(fda) # 成長データ利用のため library(fdapace) # FPCA実行のため ## 成長データ dat = t(growth$hgtm) # 人数×時点数の行列 n = ncol(growth$hgtm) # サンプルサイズ ## スパースデータ化 ## 各個体で観測時点をランダムに10時点に間引く dat_sp = Sparsify(dat, growth$age, 10) ## データplot par(mfrow=c(1,2), mar=c(4,4,1,1)) plot(dat_sp$Lt[[1]], dat_sp$Ly[[1]], type=“o”, pch=16, lty=2, xlim=range(dat_sp$Lt), ylim=range(dat_sp$Ly)) for(i in 1:n){ lines(dat_sp$Lt[[i]], dat_sp$Ly[[i]], col=i, type="o", pch=16, lty=2) } ## 参考：間引く前のオリジナルデータ matplot(growth$age, growth$hgtm, col=1:10, type="o", pch=16) Rの実行コード（1/2） dat_sp（スパースデータ） Lt：時点，Ly：子供の身長のリスト 元々行列型だったデータが Sparsify関数により リスト型に変換される （中略） dat（デンスデータ） 行：子供，列：時点の行列

36 # スパース関数主成分分析適用（カーネル関数にガウスカーネル使用） # 仕様の詳細はFPCA関数のヘルプを参照 res = FPCA(dat_sp$Ly, dat_sp$Lt, list(kernel='gauss')) # FPCA結果plot plot(res) # 主成分得点plot plot(res$xiEst[,1], res$xiEst[,2], xlab="FPC1", ylab="FPC2") abline(h=0, v=0, lty=2) # 予測曲線plot CreatePathPlot(res, showMean = T, showObs = F) # 推定共分散関数plot CreateCovPlot(res) Rの実行コード（2/2） それぞれのplotの出力結果について次頁以降で紹介

37 • 左上： 共分散関数を曲面平滑化により 推定するにあたり 時点のペア 𝑡𝑖𝛼1 , 𝑡𝑖𝛼2 の 両方で観測値がある 𝑖 の数 • この例では，4歳,5歳の両方の 時点で観測値をもつ子供は 4人以上（赤）だが， 5歳,6歳では2人のみ（青） • ここではスパースデータとして 扱っているので 対角成分は除いている （スライド17参照） 結果（1/5）: FPCA結果plot

38 • 右上：推定平均関数 Ƹ 𝜇 𝑡 • 左下：各主成分の寄与率 を示したスクリープロット このデータに対しては 第１関数主成分が90%以上の 情報を保持していることが わかる • 右下：第１～第３関数主成分 に対応する固有関数 黒：ො 𝑢1 𝑡 , 赤：ො 𝑢2 𝑡 , 緑：ො 𝑢3 𝑡 結果（2/5）: FPCA結果plot

39 • 推定された第1，第2主成分得点 Ƹ 𝑧𝑖1 , Ƹ 𝑧𝑖2 𝑖 = 1, … , 𝑛 の 散布図を描画（右上） • 第1主成分が最大となる子供（赤）と 第2主成分が最大となる子供（緑）の身長の推移を 右下の元のデータ（間引き前のもの）に対応 • 前頁にあった固有関数の特徴と照らし合わせると 第1主成分は全期間を通した身長の高さ， 第2主成分は成長期の到来の遅さを表していると 解釈できる 結果（3/5）: 主成分得点plot 主成分得点（上）と元データ（下） 固有関数（再掲） ※右の図と比べて対応する色が違ってます。ご了承ください。

40 • 下図はスパースデータ（左）から推定された主成分得点と固有関数を用いて 予測された関数データ𝑥𝑖 𝑡 （中央，黒太線は平均関数） • 間引き前のデンスデータ（右）と比べるとある程度の傾向は捉えられているが 細部での予測はうまくいってない • 各個体の予測曲線を利用する際は慎重に 結果（4/5）: 予測曲線plot スパースデータからの予測曲線 スパースデータ デンスデータ（参考）

41 • 下図左はスパースデータに対する推定共分散関数 Ƹ 𝑐 𝑠, 𝑡 （𝑠 = 𝑡 パート除く） 下図右は元のデンスデータに対して推定された標本共分散関数 • 成長期以降で共分散が大きくなるという傾向を捉えている 結果（5/5）: 推定共分散関数plot スパースデータからの推定共分散関数 デンスデータからの推定共分散関数

42 • 観測個体1つ1つがまばらな時点での観測値を得たデータ（スパースデータ） に対して関数主成分分析を行う方法であるPACEを紹介 • 各個体で関数化を行うことは困難なため，全データを利用して平均関数と 共分散関数を推定し，これらを用いて主成分分析を行う • 主成分得点の推定値は，関数データが与えられた下での事後期待値として 求められる • 主成分得点の推定値が得られれば，それをそのまま多変量データとして 扱ったり，関数回帰モデルの推定へ応用できたりする • RではfdapaceパッケージのFPCA関数でPACEを実行可能 まとめ

43 – Green, P. and Silverman, B. W. (1993). Nonparametric Regression and Generalized Linear Models. Chapman and Hall/CRC. – James, G., Hastie, T. and Sugar, C. (2000). Principal component models for sparse functional data. Biometrika 87, 587-602. – Kokoszka, P. and Reimherr, M. (2017). Introduction to Functional Data Analysis. Chapman and Hall/CRC. – Loader, C. R. (1999). Local Regression and Likelihood, Springer. – Ramsay, J. O. and Silverman, B. W. (2005). Functional Data Analysis 2nd ed. Springer. – R package “fda” https://cran.r-project.org/web/packages/fda/index.html – R package “fdapace” https://cran.r-project.org/web/packages/fdapace/index.html – Wang, J., Chiou, J. M. and Müller, H-G. (2016). Review of functional data analysis. Annual Review of Statistics and Its Application 3, 257-295. – Yao, F., Müller, H-G. and Wang, J. (2005a). Functional data analysis for sparse longitudinal data. Journal of the American Statistical Association 100, 577-590. – Yao, F., Müller, H-G. and Wang, J. (2005b). Functional linear regression analysis for longitudinal data. Annals of Statistics 33, 2873-2903. 参考文献