補完のためのロバストな比率モデルの一般化

2025年度統計関連学会連合大会（38. ロバストネス） 補完のためのロバストな 比率モデルの一般化 総務省統計研究研修所 和田 かず美

目次 1. 概要 2. 研究の背景 3. 比率モデルについて 4. ロバスト推定法 5. 比率モデルのロバスト化と一般化 6. 実装関数について

1. 概要 • 比率モデルの𝑥と相関のある不等分散誤差𝜖𝑖 を、回帰モデルのような等分散な 誤差𝜀𝑖 で再定式化することにより、M-推定のIRLSアルゴリズムによりロバスト推 定を可能とした これらの誤差項の関係性 𝜖𝑖 = 𝜀𝑖 𝑥𝑖0.5 を拡張し、右辺の𝑥𝑖 の冪乗部分 0.5 を可変(γ) として一般化した。 • γが所与の場合について、比率βをロバスト推定する関数をRで実装し、パッ ケージ化して公開した • 一般化比率モデルについて、 比率βとγを同時推定する関数も実装したが、精 度があまり良くないので、より良い推定方法を探したい • 比率モデル： 𝑦𝑖 = 𝛽𝑥𝑖 + 𝜖𝑖 , 𝜖𝑖 ~𝑁(0, 𝑥𝜎 2 ) 一般化比率モデル： 𝛾 𝑦𝑖 = 𝛽𝑥𝑖 + 𝜀𝑖 𝑥𝑖 , 𝜀𝑖 ~𝑁(0, 𝜎 2 ) 3

2. 研究の背景 • 調査統計において、欠測と同様に外れ値も避け難いことが多く、ロバ スト推定により外れ値の影響を緩和するることにより推定精度が向上 することができる • 実データは、等分散ではないことも多く、線形回帰モデルを適用する 場合に変数変換が必要になるが、変数変換する場合、（対数軸上で はなく実軸上での）平均値の推定が不安定になり、またデータ分散 が大きい場合に上ブレする • 欠測補完（imputation）において、比率モデルは回帰モデル同様に よく使用されるが、その不等分散な誤差項のために、回帰モデルのよ うなM-推定によるロバスト化が難しい 4

変数変換の問題点～補完値への影響 ➔ 実データの平均値の値の振れが大きくなる。 ※ ただし実データの平均値を安定させるためには、変数変 換しない場合でも推定時に乗率を考慮する必要がある。 ➔ 変数変換を施すと、補完値算出時に残差分散による補正が 必要になる。これを忘れると過小推計になるが、データ分散が 大きいと補正すれば過大になる。 ※ ロバスト推定する場合は、残差分散自体も補正されるた めに過大推計にはならない。 5

このとき補完値の推計式は? ln( yi ) = a + b  ln( xi ) +  i yˆ = exp( aˆ + bˆ  ln( x )) i 指数化（対数化の逆） i 現行の推計式 不偏推定量は n  2   ri   1 i =1   ˆ yˆ i = exp( aˆ + b  ln( xi ))  exp   2 (n − p)   残差r による補正 ※ 自然対数変換の場合   6

補完値の推計式の求め方（１） －変換のない単回帰の場合－ モデル式 yi = a + b  xi +  i  i : 誤差項 期待値をとる E[ yi ] = E[a + b  xi +  i ] 確率変数は誤差項だけ = a + b  xi + E[ i ] 回帰の誤差項は期待値（平均）0になるのでこの項は消える。 aとbに回帰パラメータの推計値を用いてyの推計値を計算する。 推計式 yˆi = aˆ + bˆ  xi 7

補完値の推計式の求め方（２） －対数変換した単回帰の場合－ モデル式 ln( yi ) = a + b  ln( xi ) +  i  i : 回帰誤差 ri : 回帰残差 期待値をとる E[ln( yi )] = E[a + b  ln( xi ) +  i ] 指数化（対数化の逆） μ : 誤差の平均 σ : 誤差の分散 確率変数は誤差項だけ b E[ yi ] = exp( a)  xi  E[exp(  i )] exp(ei)は対数正規分布なのでこの項は消えない。 1 2 b = exp( a)  xi  exp(  +  ) 2 誤差の平均・分散は未知なので、実測 される残差rの平均(0)と分散で推計。 推計式 1 bˆ yˆi = exp( aˆ )  xi  exp(  2 i=1 ri ) n 2 n−2 補正項がなければ過小推計になる。一方、残差が正規分布よりも裾が 長い場合はこの補正により推定値が過大推計になる可能性もある。 ロバスト回帰で推定すると この問題を解消できる 8

3. 比率モデルについて

単回帰モデル y 𝑦 = 𝛼 + 𝛽𝑥 + 𝜀 𝑉(𝜀) = 𝜎 2 誤差項𝜀の分布 𝑟𝑖 𝜎 𝑦ො𝑖 = 𝛼ො + 𝛽𝑥𝑖 𝜎 0 𝑦𝑖 𝑥𝑖 x 10

比率モデル y 誤差項𝜖の分布 𝑦 = 𝛽𝑥 + 𝜖, 𝑉(𝜖) = 𝑥𝜎 2 xの値が増えるにつれて誤 差の分散も大きくなる 𝑟𝑖 𝑥𝑖 𝜎 𝑥𝑗 𝜎 መ 𝑖 𝑦ො𝑖 = 𝛽𝑥 0 𝑥𝑖 𝑥𝑗 𝑦𝑖 x 11

比率モデルによる補完値の推定 𝑦𝑖 を目的変数（補完対象の調査項目）、𝑥𝑖 を説明変数（補完対象項 目と相関が高い調査項目）、𝛽 は𝑦𝑖 と𝑥𝑖 との比率とすると、比率補完の モデルは以下のとおり。 𝑦𝑖 = 𝛽𝑥𝑖 + 𝜖𝑖 𝜖𝑖 は誤差項で、 𝜖𝑖 ~𝑁(0, 𝜎 2 𝑥𝑖 )。← 誤差項𝜖𝑖 は𝑥𝑖 と相関のある不等分 散な誤差。欠測があれば真の𝛽は未知なので、下式により推定する。 σ𝑘∈obs 𝑦𝑖 𝛽መ = obs:欠測なく𝑥𝑖 と𝑦𝑖 がセットで得られている観測値 σ𝑘∈obs 𝑥𝑖 De Waal et al. (2011) Handbook on Statistical Data Editing and Imputation, Wiley handbooks in survey methodology, John Wiley & Sons, p. 244-245. 12

比率モデルと回帰モデルの違い 比率モデル 𝑦𝑖 = 𝛽𝑥𝑖 + 𝜖𝑖 不等分散誤差 𝜖𝑖 ~𝑁(0, 𝜎 2 𝑥𝑖 ) 回帰モデル 𝑦𝑖 = 𝛼 + 𝛽𝑥𝑖 + 𝜀𝑖 等分散誤差 𝜀𝑖 ~𝑁(0, 𝜎 2 ) boundary of 2σ boundary of 2σ 13

不等分散なデータへの回帰補完と比率補完の違い 1. 回帰補完の場合、変数変換が必要 2. 比率補完の場合、通常変換は不要 調査データは等分散ではない場合が多い 回帰の場合、平均値の推定があまり安定しない e.g. 対数変換の例 log 𝑦𝑖 = 𝛼 + 𝛽 ∙ log 𝑥𝑖 + 𝜀𝑖 1 σ𝑛𝑖=1 𝜀𝑖2 𝑦ො𝑖 = 𝑒𝑥𝑝 𝛼 + 𝛽 ∙ log 𝑥𝑖 ∙ 𝑒𝑥𝑝 ∙ 2 𝑛−𝑝 回帰モデルに変換が必要な場合、補完には比率モデル の方が良い 14

4. ロバスト推定法

Statistical Data Editing Vol.2の掲載事例 米センサス局 月次卸売調査の在庫データ（当期と前期） OLS （通常の回帰推定） IRLSによるロバスト回帰の例 16

出典: UNSE/UNECE Statistical Data Editingシリーズ 公的統計編纂のためのベストプラクティス集 Vol.1 Methods and Techniques (1994) https://unece.org/info/publications/pub/21885 Vol.2 Methods and Techniques (1997) https://unece.org/info/publications/pub/21884 Vol.3 Impact on Quality (2006) https://unece.org/info/publications/pub/21881 Bienias et al. (1997) 紹介の抵抗回帰 (Vol.2 収録) ロバストだが実用化が容易な繰返し計算を用いたM-推定量を算出する IRLS(IWLS)アルゴリズム 古典的だが収束も速く実用的 ウエイト関数にTukeyのbiweight、残差の尺度にAAD(平均絶対偏差)を使用 17

• https://unece.org/info/publications/pub/21885
• https://unece.org/info/publications/pub/21884
• https://unece.org/info/publications/pub/21881

IRLSの効果(回帰推定) デモ用データ データ数 N=102 うち100データ(黒点)は、相関0.7、分散1の 正規乱数。 残り2データ(赤点)は人工外れ値。 試算結果 OLS(普通の最小二乗法) : 赤線 切片 4.74 傾き 0.95 外れ値を除いたOLS : 橙線 切片 30.0 傾き 0.7 IRLS: 青線 切片 29.10 傾き 0.71 18

IRLSの仕組み(1) 1. 通常の最小二乗法を行う。 2. 標準化残差𝑒𝑖 とその平均絶対偏差s から、TukeyあるいはHuberの関数 に従ってウエイト𝑤𝑖 を算出する。 Huberのウエイト関数 Tukeyのbiweight関数 𝑒𝑖 2 if 𝑒𝑖 ≤ (𝑐 × 𝑠) 𝑤𝑖 = 1 − 𝑐×𝑠 if 𝑒𝑖 > (𝑐 × 𝑠) 𝑤𝑖 = 0 2 if 𝑒𝑖 ≤ (𝑘 × 𝑠) 𝑤𝑖 = 1 𝑘×𝑠 ൞ if 𝑒𝑖 > (𝑘 × 𝑠) 𝑤𝑖 = |𝑟𝑖 | データの基準幅 データの基準幅 i はデータ番号 Tukeyのc 及びHuberのk は推計のロバスト性を調節するための定数（tuning constant） 3. ウエイト付きで最小二乗法を行い、同様にウエイト付きで計算した平均絶 対偏差sの1回前の値との変化率が0.01以下になるまで2～3を繰返す。 19

IRLSの仕組み(2) 20

5. 比率モデルのロバスト化と一般化

ロバスト化への障害 通常の比率補完のモデルの場合、誤差項𝜖𝑖 の分散は、 𝜖𝑖 ~𝑁(0, 𝜎 2 𝑥𝑖 ) で𝑥𝑖 に比例する。 比率モデルによる推定をロバスト化するためには、回帰モデルと同様に、 誤差項を 𝑥𝑖 と関係を持たない 𝜀𝑖 ~𝑁 0, 𝜎 2 で定式化したい。 𝜀𝑖 = 𝜖𝑖 / 𝑥𝑖 なので、モデル式を 𝑥𝑖 で割ると 𝑦𝑖 𝑥𝑖 = 𝛽 𝑥𝑖 + 𝜀𝑖 𝑦𝑖 = 𝛽𝑥𝑖 + 𝜀𝑖 𝑥𝑖 Cochran, W. G. (1977) Sampling Techniques, 3rd ed., John Wiley & Sons. 22

比率モデルの一般化 0.5 誤差項の 𝑥𝑖 のべき乗を𝛾に置き替える モデル式 推定式 𝑦𝑖 (1−𝛾) + 𝜀𝑖 𝛾 = 𝛽𝑥𝑖 𝑥𝑖 1−2𝛾 σ 𝑦𝑖 𝑥𝑖 𝛽መ = 2(1−𝛾) σ 𝑥𝑖 𝑥𝑖 と関係を持たない誤差項 ※ 𝛾は任意の定数 23

γ=1のとき: 大昔の比率モデル 𝑦𝑖 = 𝛽 + 𝜀𝑖 , 𝑥𝑖 𝑦𝑖 𝜀𝑖 = − 𝛽~ 𝑁 0, 𝜎 2 𝑥𝑖 1 𝑦𝑖 መ 𝑦𝑖 = 𝛽𝑥𝑖 + 𝜀𝑖 𝑥𝑖 , 𝛽= ෍ 𝑛 𝑥𝑖 A’ γ =1/2のとき: 通常の比率モデル 𝑦𝑖 𝑦𝑖 = 𝛽 𝑥𝑖 + 𝜀𝑖 , 𝜀𝑖 = − 𝛽 𝑥𝑖 ~ 𝑁(0, 𝜎 2 ) 𝑥𝑖 𝑥𝑖 σ 𝑦𝑖 መ 𝑦𝑖 = 𝛽𝑥𝑖 + 𝜀𝑖 𝑥𝑖 , 𝛽 = B’ σ 𝑥𝑖 γ =0のとき: 切片のない単回帰モデル 𝑦𝑖 = 𝛽𝑥𝑖 + 𝜀𝑖 , 𝜀𝑖 = 𝑦𝑖 − 𝛽𝑥𝑖 ~ 𝑁(0, 𝜎 2 ) σ 𝑦𝑖 𝑥𝑖 𝛽መ = σ 𝑥𝑖2 C’ 24

一般化比率モデルのロバスト化 1−2𝛽 σ 𝑦𝑖 𝑥𝑖 𝑟Ƹ = 2(1−𝛽) σ 𝑥𝑖 σ 𝑤𝑖 𝑦𝑖 𝑤𝑖 𝑥𝑖 1−2𝛽 𝑟Ƹ = σ 𝑤𝑖 𝑥𝑖 2 1−𝛽 γ=1のとき: 1 𝑤𝑖 𝑦𝑖 𝒓ො 𝑟𝑜𝑏𝐴 = ෍ 𝑛 𝑤𝑖 𝑥𝑖 A σ 𝑤𝑖 𝑦𝑖 𝒓ො 𝑟𝑜𝑏𝐵 = σ 𝑤𝑖 𝑥𝑖 B γ =1/2のとき: 25

ウェイト関数について 赤：c=4 青：c=6 黒：c=8    1 −  e   w( e) =    c    0  2 2 赤：k=1.15 青：k=1.72 黒：k=2.30 𝑒 𝑒 | e | c | e | c ある程度中心部から遠い観測値の影 響を完全排除できる 1  w( e) =  k  |e| | e | k | e | k 中心部から非常に遠い観測値でも、 その影響を完全には排除しない 26

計算アルゴリズム 回帰推定に用いる繰返し加重最小二乗法（IRLS: Iterative Reweighted Least Squares）の仕組みを、回帰の最小二乗 法ではなく比率の推定に適用 ✓ 計算が簡単で非常に収束が速い ✓ 推定パラメータが一つなので、ウェイト関数にTukeyの biweightを使ってもループしない 27

n=1000 正常値の相関は0.6 人為的に3割の外れ値を加えても、正常値のみによる比推計に近い結果を得ることができる 28

n=15 ・データ量が少ない場合でも、推定量AとBは外れ値の影響を受けにくい ・推定量A’よりもB’の方が規模の大きい外れ値の影響を受けやすい 29 29

n=15 かなり極端な外れ値を添加しても、推定量AとBは影響を受けにくい 30

6. Rによる実装

robRatioパッケージ Github [https://github.com/kazwd2008/robRatio]公開済み CRANへは現在投稿中 収録関数一覧 モデル 比率 比率 比率 比率 回帰 回帰 回帰 回帰 比率 親関数 子関数 robRatio robRatio robRatio robRatio robReg robReg robReg robReg robGR RrT.aad RrT.mad RrH.aad RrH.mad Tirls.aad Tirls.mad Hirls.aad Hirls.mad - ウェイト関数 残差の尺度 Tukey’s biweight 平均絶対偏差 (AAD) Tukey’s biweight 中央絶対偏差 (MAD) Huber weight AAD Huber weight MAD Tukey’s biweight AAD Tukey’s biweight MAD Huber weight AAD Huber weight MAD Tukey’s biweight AAD γの値 所与 所与 所与 所与 － － － － 同時推定 32

• https://github.com/kazwd2008/robRatio

親関数と子関数の関係について • 子関数は、回帰モデル・比率モデルそれぞれについて、ウェイト関数2 種類、尺度基準二種類の組み合わせで関数4つずつ作成 • ウェイト関数と尺度の組み合わせで、ロバスト度をユーザーがコント ロールする調整定数（tuning constant, 引数c1）の値が異なる • 親関数は、直接調整定数を指定する引数c1ではなく、４, 6または8の 値のみ指定すれば良い引数tpを用意し、子関数による調整定数の違 いをユーザーが意識せずに済むようにした • robRatio及びその子関数については、γの値は0、0.5または1でユー ザーが指定する • 現在実装準備中のrobGRは、βとγを同時推定するもの 33

ウェイト関数と尺度の組み合わせと調整定数 ウェイト関数 尺度基準 tp=4 Tukey Tukey Tukey Huber Huber Huber SD AAD MAD SD AAD MAD c1=5.01 c1=4 c1=3.38 c1=1.44 c1=1.15 c1=0.97 tp=6 tp=8 c1=7.52 c1=6 c1=5.07 c1=2.16 c1=1.72 c1=1.46 c1=10.03 c1=8 c1=6.76 c1=2.88 c1=2.30 c1=1.94 ※ SDは標準偏差（Standard Deviation） 34

5.4 rogGR関数 二段階最小二乗法 (Weighted two-stage least squares) 1. 𝛽෠ = σ 𝑦𝑖 𝑥𝑖1−2𝛾 ൗσ 𝑥𝑖 2 1−𝛾 を、適切な初期値 𝛾 0について推定 ෠ 2. 求めた𝛽の値を使って、γを推定する（非ロバスト推定） 1 ⊤ −1 ⊤ i. ステップ1で推定した𝛽෠ を用いて、 𝑋 𝑋 𝑋 𝑌 により 𝛾 を推定。 ここで、 𝑋 = ⋮ 1 最小二乗法 ෠ 1 log 𝑦1 − 𝛽𝑥 log 𝑥1 ⋮ , 𝑌= ⋮ ෠ 𝑛 log(𝑥𝑛 ) log 𝑦𝑛 − 𝛽𝑥 ii. 推定された 𝛾ො を用いて新たな 𝛽෠ を推定 3. 収束条件 1−2𝛾 ෡ = σ 𝑤𝑖 𝑦𝑖 𝑤𝑖𝑥𝑖 ステップ2で得られた 𝛾ො を用いて、 𝜷 によるロバストな推定を標準化残差の尺度𝜎の変化が十分小さくなるまで繰 ො σ 𝑤𝑖 𝑥𝑖 2 1−𝛾 ෠ 𝑖 )ൗ𝑥 𝛾 を使用し、ウェイト関数に基づい り返す。ここで、ウェイト 𝑤𝑖 は、直近の繰り返し計算での準残差（quasi residuals） 𝑟෕𝑖 = (𝑦𝑖 − 𝛽𝑥 𝑖 て算出する 加重最小二乗法 4. 𝛾ො を 𝑋 ⊤𝑊𝑋 −1𝑋 ⊤𝑊𝑌 によりロバスト推定する。ここで、 𝑤 = 𝑤1 0 ⋱ 0 5. 以下によるロバストな𝜷 と𝜸の推定を収束条件を満たすまで繰り返す (𝑘−1) i. 直近の 𝛾 (𝑘−1) と 𝑤𝑖 ii. 新たなウェイト 𝑤𝑖𝑘 を用いて 𝛾 𝑘 を推定 を用いて、 𝛽෠ 𝑘 = 𝑤は対角成分以外は全てゼロ σ 𝑤𝑖 𝑦𝑖 𝑤𝑖 𝑥𝑖 1−2𝛾 を推定 σ 𝑤𝑖 𝑥𝑖 2 1−𝛾 𝑤𝑛 使用するウェイト関数はTukeyのbiweight, 残差の尺度はAAD 35

二段階最小二乗法 𝛾 2 𝜎𝑥𝑖 からγを推定する Var 𝜖𝑖 = 1) 誤差𝜖𝑖 は観測できないが、代わりに残差𝑟𝑖 を使用し、𝑟𝑖2 = により求める 2 ∙ log 𝑟𝑖 = 2 ∙ log 𝜎 + 2𝛾 ∙ log 𝑥𝑖 y = a + ⊤ −1 γ x ⊤ 通常の最小二乗法 𝛾= 𝑋 𝑋 𝑋 𝑦 加重最小二乗法 𝛾 = 𝑋 ⊤ 𝑊𝑋 −1 𝑋 ⊤ 𝑊𝑦 2 መ 𝑦𝑖 − 𝛽𝑥𝑖 γが回帰推定できる 1 𝑋= ⋮ 1 log 𝑥1 ⋮ , log(𝑥𝑛 ) 𝑤1 𝑤= 0 𝑤は対角成分以外は全てゼロ ⋱ 0 መ 1 log 𝑦1 − 𝛽𝑥 𝑌= ⋮ መ 𝑛 log 𝑦𝑛 − 𝛽𝑥 𝑤𝑛 36

シミュレーションによる推定結果(1) 初期コード β=2, γ=0, 0.5の場合の乱数データに よるシミュレーション結果1000回分 (パッケージのtestthatフォルダ収録) あまり良い推定にならない。 γの推定に偏りがみられる。 また、γの推定が荒れるために3～4割 収束しない。 Data A （真値β=2, γ=0） βとγの平均 2.000258 -4.396863 βとγの中央値 1.999786 -3.526448 Data B（真値β=2, γ=0.5） βとγの平均 1.998948 -4.283506 βとγの中央値 1.996778 -3.400593 37

シミュレーションによる推定結果(2) 改良版 【コードの修正点】 繰り返し計算の中でγの推定が入る と尺度基準が膨れることが多いため、 尺度が大きくなるγの推定は捨てる ルーチンを挿入した。 【結果の概要】 γのマイナス値が大幅に減ったが、恐 らく2SLSの推定精度が悪く尺度基準 𝜎ො を縮小できないために、γが更新さ れず初期値のまま収束してしまう（こ の例では、γの初期値は0） 。 38

robGR 新旧コードの比較 初期値 データ データA γ=0 平均値 β 中央値 γ β 標準偏差 γ β γ 旧コード（全推計） 1.9996 0.0805 2.0000 0.3071 0.0225 3.1806 β=2 旧コード（収束分のみ） 2.0002 -1.1557 2.0000 -0.9200 0.0082 2.7240 γ=0.5 新コード 2.0000 3.0497 1.9998 1.4150 0.0076 3.8958 旧コード（全推計） 1.9962 0.4418 1.9986 0.6836 0.1707 3.0382 β=2 旧コード（収束分のみ） 2.0034 -0.9421 1.9995 0.1936 0.0636 2.6198 γ=1 新コード 1.9994 3.5428 1.9988 2.3395 0.0584 4.1149 旧コード（全推計） 1.9990 0.0584 2.0000 0.2444 0.0237 3.1741 β=2 旧コード（収束分のみ） 1.9998 -1.1249 2.0000 -0.8920 0.0086 2.7448 γ=0.5 新コード 1.9999 5.6053 1.9997 5.0000 0.0078 2.2088 旧コード（全推計） 1.9947 0.4484 2.0000 0.7609 0.1782 3.0556 β=2 旧コード（収束分のみ） 2.0030 -0.8714 1.9999 5.0000 0.0575 2.6682 γ=1 新コード 1.9989 5.8992 1.9985 5.0000 0.0608 2.8610 データB データA γ=5 使用コードと収束 データB - 新・旧コードいずれもγの推定値が真値から外れる - 改良コードは全てのデータについて計算が収束し、βの推定精度は旧コードよりもわずかに良い。ただし、γの 推定は改善したとはいえない - とりあえず新コードをパッケージに収録した 39

今後の予定 • CRAN掲載 • できればrogGRの推定を改善するアルゴリズムを探す 40

関連資料 IRLSアルゴリズムによるM-推定（回帰） • 和田(2012)多変量外れ値の検出～繰返し加重最小二乗（IRLS）法による欠測値の補完方法～、 統計研究彙報第６９号、pp.23-52. https://www.stat.go.jp/training/2kenkyu/ihou/69/pdf/2-2-692.pdf • 和田・野呂(2019)ロバスト回帰推定へのウェイト関数や残差尺度の影響について、統計研究彙報 第76号、pp.101-114. https://www.stat.go.jp/training/2kenkyu/ihou/69/pdf/2-2-692.pdf ロバストな比率の推定 • 2016年度統計関連学会連合大会発表スライド「平成28年経済センサス‐活動調査のためのロバ ストな比率補完の方法について」 https://www.nstac.go.jp/sys/files/static/services/society_paper/28_03_02_PowerPoint.pdf • 経済センサス-活動調査に関する研究会「経済センサス‒活動調査のためのロバストな比率補完 の方法について」 https://www.nstac.go.jp/sys/files/static/services/society_paper/28_03_02.pdf • Wada, Sakashita and Tsubaki (2021). Robust Estimation for a Generalised Ratio Model. Austrian Journal of Statistics, 50(1), 74–87. https://doi.org/10.17713/ajs.v50i1.994 • Wada and Takata (2019) An Algorithm of Generalized Robust Ratio Model Estimation for Imputation, In JSM proceedings, pp.3120-3128. https://ww2.amstat.org/meetings/proceedings/2019/data/assets/pdf/1199680.pdf 41

• https://www.stat.go.jp/training/2kenkyu/ihou/69/pdf/2-2-692.pdf
• https://www.nstac.go.jp/sys/files/static/services/society_paper/28_03_02_PowerPoint.pdf
• https://www.nstac.go.jp/sys/files/static/services/society_paper/28_03_02.pdf
• https://doi.org/10.17713/ajs.v50i1.994
• https://ww2.amstat.org/meetings/proceedings/2019/data/assets/pdf/1199680.pdf