一般化線型モデルと乱数

ちょっとだけ勉強したのでメモしておこう。

一般化線型モデル

まず、確率変数ベクトル $\mathbf{Y}$ の期待値ベクトルを
\[
\boldsymbol\theta=\mathrm{E}(\mathbf Y)
\] と書いておきます。
一般化線型モデルと呼ばれている統計モデルでは、
\[
g(\boldsymbol\theta)=\mathbf{X}\boldsymbol{\beta}
\] と呼ばれる期待値の構造を仮定していたと思います。
ここで、リンク関数 $g$ には、逆関数がちゃんと存在する常識的なものを用います。
理論的な話をきちんと知りたい場合は、位置/尺度パラメータを持つ指数型分布族における正準リンク関数 (canonical link function) について勉強する必要があります (McCullough and Nelder 1989 など)。
逆関数があれば、
\[
\boldsymbol{\theta}=g^{-1}(\mathbf{X}\boldsymbol{\beta})
\] こう書けるわけですね。
そして、普通は確率変数ベクトル $\mathbf{Y}$ が、この期待値ベクトルを位置パラメータとして持つ指数型分布族
\[
p(\mathbf{Y} | \boldsymbol{\theta}, \boldsymbol{\phi}),
\] に従うと仮定したものが扱われているはずです。
たぶん。
ただし、$\boldsymbol{\phi}$ は期待値ベクトル以外のパラメータとします (分散などの尺度パラメータとか)。
各 $Y_i$ については独立で、$\mathbf{x}_i$ に伴って、期待値の構造が変わるのでベクトル表記にしています。

モデルの性質とか、指数型分布族とか、色々勉強しないとなぁ。

乱数の生成

必要なもの

$\mathbf{Y}$ が従う分布族
期待値ベクトルと必要に応じてその他のパラメータ
逆リンク関数

logistic model (ベルヌーイ分布、ロジットリンク)

data logistic;
  beta1 = 1;
  do r = 1 to 100;
  do i = 1 to 100;
    do x1 = 0 to 1;
      logit = beta1 * x1;
      theta = 1 / (1 + exp(-logit));
      y = rand('Bernoulli', theta);
      output;
    end;
  end; end;
run;
ods exclude all;
proc genmod data = logistic;
  by r;
  class x1 / param = ref;
  model y = x1 / noscale d = b link = logit;
  ods output ParameterEstimates = sim_x1;
run;
ods select all;
proc means data = sim_x1;
  where Parameter = "x1";
  var Estimate;
run;

ベルヌーイ分布 y = rand('Bernoulli', theta);
期待値ベクトル logit = beta1 * x1;
- その他のパラメータ無し
逆ロジット関数 theta = 1 / (1 + exp(-logit));

シミュレーションを想定しているので、r は繰り返しの回数です。
x1 が適当な線形関数、i が各 x1 毎のデータ数です。
二項分布で考えて r/n 型にも書けます。

normal linear model (正規分布、恒等リンク)

data linear;
  beta1 = 1;
  do r = 1 to 100;
  do i = 1 to 100;
    do x1 = 0 to 1;
      id = beta1 * x1;
      theta = id;
      phi = 1;
      y = rand('Normal', theta, phi);
      output;
    end;
  end; end;
run;
ods exclude all;
proc genmod data = linear;
  by r;
  class x1 / param = ref ref = first;
  model y = x1 / d = n link = identity noint;
  ods output ParameterEstimates = sim_x1;
run;
ods select all;
proc means data = sim_x1;
  where Parameter = "x1";
  var Estimate;
run;

正規分布 rand('Normal', theta, phi);
期待値ベクトル id = beta1 * x1;
- その他のパラメータ分散 phi
恒等リンク theta = id;