SASで「データ解析のための統計モデリング入門」を読む 第3章

まえおき

今回もSASを使って「データ解析のための統計モデリング入門」を読んでいきます。
実行結果の解説は特に行わないので、内容については書籍をご覧下さい。

第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;
実行結果

f:id:Prunus1350:20150203003938p:plain

P46 図3.3 植物の種子数の分布を、施肥処理でグループわけした箱ひげ図
proc sgplot data = d;
    vbox y / category = f;
run;
実行結果

f:id:Prunus1350:20150203003939p:plain

3.4 ポアソン回帰の統計モデル

P48 図3.4 個体iの平均種子数{\displaystyle \lambda_i}と体サイズ{\displaystyle x_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;
実行結果

f:id:Prunus1350:20150203003940p:plain

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;
実行結果

f:id:Prunus1350:20150203003941p:plain

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;
実行結果

f:id:Prunus1350:20150203003942p:plain

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;
実行結果

f:id:Prunus1350:20150203003943p:plain

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;
実行結果

f:id:Prunus1350:20150203003944p:plain