SASで「データ解析のための統計モデリング入門」を読む 第3章
まえおき
今回もSASを使って「データ解析のための統計モデリング入門」を読んでいきます。
実行結果の解説は特に行わないので、内容については書籍をご覧下さい。
データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC (確率と情報の科学)
- 作者: 久保拓弥
- 出版社/メーカー: 岩波書店
- 発売日: 2012/05/19
- メディア: 単行本
- 購入: 16人 クリック: 163回
- この商品を含むブログ (21件) を見る
第3章 一般化線形モデル(GLM) ―ポアソン回帰―
3.2 観測されたデータの概要を調べる
データ読込
%let home_dir = /folders/myfolders/midoribon/chapter03/; data d; infile "&home_dir.data3a.csv" dsd missover firstobs = 2; attrib y length = 8 label = "種子数" x length = 8 label = "体サイズ" f length = $1 label = "施肥処理" ; input y x f; run;
P42 全100個体ぶんのデータを表示
proc print data = d; run;
実行結果
OBS | y | x | f |
---|---|---|---|
1 | 6 | 8.31 | C |
2 | 6 | 9.44 | C |
3 | 6 | 9.50 | C |
...(中略)...
99 | 7 | 10.86 | T |
---|---|---|---|
100 | 9 | 9.97 | T |
P42 列ごとにデータを表示:x列
proc print data = d; var x; run;
実行結果
OBS | x |
---|---|
1 | 8.31 |
2 | 9.44 |
3 | 9.50 |
...(中略)...
99 | 10.86 |
---|---|
100 | 9.97 |
P43 列ごとにデータを表示:y列
proc print data = d; var y; run;
実行結果
OBS | y |
---|---|
1 | 6 |
2 | 6 |
3 | 6 |
...(中略)...
99 | 7 |
---|---|
100 | 9 |
P43 列ごとにデータを表示:f列
proc print data = d; var f; run;
実行結果
OBS | f |
---|---|
1 | C |
2 | C |
3 | C |
...(中略)...
99 | T |
---|---|
100 | T |
P44 データセットの概要:数値変数
proc means data = d min q1 median mean q3 max; run;
実行結果
変数 | ラベル | 最小値 | 下側四分位点 | 中央値 | 平均 | 上側四分位点 | 最大値 |
---|---|---|---|---|---|---|---|
y | 種子数 | 2.0000000 | 6.0000000 | 8.0000000 | 7.8300000 | 10.0000000 | 15.0000000 |
x | 体サイズ | 7.1900000 | 9.4250000 | 10.1550000 | 10.0891000 | 10.7000000 | 12.4000000 |
P44 データセットの概要:文字変数
proc freq data = d; tables f / missing nocol norow nopercent; run;
実行結果
施肥処理 | ||
---|---|---|
f | 度数 | 累積 度数 |
C | 50 | 50 |
T | 50 | 100 |
3.3 統計モデリングの前にデータを図示する
P46 図3.2 例題の架空データの図示
proc sgplot data = d; scatter x = x y = y / group = f markerattrs = (symbol=circlefilled); run;
実行結果
P46 図3.3 植物の種子数の分布を、施肥処理でグループわけした箱ひげ図
proc sgplot data = d; vbox y / category = f; run;
実行結果
3.4 ポアソン回帰の統計モデル
P48 図3.4 個体iの平均種子数と体サイズの関係
data data3_4; do x = -4 to 5 by 0.1; lambda1 = exp(-2 - 0.8 * x); lambda2 = exp(-1 + 0.4 * x); output; end; run; proc sgplot data = data3_4; series x = x y = lambda1 / lineattrs = (pattern=dash); series x = x y = lambda2; refline 0 / axis = x lineattrs = (pattern=dot); xaxis label = "個体iの体のサイズx_i"; yaxis label = "個体iのλ_i"; run;
実行結果
P50 ポアソン回帰
proc genmod data = d; model y = x / dist = poisson link = log; ods output parameterestimates = param_est; run;
実行結果(抜粋)
最大尤度パラメータ推定値の分析 | |||||||
---|---|---|---|---|---|---|---|
パラメータ | 自由度 | 推定値 | 標準誤差 | Wald 95% 信頼限界 | Wald カイ 2 乗 |
Pr > ChiSq | |
Intercept | 1 | 1.2917 | 0.3637 | 0.5789 | 2.0045 | 12.61 | 0.0004 |
x | 1 | 0.0757 | 0.0356 | 0.0059 | 0.1454 | 4.52 | 0.0336 |
Scale | 0 | 1.0000 | 0.0000 | 1.0000 | 1.0000 |
P52 図3.6 パラメーター推定値のばらつきの評価
data _null_; set param_est; if parameter eq "Intercept" then do; call symputx("beta1", estimate); call symputx("std_err1", stderr); end; if parameter eq "x" then do; call symputx("beta2", estimate); call symputx("std_err2", stderr); end; run; data param3_6; do i = -0.1 to 1.6 by 0.005; beta1 = pdf("normal", i, &beta1., &std_err1.); beta2 = pdf("normal", i, &beta2., &std_err2.); output; end; run; proc sgplot data = param3_6; series x = i y = beta1; series x = i y = beta2; refline 0 / axis = x; run;
実行結果
P54 図3.7 平均種子数λの予測
proc sql noprint; select min(x), max(x) into :minx, :maxx from d; quit; data line3_7; do x1 = &minx. to &maxx. by 0.1; y1 = exp(&beta1. + &beta2. * x1); output; end; run; data d3_7; set d line3_7; run; proc sgplot data = d3_7; scatter x = x y = y / group = f markerattrs = (symbol=circlefilled); series x = x1 y = y1 / lineattrs=(color=orange); run;
実行結果
3.5 説明変数が因子型の統計モデル
P56 説明変数が因子型の統計モデル
proc genmod data = d; class f(ref="C"); model y = f / dist = poisson link = log; run;
実行結果(抜粋)
適合度評価の基準 | |||
---|---|---|---|
基準 | 自由度 | 値 | 値/自由度 |
デビアンス | 98 | 89.4750 | 0.9130 |
Scaled デビアンス | 98 | 89.4750 | 0.9130 |
Pearson カイ 2 乗 | 98 | 87.0794 | 0.8886 |
Scaled Pearson カイ 2 乗 | 98 | 87.0794 | 0.8886 |
対数尤度 | 828.4006 | ||
完全対数尤度 | -237.6273 | ||
AIC (小さいほどよい) | 479.2545 | ||
AICC (小さいほどよい) | 479.3782 | ||
BIC (小さいほどよい) | 484.4649 |
最大尤度パラメータ推定値の分析 | ||||||||
---|---|---|---|---|---|---|---|---|
パラメータ | 自由度 | 推定値 | 標準誤差 | Wald 95% 信頼限界 | Wald カイ 2 乗 |
Pr > ChiSq | ||
Intercept | 1 | 2.0516 | 0.0507 | 1.9522 | 2.1509 | 1637.26 | <.0001 | |
f | T | 1 | 0.0128 | 0.0715 | -0.1273 | 0.1529 | 0.03 | 0.8582 |
f | C | 0 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | . | . |
Scale | 0 | 1.0000 | 0.0000 | 1.0000 | 1.0000 |
3.6 説明変数が数量型+因子型の統計モデル
P58 説明変数が数量型+因子型の統計モデル
proc genmod data = d; class f(ref="C"); model y = x f / dist = poisson link = log; run;
実行結果(抜粋)
適合度評価の基準 | |||
---|---|---|---|
基準 | 自由度 | 値 | 値/自由度 |
デビアンス | 97 | 84.8079 | 0.8743 |
Scaled デビアンス | 97 | 84.8079 | 0.8743 |
Pearson カイ 2 乗 | 97 | 83.8215 | 0.8641 |
Scaled Pearson カイ 2 乗 | 97 | 83.8215 | 0.8641 |
対数尤度 | 830.7341 | ||
完全対数尤度 | -235.2937 | ||
AIC (小さいほどよい) | 476.5874 | ||
AICC (小さいほどよい) | 476.8374 | ||
BIC (小さいほどよい) | 484.4029 |
最大尤度パラメータ推定値の分析 | ||||||||
---|---|---|---|---|---|---|---|---|
パラメータ | 自由度 | 推定値 | 標準誤差 | Wald 95% 信頼限界 | Wald カイ 2 乗 |
Pr > ChiSq | ||
Intercept | 1 | 1.2631 | 0.3696 | 0.5386 | 1.9876 | 11.68 | 0.0006 | |
x | 1 | 0.0801 | 0.0370 | 0.0075 | 0.1527 | 4.67 | 0.0306 | |
f | T | 1 | -0.0320 | 0.0744 | -0.1778 | 0.1138 | 0.19 | 0.6670 |
f | C | 0 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | . | . |
Scale | 0 | 1.0000 | 0.0000 | 1.0000 | 1.0000 |
3.7 「何でも正規分布」「何でも直線」には無理がある
データ読込
data d0; infile "&home_dir.d0.csv" dsd missover firstobs = 2; attrib x length = 8 label = "" y length = 8 label = "" ; input x y; run;
P61 図3.9(A) 正規分布・恒等リンク関数の統計モデル
proc genmod data = d0; model y = x / dist = normal link = id; ods output parameterestimates = paramest1; run; data _null_; set paramest1; if parameter eq "Intercept" then call symputx("beta3", estimate); if parameter eq "x" then call symputx("beta4", estimate); run; data line3_9a; do x1 = 0 to 2 by 0.1; y1 = &beta3. + &beta4. * x1; output; end; run; data d0_3_9a; set d0 line3_9a; run; proc sgplot data = d0_3_9a; scatter x = x y = y; series x = x1 y = y1 / lineattrs = (pattern=dash); refline 0 / axis = y; yaxis min = -2; run;
実行結果
P61 図3.9(B) ポアソン分布・対数リンク関数の統計モデル
proc genmod data = d0; model y = x / dist = poisson link = log; ods output parameterestimates = paramest2; run; data _null_; set paramest2; if parameter eq "Intercept" then call symputx("beta5", estimate); if parameter eq "x" then call symputx("beta6", estimate); run; data line3_9b; do x1 = 0 to 2 by 0.1; y1 = exp(&beta5. + &beta6. * x1); output; end; run; data d0_3_9b; set d0 line3_9b; run; proc sgplot data = d0_3_9b; scatter x = x y = y; series x = x1 y = y1 / lineattrs = (pattern=dash); refline 0 / axis = y; yaxis min = -2; run;
実行結果