Zenn
🏋️‍♂️

Stata と R の weight(s) オプション

2023/05/13に公開
R
データ分析
stata
ウエイト
tech

本記事の目的

この記事では、データ分析を行う際に考慮するウエイト（weight）について整理します。基本的には回帰分析を想定しており、後半では単回帰分析の例もとに、 Stata と R の結果を比較しています。説明は Dupraz (2013) を参考にしました。色々と省略している箇所もあるので、詳細は必要に応じて元論文をご覧ください。

Stata で選択可能な 4 つの *weight オプション

Stata には 4 つの weight オプションがあります。

  • fweight: frequency weights
  • aweight: analytic weights
  • pweight: probability weights
  • iweight: importance weights

（iweight は特殊なオプションとのことで、この記事ではそれ以外の 3 つを扱います）

各ウエイトの説明

結論、どのオプションでも点推定値には違いがありません。異なるのは推定値の標準誤差です。

fweight

fweight が想定するのは、1 つのデータセットに重複したデータが入っているようなケースです。例えば、以下のケースでは i 行目が重複しています。

\boldsymbol{X}=\begin{bmatrix} 1 & x_{11} & \cdots & x_{1k} \\ \vdots & \vdots & \vdots & \vdots \\ 1 & x_{i1} & \cdots & x_{ik} \\ 1 & x_{i1} & \cdots & x_{ik} \\ \vdots & \vdots & \vdots & \vdots \\ 1 & x_{n1} & \cdots & x_{nk} \end{bmatrix},\ \boldsymbol{y}=\begin{bmatrix} y_1\\ \vdots\\ y_i\\ y_i\\ \vdots\\ y_{n} \end{bmatrix}

fweight で推定した回帰係数の分散共分散行列は以下のとおりです。n はデータセットの行数、k は説明変数の数です。残差分散 \hat{u}^\mathrm{t}\hat{u} を n-(k+1) で割っています

\sum=\frac{\hat{u}^\mathrm{t}\hat{u}}{n-(k+1)}(\boldsymbol{X^\mathrm{t}\boldsymbol{X}})^{-1}

aweight

ここでは以下のモデルを推定することを想定します。

y_{ij}=T_j\tau+x_j\beta+u_{ij}

i は 個人、j はグループです。y_{ij} と u_{ij} は個人単位で異なりますが、説明変数は添え字が j でグループ単位です。またここでは、計算効率およびランダムなのはグループ単位であるという理由から、以下の式を推定することを考えます。

\bar{y}_j=T_j\tau+x_j\beta+\bar{u}_{i}

上述のとおり、このモデルを aweight で推定すると点推定値は fweight と一致します。異なるのは標準誤差で、分散共分散行列は以下の式で表されます。ここでの m はグループやクラスタの数です。\tilde{u} の上についている記号（チルダ）は重みづけされていることを示しています。

\sum=\frac{\tilde{\hat{u}}^\mathrm{t}\tilde{\hat{u}}}{m-(k+1)}(\tilde{\boldsymbol{X}}^\mathrm{t}\tilde{\boldsymbol{X}})^{-1}

pweight

pweight は母集団からの抽出確率を反映するウエイトです。参考にした論文の例では、都市圏と地方の調査対象者がおり、都市圏では 2 人に 1 人が、地方では 10 人に 1 人が抽出された状況が想定されています。この場合、都市圏では 1 人の対象者が 2 人分の、地方では 10 人分の重みもちます。この重みを考慮しない分析を行うと、分析結果に都市圏の意見が実態より大きく反映されることになります。

pweight で実行した結果も回帰係数の点推定値は aweight、fweight と同じですが、ここでもやはり標準誤差の結果が異なります。pweight の標準誤差は以下のとおりです（Sandwich Estimater）。\boldsymbol{W} は m \times m の対角行列で、対角成分には各グループにおける残差の 2 乗（\hat{u}_j^2）が入っています。

\sum=\frac{m}{m-(k+1)}(\tilde{\boldsymbol{X}}^\mathrm{t}\tilde{\boldsymbol{X}})^{-1}(\tilde{\boldsymbol{X}}^\mathrm{t}\boldsymbol{W}\tilde{\boldsymbol{X}})(\tilde{\boldsymbol{X}}^\mathrm{t}\tilde{\boldsymbol{X}})^{-1}

データで再現

aweight と pweight を使って実際に分析してみます。

扱うデータ

PISA 2018（OECD による生徒の学習到達度調査）のデータセットを PISA 2018 Database から取得し、dat_pisa というデータを作成しています。

扱った変数はデータセットのうちの一部で、数学的リテラシー、読解力、生徒のウエイトです。ウエイト（W_FSTUWT）を考慮し、読解力（PV1READ）を独立変数、数学的リテラシー（PV1MATH）を従属変数として回帰分析を行います。

R
> dat_pisa %>% head()
# A tibble: 6 × 3
  PV1MATH PV1READ W_FSTUWT
    <dbl>   <dbl>    <dbl>
1    704.    705.     166.
2    581.    570.     166.
3    692.    648.     166.
4    628.    672.     166.
5    682.    672.     166.
6    784.    770.     166.

Stata

aweightと pweight で回帰係数（Coefficient）の値が一致していることがわかります。一方、標準誤差（Std. err.）の値は pweight の方が大きくなっています。論文では aweight かつ robust standard error（ロバスト標準誤差）を実行することで pweight の結果と一致すると記載がありましたので、そちらも試しています。

aweight

実行コードはこちら。

Stata
reg pv1math pv1read [aweight=w_fstuwt]

結果はこちら。

Stata
. reg pv1math pv1read [aweight=w_fstuwt]
(sum of wgt is 1,078,921.3259125)
------------------------------------------------------------------------------
     pv1math | Coefficient  Std. err.      t    P>|t|     [95% conf. interval]
-------------+----------------------------------------------------------------
     pv1read |   .7102128    .006849   103.70   0.000     .6967865    .7236392
       _cons |   171.1731   3.510538    48.76   0.000     164.2912    178.0549
------------------------------------------------------------------------------

pweight

実行コードはこちら。

Stata
reg pv1math pv1read [pweight=w_fstuwt]

結果はこちら。

Stata
. reg pv1math pv1read [pweight=w_fstuwt]
------------------------------------------------------------------------------
             |               Robust
     pv1math | Coefficient  std. err.      t    P>|t|     [95% conf. interval]
-------------+----------------------------------------------------------------
     pv1read |   .7102128    .007116    99.81   0.000      .696263    .7241627
       _cons |   171.1731   3.714015    46.09   0.000     163.8923    178.4538
------------------------------------------------------------------------------

aweight + robust standard error

実行コードはこちら。

Stata
reg pv1math pv1read [aweight=w_fstuwt], vce(robust)

結果はこちら。

Stata
. reg pv1math pv1read [aweight=w_fstuwt], vce(robust)
(sum of wgt is 1,078,921.3259125)
------------------------------------------------------------------------------
             |               Robust
     pv1math | Coefficient  std. err.      t    P>|t|     [95% conf. interval]
-------------+----------------------------------------------------------------
     pv1read |   .7102128    .007116    99.81   0.000      .696263    .7241627
       _cons |   171.1731   3.714015    46.09   0.000     163.8923    178.4538
------------------------------------------------------------------------------

R

R では lm 関数で weights という引数を設定できます。W_FSTUWT を指定して分析すると結果は以下のとおりです。Stata で aweight を指定したときと同じ結果になっています。

実行コードはこちら。

R
model_pisa <- dat_pisa %>% 
  lm(data = ., 
     formula = PV1MATH ~ PV1READ, weights = W_FSTUWT)

summary(model_pisa)

結果はこちら（一部省略）。

R
Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
PV1READ       0.710213   0.006849  103.70   <2e-16 ***
(Intercept) 171.173051   3.510538   48.76   <2e-16 ***

{lmtest} パッケージを使うとロバスト標準誤差が計算できます。 Replicating Stata's "robust" option in R を見ると、type 引数で HC1 を設定すると Stata の結果と一致するようです。

実行コードはこちら。

R
lmtest::coeftest(x = model_pisa, vcov. = vcovHC, type = "HC1")

結果はこちら。pweight の結果と一致しています。

R
t test of coefficients:

              Estimate Std. Error t value  Pr(>|t|)    
PV1READ       0.710213   0.007116  99.805 < 2.2e-16 ***
(Intercept) 171.173051   3.714015  46.088 < 2.2e-16 ***

今後

今回は従属変数が連続変数のケースを扱いましたが、2 値のケースも書くかもしれません。

参考

Dupraz, Y. (2013). Using weights in Stata. Accessed on August.

GitHubで編集を提案

Discussion