SASによる統計解析

SASによる統計解析#

ここでは、SASにおける様々な基本的な統計検定の実行方法について説明します。

  • 比率検定

  • カイ二乗検定

  • フィッシャーの正確確率検定

  • 相関係数

  • t検定/順位和検定

  • 一元配置分散分析 (ANOVA)/Kruskal-Wallis検定

  • 線形回帰分析

  • ロジスティック回帰分析

  • ポアソン回帰分析

注: このコースでは、これらの検定の統計理論や「公式」については詳しく説明しません。 オンラインには、これらの検定について詳しく解説されているものがたくさんあるので、ここで扱っていない分野を勉強したい場合はそちらを参照してください。ここでは、検定に適合するコードを書くことだけが必要とされ、結果の解釈は求められません。

比率検定#

１標本の比率検定を行うには、PROC FREQ を使用できます。この検定を行うには、TABLES ステートメントで BINOMIAL オプションを使用します。BINOMIAL オプションには以下のようなものを指定できます。

  • p= - 仮説検定の帰無仮説値

  • level= - “成功” とみなすグループを指定

  • CORRECT - 連続性補正を使用して p 値を計算します (標本サイズが小さい場合に有用)

  • CL= - Wald、EXACT、LOGIT などの異なるタイプの信頼区間を選択できます。

例#

以下の例では、”successes” のカウントと “failures” のカウントを持つ集計データセットを使用しています。ここでは、喫煙者の比率に興味があるので、喫煙者数と非喫煙者数をカウントしています。

data smoke;
  input smkstatus $ count;
  datalines;
Y 15
N 17
;
run;

proc freq data = smoke;
  tables smkstatus / binomial(p = 0.5 level = "Y" correct) alpha = 0.05;
  weight count;
run;
SAS 出力

SAS システム

FREQ プロシジャ

smkstatus 度数 パーセント 累積
度数		 累積
パーセント
N 17 53.13 17 53.13
Y 15 46.88 32 100.00
二項分布の比率
smkstatus = Y
比率 0.4688
漸近標準誤差 0.0882
95% 信頼下限 0.2802
95% 信頼上限 0.6573
   
正確な信頼限界  
95% 信頼下限 0.2909
95% 信頼上限 0.6526
H0: 母比率 = 0.5 に対する検定
漸近の信頼区間と検定には
連続修正が含まれます。
帰無仮説が正しいもとでの漸近標準誤差 0.0884
Z -0.1768
片側 Pr < Z 0.4298
両側 Pr > |Z| 0.8597

標本サイズ = 32

WEIGHT ステートメントを使用して Y と N のカウントを指定することに注意してください。このステートメントがない場合、 Y が 1、N が 1 オブザベーションと読み込んでしまいます。
推定比率は 0.4688 です。(漸近)95%信頼区間は (0.2802, 0.6573) であり、両側 (連続性補正) p 値は、\(H_0: p=0.5\) vs \(H_a: p\neq 0.5\) を検定するもので、0.8597 です。
あるいは、以下のように、個々のデータが列挙されている形でも構いません。

data smoke2;
  do i = 1 to 15;
    smkstatus = "Y";
    output;
  end;
  do i = 1 to 17;
    smkstatus = "N";
    output;
  end;
  drop i;
run;

proc freq data = smoke2;
  tables smkstatus / binomial(p = 0.5 level = "Y" correct) alpha = 0.05;
run;
SAS 出力

SAS システム

FREQ プロシジャ

smkstatus 度数 パーセント 累積
度数		 累積
パーセント
N 17 53.13 17 53.13
Y 15 46.88 32 100.00
二項分布の比率
smkstatus = Y
比率 0.4688
漸近標準誤差 0.0882
95% 信頼下限 0.2802
95% 信頼上限 0.6573
   
正確な信頼限界  
95% 信頼下限 0.2909
95% 信頼上限 0.6526
H0: 母比率 = 0.5 に対する検定
漸近の信頼区間と検定には
連続修正が含まれます。
帰無仮説が正しいもとでの漸近標準誤差 0.0884
Z -0.1768
片側 Pr < Z 0.4298
両側 Pr > |Z| 0.8597

標本サイズ = 32

カイ二乗検定#

2 つのカテゴリ変数間の関連を検定するには、カイ二乗独立性検定を行うことができます。ここでも、PROC FREQ を使用して tables ステートメントを実行します。2x2 の表の場合は、カイ二乗検定が自動的に実行されますが、より大きな 表の場合は、tables ステートメントに CHISQ オプションを指定することで出力することができます。また、有用なオプションとして、帰無仮説である独立性の下での期待されるセルカウントを提供する EXPECTED オプションもあります。これらの期待されるセルカウントは、カイ二乗検定が適切かどうかを判断するために必要です。

例#

次の例では、Kaggle の中古車オークションのデータセットを使用して、オンラインでの販売と損をしているかどうかの関連を検定しています。

filename cardata '/folders/myfolders/SAS_Notes/data/kaggleCarAuction.csv';

proc import datafile = cardata out = cars dbms = csv replace;
  getnames = yes;
  guessingrows = 1000;
run;

proc freq data = cars;
  tables isbadbuy*isonlinesale / chisq expected;
run;
SAS 出力

SAS システム

FREQ プロシジャ

度数
期待
パーセント
行のパーセント
列のパーセント
表 : IsBadBuy * IsOnlineSale
IsBadBuy IsOnlineSale
0 1 合計
0
62375
62389
85.47
97.45
87.68
1632
1618.1
2.24
2.55
88.46
64007
 
87.70
 
 
1
8763
8749.1
12.01
97.63
12.32
213
226.91
0.29
2.37
11.54
8976
 
12.30
 
 
合計
71138
97.47
1845
2.53
72983
100.00

IsBadBuy * IsOnlineSale の統計量

統計量 自由度 値 p 値
カイ 2 乗値 1 0.9978 0.3178
尤度比カイ 2 乗値 1 1.0154 0.3136
連続性補正カイ 2 乗値 1 0.9274 0.3356
Mantel-Haenszel のカイ 2 乗値 1 0.9978 0.3179
ファイ係数   -0.0037  
一致係数   0.0037  
Cramer の V 統計量   -0.0037  
Fisher の正確検定
セル (1,1) 度数 (F) 62375
左側 Pr <= F 0.1679
右側 Pr >= F 0.8498
   
表の確率 (P) 0.0177
両側 Pr <= P 0.3324

標本サイズ = 72983

カイ二乗検定の結果、p 値は 0.3178 です。連続性補正付きのカイ二乗検定を使用すると、p 値は 0.3356 になります。

2x2 の場合、この例のように、リスク差、相対リスク、オッズ比などの効果量を求めることもできます。これらは、RISKDIFF、RELRISK、OR オプションを使用すると、信頼区間とともに 3 つすべてを求めることができます。

proc freq data = cars;
  tables isbadbuy*isonlinesale / riskdiff relrisk or;
run;
SAS 出力

SAS システム

FREQ プロシジャ

度数
パーセント
行のパーセント
列のパーセント
表 : IsBadBuy * IsOnlineSale
IsBadBuy IsOnlineSale
0 1 合計
0
62375
85.47
97.45
87.68
1632
2.24
2.55
88.46
64007
87.70
 
 
1
8763
12.01
97.63
12.32
213
0.29
2.37
11.54
8976
12.30
 
 
合計
71138
97.47
1845
2.53
72983
100.00

IsBadBuy * IsOnlineSale の統計量

列 1 リスクの推定値
  リスク ASE 95%
信頼限界		 正確 95%
信頼限界
行 1 - 行 2 の差
行 1 0.9745 0.0006 0.9733 0.9757 0.9733 0.9757
行 2 0.9763 0.0016 0.9731 0.9794 0.9729 0.9793
合計 0.9747 0.0006 0.9736 0.9759 0.9736 0.9758
差 -0.0018 0.0017 -0.0051 0.0016    
列 2 リスクの推定値
  リスク ASE 95%
信頼限界		 正確 95%
信頼限界
行 1 - 行 2 の差
行 1 0.0255 0.0006 0.0243 0.0267 0.0243 0.0267
行 2 0.0237 0.0016 0.0206 0.0269 0.0207 0.0271
合計 0.0253 0.0006 0.0241 0.0264 0.0242 0.0264
差 0.0018 0.0017 -0.0016 0.0051    
オッズ比と相対リスク
統計量 値 95% 信頼限界
オッズ比 0.9290 0.8040 1.0735
相対リスク (列 1) 0.9982 0.9947 1.0016
相対リスク (列 2) 1.0745 0.9331 1.2373

標本サイズ = 72983

リスク差については、2 つの表が出力されます。1 つ目の表は、第 1 列の条件付き行比率と第 2 列の条件付き行比率を比較したものです。同様に、相対リスクについては、第 1 列と第 2 列の相対リスクが得られます。これにより、目的のアウトカムに対応する列を選択することができます。

フィッシャーの正確確率検定#

2 つのカテゴリ変数間の関連を検定する別の方法として、フィッシャーの正確確率検定があります。この検定は、ランダムサンプリング以外には何の仮定も置かないノンパラメトリック検定です。ただし、注意しなければならないのは、変数の水準の数と観察回数が増えるにつれて、検定の実行に必要な計算時間が長くなることです。2x2 表の場合ではこの検定の実行時間は非常に短いですが、5x5表の場合は、データ量や使用するコンピュータによっては、完了までに数時間かかる場合もあります。

2x2 表の場合はこの検定が出力され、それより大きい表の場合はTableステートメントにFISHERオプションを使用すると検定を行うことができます。

例#

次のプログラムは、フィッシャーの正確確率検定を使用して、車の購入が損であったかと、その車をオンラインで購入することの関連について検定します。

proc freq data = cars;
  tables isbadbuy*isonlinesale / fisher;
run;
SAS 出力

SAS システム

FREQ プロシジャ

度数
パーセント
行のパーセント
列のパーセント
表 : IsBadBuy * IsOnlineSale
IsBadBuy IsOnlineSale
0 1 合計
0
62375
85.47
97.45
87.68
1632
2.24
2.55
88.46
64007
87.70
 
 
1
8763
12.01
97.63
12.32
213
0.29
2.37
11.54
8976
12.30
 
 
合計
71138
97.47
1845
2.53
72983
100.00

IsBadBuy * IsOnlineSale の統計量

統計量 自由度 値 p 値
カイ 2 乗値 1 0.9978 0.3178
尤度比カイ 2 乗値 1 1.0154 0.3136
連続性補正カイ 2 乗値 1 0.9274 0.3356
Mantel-Haenszel のカイ 2 乗値 1 0.9978 0.3179
ファイ係数   -0.0037  
一致係数   0.0037  
Cramer の V 統計量   -0.0037  
Fisher の正確検定
セル (1,1) 度数 (F) 62375
左側 Pr <= F 0.1679
右側 Pr >= F 0.8498
   
表の確率 (P) 0.0177
両側 Pr <= P 0.3324

標本サイズ = 72983

フィッシャーの正確確率検定によるP値は0.3324です。

相関分析#

CORR プロシージャは、パラメトリックなピアソン相関係数とノンパラメトリックなスピアマン順位相関係数、および仮説検定の両方を実行して、相関分析を行うことができます。デフォルトの相関出力はピアソン相関係数です。スピアマン順位相関係数を要求するには、PROC CORR ステートメントに SPEARMAN オプションを追加します。

例#

Charm City 循環バス乗車率データセットを使用して、PROC CORR を用いたいくつかの例を見てみましょう。次のプログラムは、オレンジラインとパープルラインのバス間の平均乗車人数の相関関係と仮説検定結果を出力します。

filename busdata '/folders/myfolders/SAS_Notes/data/Charm_City_Circulator_Ridership.csv';

proc import datafile = busdata out = circ dbms = csv replace;
  getnames = yes;
  guessingrows = 1000;
run;

proc corr data = circ;
  var orangeaverage purpleaverage;
run;
SAS 出力

SAS システム

CORR プロシジャ

2 変数 : orangeAverage purpleAverage
単純統計量
変数 N 平均 標準偏差 合計 最小値 最大値
orangeAverage 1136 3033 1228 3445671 0 6927
purpleAverage 993 4017 1407 3988816 0 8090
Pearson の相関係数
H0: Rho=0 に対する Prob > |r|
オブザベーション数
  orangeAverage purpleAverage
orangeAverage
1.00000
 
1136
0.91954
<.0001
993
purpleAverage
0.91954
<.0001
993
1.00000
 
993

例#

また、複数の変数の相関行列を一度に取得することもできます。次の例でも、相関係数を計算する際にペアワイズ完全ではなく完全なオブザベーションのみを使用するために、NOMISS オプションを使用しています。ここでは、オレンジ、パープル、バナー、グリーンの 4 つのバス路線間の平均乗車人数間の相関行列を算出します。

proc corr data = circ nomiss;
  var orangeAverage purpleAverage greenAverage bannerAverage;
run;
SAS 出力

SAS システム

CORR プロシジャ

4 変数 : orangeAverage purpleAverage greenAverage bannerAverage
単純統計量
変数 N 平均 標準偏差 合計 最小値 最大値
orangeAverage 270 3859 1095 1041890 0 6927
purpleAverage 270 4552 1297 1228935 0 8090
greenAverage 270 2090 556.00353 564213 0 3879
bannerAverage 270 827.26852 436.04872 223363 0 4617
Pearson の相関係数, N = 270
H0: Rho=0 に対する Prob > |r|
  orangeAverage purpleAverage greenAverage bannerAverage
orangeAverage
1.00000
 
0.90788
<.0001
0.83958
<.0001
0.54470
<.0001
purpleAverage
0.90788
<.0001
1.00000
 
0.86656
<.0001
0.52135
<.0001
greenAverage
0.83958
<.0001
0.86656
<.0001
1.00000
 
0.45334
<.0001
bannerAverage
0.54470
<.0001
0.52135
<.0001
0.45334
<.0001
1.00000
 

すべてのペアワイズ相関ではなく、特定のペアのみが必要な場合は、次の例のように WITH ステートメントを使用できます。

proc corr data = circ nomiss;
  var orangeAverage purpleAverage;
  with greenAverage bannerAverage;
run;
SAS 出力

SAS システム

CORR プロシジャ

2 With 変数 : greenAverage bannerAverage
2 変数 : orangeAverage purpleAverage
単純統計量
変数 N 平均 標準偏差 合計 最小値 最大値
greenAverage 270 2090 556.00353 564213 0 3879
bannerAverage 270 827.26852 436.04872 223363 0 4617
orangeAverage 270 3859 1095 1041890 0 6927
purpleAverage 270 4552 1297 1228935 0 8090
Pearson の相関係数, N = 270
H0: Rho=0 に対する Prob > |r|
  orangeAverage purpleAverage
greenAverage
0.83958
<.0001
0.86656
<.0001
bannerAverage
0.54470
<.0001
0.52135
<.0001

ピアソン相関係数ではなくスピアマン順位相関係数を取得するには、PROC CORR ステートメントに SPEARMAN オプションを追加します。

例#

次のプログラムは、オレンジラインとパープルラインの平均乗車人数間の相関係数が 0 であるという仮説検定に関連するスピアマン順位相関係数と p 値を出力します。

proc corr data = circ spearman;
  var orangeAverage purpleAverage;
run;
SAS 出力

SAS システム

CORR プロシジャ

2 変数 : orangeAverage purpleAverage
単純統計量
変数 N 平均 標準偏差 中央値 最小値 最大値
orangeAverage 1136 3033 1228 2968 0 6927
purpleAverage 993 4017 1407 4223 0 8090
Spearman の相関係数
H0: Rho=0 に対する Prob > |r|
オブザベーション数
  orangeAverage purpleAverage
orangeAverage
1.00000
 
1136
0.91455
<.0001
993
purpleAverage
0.91455
<.0001
993
1.00000
 
993

t検定#

SASでは、TTEST プロシージャを使用して以下の t 検定を実行できます。

  • 一標本 t 検定

  • 対応のある標本 t 検定

  • 二標本 t 検定

例#

この例では、一標本 t 検定を使用して、オレンジラインの平均乗車人数が 3000 を超えているかどうかを検定します。

proc ttest data = circ h0 = 3000 side = u;
  var orangeAverage;
run;
SAS 出力

SAS システム

TTEST プロシジャ

 

変数 : orangeAverage

N 平均 標準偏差 標準誤差 最小値 最大値
1136 3033.2 1227.6 36.4217 0 6926.5
平均 平均の
95% 信頼限界		 標準偏差 標準偏差の
95% 信頼限界
3033.2 2973.2 Infty 1227.6 1179.1 1280.3
自由度 t 値 Pr > t
1135 0.91 0.1814
orangeAverage の要約パネル
orangeAverage の QQ プロット

H0= オプションは t 検定における帰無仮説の値を指定し、SIDE= オプションは片側 (L)、両側 (U)、または不等号 (2) の検定を行うかどうかを指定します。デフォルト値は、帰無仮説の値が 0、片側検定が 2 (両側) です。出力には、いくつかの要約統計量、検定の p 値、信頼区間、ヒストグラム、および QQ プロット (正規性の仮定を評価するため) が含まれています。

出力を見ると、p 値は 0.1814 であることがわかります。片側検定を要求したので、片側の信頼区間が得られます。通常の (両側) 信頼区間を得るには、両側検定を要求する必要があります。

二標本 t 検定を行うには、データが 2 列の形式になっている必要があります。

  • 両方のグループの定量データを含むデータ列

  • グループを示すグループ変数列

PROC TTEST では、データ変数を VAR ステートメントに、グループ変数を CLASS ステートメントに指定して、二標本 t 検定を行います。

例#

次のプログラムでは、オレンジラインとパープルラインの平均乗車人数間で二標本 t 検定を実行します。最初に、データを PROC TTEST に必要なデータ形式に変換する必要があります。

data circ_sub;
  set circ;
  count = orangeaverage;
  group = "orange";
  output;
  count = purpleaverage;
  group = "purple";
  output;
  keep count group;
run;

proc ttest data = circ_sub;
  var count;
  class group;
run;
SAS 出力

SAS システム

TTEST プロシジャ

 

変数 : count

group 手法 N 平均 標準偏差 標準誤差 最小値 最大値
orange   1136 3033.2 1227.6 36.4217 0 6926.5
purple   993 4016.9 1406.7 44.6388 0 8089.5
Diff (1-2) Pooled   -983.8 1314.1 57.0906    
Diff (1-2) Satterthwaite   -983.8   57.6122    
group 手法 平均 平均の
95% 信頼限界		 標準偏差 標準偏差の
95% 信頼限界
orange   3033.2 2961.7 3104.6 1227.6 1179.1 1280.3
purple   4016.9 3929.3 4104.5 1406.7 1347.4 1471.4
Diff (1-2) Pooled -983.8 -1095.7 -871.8 1314.1 1275.8 1354.9
Diff (1-2) Satterthwaite -983.8 -1096.8 -870.8      
手法 分散 自由度 t 値 Pr > |t|
Pooled Equal 2127 -17.23 <.0001
Satterthwaite Unequal 1984 -17.08 <.0001
等分散性
手法 分子の自由度 分母の自由度 F 値 Pr > F
Folded F 992 1135 1.31 <.0001
count の要約パネル
count の QQ プロット

出力には、各グループの要約統計量、各グループ平均の信頼区間、2 つの平均の差の信頼区間、2 つの平均の差の仮説検定、分散の等性に関する F 検定が含まれます。Pooled 行は、2 標本 t 検定 (2 つのグループ間の母集団分散が等しいことを仮定) に対応し、Satterthwaite は、母集団分散が等しくないことを仮定しています。

ここで使用しているデータは実際には対応のあるデータであることに注意してください。これは、2 つのバス路線間の日付ごとの平均乗車人数の行が一致しているためです。次は、対応のある標本 t 検定について見ていきます。

対応のある標本 t 検定を行うには、PAIRED ステートメントを使用する必要があります。この場合、各グループからのデータが 2 つの別々の列にあり、同じ行のオブザベーションが一致するペアと仮定されます。

例#

以下のプログラムはオレンジとパープルのバス路線の平均乗客数で対応のあるt検定を行います。

proc ttest data = circ;
  paired orangeAverage*purpleAverage;
run;
SAS 出力

SAS システム

TTEST プロシジャ

 

差 : orangeAverage - purpleAverage

N 平均 標準偏差 標準誤差 最小値 最大値
993 -764.1 572.3 18.1613 -2998.0 2504.5
平均 平均の
95% 信頼限界		 標準偏差 標準偏差の
95% 信頼限界
-764.1 -799.8 -728.5 572.3 548.2 598.6
自由度 t 値 Pr > |t|
992 -42.08 <.0001
orangeAverage と purpleAverage の差の要約パネル
orangeAverage と purpleAverage のプロファイルプロット
orangeAverage と purpleAverage の一致度プロット
orangeAverage と purpleAverage の差の QQ プロット

ノンパラメトリック t 検定の代替#

サンプルサイズが小さく、データが正規分布に従っていないと仮定できない場合、ノンパラメトリック検定を使用する必要があります。t 検定の代替として、以下の検定が利用できます。

  • 一標本 t 検定または対応のある t 検定の代替として、符号検定またはウィルコクソン符号順位検定

  • 二標本 t 検定の代替として、ウィルコクソン順位和検定

ウィルコクソン順位和検定を行うには、PROC NPAR1WAY を使用します。

例#

以下の例では、PROC NPAR1WAY を使用してウィルコクソン順位和検定を行い、オレンジとパープルのバス路線の乗客数の中央値を比較します。

proc npar1way data = circ_sub wilcoxon;
  var count;
  class group;
run;
SAS 出力

SAS システム

NPAR1WAY プロシジャ

変数 count に対する Wilcoxon スコア (順位和)
分類変数 : group
group N スコアの
合計		 H0 のもとでの
期待値		 H0 のもとでの
標準偏差		 平均
スコア
同順位には平均スコアを使用しました。
orange 1136 982529.50 1209840.0 14150.2115 864.90273
purple 993 1284855.50 1057545.0 14150.2115 1293.91289
Wilcoxon の順位和検定 (2 標本)
統計量 Z Pr > Z Pr > |Z| t 近似
Pr > Z Pr > |Z|
Z には 0.5 の連続性の補正が含まれています。
1284856 16.0641 <.0001 <.0001 <.0001 <.0001
Kruskal-Wallis 検定
カイ 2 乗 自由度 Pr > ChiSq
258.0555 1 <.0001
group で分類された count に対する Wilcoxon スコアの箱ひげ図

符号検定またはウィルコクソン符号順位検定を行うには、まず対応のあるペア間の差分を計算し、その差分を PROC UNIVARIATE に渡す必要があります。

例#

以下のプログラムは、PROC UNIVARIATE を使用して、オレンジとパープルのバス路線の乗客数の差の中央値を検定するための符号検定とウィルコクソン符号順位検定の p 値を算出します。

data circ_diff;
  set circ;
  diff = orangeAverage - purpleAverage;
run;

proc univariate data = circ_diff;
  var diff;
run;
SAS 出力

SAS システム

UNIVARIATE プロシジャ

変数 : diff

モーメント
N 993 重み変数の合計 993
平均 -764.14401 合計 -758795
標準偏差 572.297901 分散 327524.888
歪度 0.73397459 尖度 2.88082288
無修正平方和 904733341 修正済平方和 324904688
変動係数 -74.893985 平均の標準誤差 18.1613249
基本統計量
位置 ばらつき
平均 -764.144 標準偏差 572.29790
中央値 -788.000 分散 327525
最頻値 0.000 範囲 5503
    四分位範囲 732.50000
 位置の検定 H0: Mu0=0
検定 統計量 p 値
Student の t 検定 t -42.0753 Pr > |t| <.0001
符号検定 M -424.5 Pr >= |M| <.0001
符号付順位検定 S -227467 Pr >= |S| <.0001
 分位点 (定義 5)
水準 分位点
100% 最大値 2504.5
99% 912.5
95% 131.5
90% -72.0
75% Q3 -426.5
50% 中央値 -788.0
25% Q1 -1159.0
10% -1433.0
5% -1596.0
1% -1980.0
0% 最小値 -2998.0
極値
最小値 最大値
値 Obs 値 Obs
-2998.0 552 1360.5 171
-2649.5 999 1418.5 672
-2365.0 635 1933.5 743
-2300.5 553 2484.0 674
-2049.5 188 2504.5 709
欠損値
欠損値 カウント パーセント
全体 欠損値
. 153 13.35 100.00

PROC UNIVARIATE はデフォルトで多くの出力をします。符号検定と符号順位検定のp値は位置母数検定表で見つけることができます。

一元配置分散分析 (ANOVA) とクラスカル-ウォリス検定#

2 つを超える独立なグループ間の平均値を比較したい場合は、分散分析 (ANOVA) またはサンプルサイズが少ない場合はクラスカル-ウォリス検定を行うことができます。SASでの一元配置分散分析は、PROC GLM を使用して実行できます。

例#

以下のプログラムは、オレンジ、パープル、グリーンの 3 つのバス路線間の乗客数の平均値が等しいことを検定するために、一元配置分散分析を行います。

data circ_aov;
  set circ;
  count = orangeaverage;
  group = "orange";
  output;
  count = purpleaverage;
  group = "purple";
  output;
  count = greenaverage;
  group = "green";
  output;
  keep count group;
run;

proc glm data = circ_aov;
  class group;
  model count = group;
run;
SAS 出力

SAS システム

GLM プロシジャ

分類変数の水準の情報
分類 水準 値
group 3 green orange purple
読み込んだオブザベーション数 3438
使用されたオブザベーション数 2614

SAS システム

GLM プロシジャ

 

従属変数 : count

要因 自由度 平方和 平均平方 F 値 Pr > F
Model 2 1442596863 721298432 490.02 <.0001
Error 2611 3843371488 1471992    
Corrected Total 2613 5285968351      
R2 乗 変動係数 Root MSE count の平均
0.272911 37.82740 1213.257 3207.349
要因 自由度 Type I
平方和		 平均平方 F 値 Pr > F
group 2 1442596863 721298432 490.02 <.0001
要因 自由度 Type III
平方和		 平均平方 F 値 Pr > F
group 2 1442596863 721298432 490.02 <.0001
group ごとの count の当てはめプロット

代わりにクラスカル-ウォリス検定を行う場合は、PROC NPAR1WAY を使用します。

proc npar1way data = circ_aov wilcoxon;
  class group;
  var count;
run;
SAS 出力

SAS システム

NPAR1WAY プロシジャ

変数 count に対する Wilcoxon スコア (順位和)
分類変数 : group
group N スコアの
合計		 H0 のもとでの
期待値		 H0 のもとでの
標準偏差		 平均
スコア
同順位には平均スコアを使用しました。
orange 1136 1395578.00 1485320.00 19128.0882 1228.50176
purple 993 1720911.50 1298347.50 18728.8588 1733.04280
green 485 301315.50 634137.50 15000.4361 621.26907
Kruskal-Wallis 検定
カイ 2 乗 自由度 Pr > ChiSq
729.0668 2 <.0001
group で分類された count に対する Wilcoxon スコアの箱ひげ図

線形回帰#

線形回帰モデルを使用するには 2 つのプロシージャを使用できます。

  • PROC REG

  • PROC GLM

PROC REG は標準出力の大部分を生成しますが、MODELステートメントでは相互作用項などすべての変数を事前にデータステップで準備する必要があります。一方PROC GLM は、MODEL ステートメント内で相互作用項を計算できます。一般的に、PROC GLM を好まれますが、どちらのプロシージャでも標準的な回帰出力を得ることができます。

PROC REG を PROC GLM よりも好まないもう一つの理由は、PROC REG には CLASS ステートメントがないため、PROC REG を使用するときは、カテゴリ変数を0/1にするダミー変数化をデータステップで手動で行う必要があるからです。

以下に、両方のプロシージャを使用したいくつかの例を示します。

例#

最初の例では、単一の２値の予測変数を持つ単純な線形回帰モデルを使用します。この場合、傾きに関する t 検定は二標本 t 検定と同等であることに注意してください。

data circ_sub;
  set circ_sub;
  if group = "orange" then grp_bin = 1;
  else grp_bin = 0;
run;

proc reg data = circ_sub;
  model count = grp_bin;
run;
SAS 出力

SAS システム

REG プロシジャ

モデル : MODEL1

従属変数 : count

読み込んだオブザベーション数 2292
使用されたオブザベーション数 2129
欠損値を含むオブザベーション数 163
分散分析
要因 自由度 平方和 平均平方 F 値 Pr > F
Model 1 512793031 512793031 296.93 <.0001
Error 2127 3673232559 1726955    
Corrected Total 2128 4186025589      
Root MSE 1314.13647 R2 乗 0.1225
従属変数の平均 3492.00892 調整済み R2 乗 0.1221
変動係数 37.63268    
パラメータの推定
変数 自由度 パラメータ
推定値		 標準誤差 t 値 Pr > |t|
Intercept 1 4016.93454 41.70286 96.32 <.0001
grp_bin 1 -983.77345 57.09059 -17.23 <.0001

SAS システム

REG プロシジャ

モデル : MODEL1

従属変数 : count

count に対する適合度診断パネル
grp_bin ごとの count に対する残差の散布図
Scatterplot of count by grp_bin overlaid with the fit line, a 95% confidence band and lower and upper 95% prediction limits. 
proc glm data = circ_sub plots=diagnostics;
  class group(ref = 'purple');
  model count = group / solution;
run;
SAS 出力

SAS システム

GLM プロシジャ

分類変数の水準の情報
分類 水準 値
group 2 orange purple
読み込んだオブザベーション数 2292
使用されたオブザベーション数 2129

SAS システム

GLM プロシジャ

 

従属変数 : count

要因 自由度 平方和 平均平方 F 値 Pr > F
Model 1 512793031 512793031 296.93 <.0001
Error 2127 3673232559 1726955    
Corrected Total 2128 4186025589      
R2 乗 変動係数 Root MSE count の平均
0.122501 37.63268 1314.136 3492.009
要因 自由度 Type I
平方和		 平均平方 F 値 Pr > F
group 1 512793030.6 512793030.6 296.93 <.0001
要因 自由度 Type III
平方和		 平均平方 F 値 Pr > F
group 1 512793030.6 512793030.6 296.93 <.0001
パラメータ 推定値   標準誤差 t 値 Pr > |t|
Intercept 4016.934542 B 41.70286032 96.32 <.0001
group orange -983.773450 B 57.09058700 -17.23 <.0001
group purple 0.000000 B . . .

Note:X'Xは特異行列です。w正規方程式には一般化逆行列が使用されています。w文字'B'が付けれられた推定値は一意的な推定値ではありません。

count に対する適合度診断パネル
group ごとの count の当てはめプロット

両方のプロシージャから同じ出力を得られることがわかりますが、PROC REG ではグループ変数を手動で 0/1 のダミー変数としてコード化する必要がありました。一方、PROC GLM では、CLASS ステートメントと ref= ステートメントを使用して参照カテゴリを選択できました。また、PROC GLM では、残差診断プロットはデフォルト出力ではないため、指定する必要があります。

PROC GLM と PROC REG の両方から同じ出力を得る方法を確認したので、以降の例では PROC GLM を使用して、ダミー変数と相互作用項を計算するためのデータステップは使用しません。

例#

以下のプログラムでは、Kaggle の中古車オークションデータセットを使用して、複数の予測変数を持つ線形回帰モデルを使用します。まず、単純な線形回帰モデルを使用し、変数を追加していきます。

proc glm data = cars;
  model vehodo = vehicleage;
run;
SAS 出力

SAS システム

GLM プロシジャ

読み込んだオブザベーション数 72983
使用されたオブザベーション数 72983

SAS システム

GLM プロシジャ

 

従属変数 : VehOdo

要因 自由度 平方和 平均平方 F 値 Pr > F
Model 1 1.5863773E12 1.5863773E12 8313.88 <.0001
Error 72981 1.3925561E13 190810767.21    
Corrected Total 72982 1.5511938E13      
R2 乗 変動係数 Root MSE VehOdo の平均
0.102268 19.31948 13813.43 71500.00
要因 自由度 Type I
平方和		 平均平方 F 値 Pr > F
VehicleAge 1 1.5863773E12 1.5863773E12 8313.88 <.0001
要因 自由度 Type III
平方和		 平均平方 F 値 Pr > F
VehicleAge 1 1.5863773E12 1.5863773E12 8313.88 <.0001
パラメータ 推定値 標準誤差 t 値 Pr > |t|
Intercept 60127.24037 134.8017988 446.04 <.0001
VehicleAge 2722.94117 29.8632078 91.18 <.0001

次に、別の変数「IsBadBuy」を追加します。この変数は既に 0/1 のダミー変数であるため、クラス ステートメントに入れる必要はありません (ただし、一致する参照カテゴリを選択すれば、CLASS ステートメントに指定しても出力を得ることができます)。

proc glm data = cars;
  model VehOdo = VehicleAge IsBadBuy;
run;
SAS 出力

SAS システム

GLM プロシジャ

読み込んだオブザベーション数 72983
使用されたオブザベーション数 72983

SAS システム

GLM プロシジャ

 

従属変数 : VehOdo

要因 自由度 平方和 平均平方 F 値 Pr > F
Model 2 1.5998928E12 799946379347 4196.37 <.0001
Error 72980 1.3912045E13 190628187.46    
Corrected Total 72982 1.5511938E13      
R2 乗 変動係数 Root MSE VehOdo の平均
0.103139 19.31023 13806.82 71500.00
要因 自由度 Type I
平方和		 平均平方 F 値 Pr > F
VehicleAge 1 1.5863773E12 1.5863773E12 8321.84 <.0001
IsBadBuy 1 13515481118 13515481118 70.90 <.0001
要因 自由度 Type III
平方和		 平均平方 F 値 Pr > F
VehicleAge 1 1.4941599E12 1.4941599E12 7838.08 <.0001
IsBadBuy 1 13515481118 13515481118 70.90 <.0001
パラメータ 推定値 標準誤差 t 値 Pr > |t|
Intercept 60141.77139 134.7483412 446.33 <.0001
VehicleAge 2680.32758 30.2749125 88.53 <.0001
IsBadBuy 1329.00242 157.8350951 8.42 <.0001

MODEL ステートメントで複数の予測変数を追加する場合、 + 記号ではなくスペースで区切ります。相互作用項を追加するには、掛け算を使用して個々の相互作用項を作成して主効果も含めるか、または省略記号 | を使用して主効果と相互作用項を一度にすべて指定することができます。

proc glm data = cars;
  model VehOdo = VehicleAge IsBadBuy VehicleAge*IsBadBuy;
  *MODEL VehOdo = VehicleAge|IsBadBuy;
run;
SAS 出力

SAS システム

GLM プロシジャ

読み込んだオブザベーション数 72983
使用されたオブザベーション数 72983

SAS システム

GLM プロシジャ

 

従属変数 : VehOdo

要因 自由度 平方和 平均平方 F 値 Pr > F
Model 3 1.5998931E12 533297702109 2797.54 <.0001
Error 72979 1.3912045E13 190630794.79    
Corrected Total 72982 1.5511938E13      
R2 乗 変動係数 Root MSE VehOdo の平均
0.103139 19.31037 13806.91 71500.00
要因 自由度 Type I
平方和		 平均平方 F 値 Pr > F
VehicleAge 1 1.5863773E12 1.5863773E12 8321.73 <.0001
IsBadBuy 1 13515481118 13515481118 70.90 <.0001
VehicleAge*IsBadBuy 1 347632.54102 347632.54102 0.00 0.9659
要因 自由度 Type III
平方和		 平均平方 F 値 Pr > F
VehicleAge 1 1.2937202E12 1.2937202E12 6786.52 <.0001
IsBadBuy 1 1662398367 1662398367 8.72 0.0031
VehicleAge*IsBadBuy 1 347632.54068 347632.54068 0.00 0.9659
パラメータ 推定値 標準誤差 t 値 Pr > |t|
Intercept 60139.69756 143.2332682 419.87 <.0001
VehicleAge 2680.83719 32.5421914 82.38 <.0001
IsBadBuy 1347.28218 456.2338897 2.95 0.0031
VehicleAge*IsBadBuy -3.78953 88.7403782 -0.04 0.9659

残差と予測値を算出するには、OUTPUT ステートメントを使用します。新しいデータに行を追加して応答変数を欠損値に設定することでも、新しい予測値を算出できます。

例#

以下の例ではOUTPUT ステートメントを使用し残差と予測値を算出します。また予測値を算出するためを新しい列をデータセットに追加します。

data new;
  input VehOdo VehicleAge IsBadBuy;
  datalines;
. 6 1
. 5 0
;
RUN;

data cars_new;
  set cars(keep = VehOdo VehicleAge IsBadBuy) 
      new;
run;

proc glm data = cars_new noprint;
  model VehOdo = VehicleAge IsBadBuy VehicleAge*IsBadBuy;
  output out=res_pred residuals = resid predicted = fitted;
run;

proc sgplot data = res_pred;
  histogram resid;
run;

proc print data = res_pred (firstobs=72984);
run;
SAS 出力
SGPlot プロシジャ

SAS システム

OBS IsBadBuy VehicleAge VehOdo resid fitted
72984 1 6 . . 77549.27
72985 0 5 . . 73543.88

追加のオブザベーションにおける応答変数「VehOdo」の欠損値は、これらの行をモデルに使用できなくなりますが、モデル内の全ての予測変数に値があるため、OUTPUTデータセットでは予測値が計算されます。予測値は、予測変数の値を使用した回帰式に代入することで求められることを思い出してください。例として、追加データの最初のオブザベーションの場合は次のようになります:
$\(\widehat{y}=60139.7 + 1347.28 + 2680.84*6 -3.79*6=77549.28\)$

ロジスティック回帰#

一般化線形モデル (GLM) は、連続/正規でない従属変数の回帰分析を行うことができます。glm の構文は、lm コマンドに似ています。ロジスティック回帰はその一例です。

(単純な) ロジスティック回帰モデルでは、0/1の応答変数 Y と予測変数 x を持ちます。予測変数 x が与えられたとき、\(Y\sim\text{Bernoulli(p(x))}\)、ここで \(p(x)=P(Y=1|x)\) であり、

\[\log\left(\dfrac{P(Y=1|x)}{1-P(Y=1|x)}\right)=\beta_0+\beta_1x\]

となります。つまり、成功の対数オッズが x に従って線形に変化するということです。したがって、\(e^{\beta_1}\) は、x が 1 単位増加したときの成功のオッズ比であることがわかります。

SAS では、ロジスティック回帰モデルに 2 つのプロシージャを使用できます。

  • PROC LOGISTIC

  • PROC GENMOD

一般的には、PROC LOGISTIC を使用します。PROC LOGISTIC はロジスティック回帰専用のプロシージャであり、PROC GENMOD にはない多くの機能を提供するからです。

例#

以下の例では、PROC LOGISTIC を使用して、「IsBadBuy」 を0/1応答変数、「VehOdo」 と 「VehicleAge」 を予測変数とするロジスティック回帰モデルを使用します。

proc logistic data = cars;
  model isbadbuy(event='1') = vehodo vehicleage / clparm=wald clodds=wald;
run;
SAS 出力

SAS システム

LOGISTIC プロシジャ

モデルの情報
データセット WORK.CARS
応答変数 IsBadBuy
応答の水準数 2
モデル binary logit
最適化の手法 Fisher's scoring
読み込んだオブザベーション数 72983
使用されたオブザベーション数 72983
反応プロファイル
順番 IsBadBuy 度数の合計
1 0 64007
2 1 8976

モデルの確率基準は IsBadBuy='1' です。

モデル収束状態
収束基準(GCONV=1E-8)は満たされました。
モデルの適合度統計量
基準 切片のみ 切片と共変量
AIC 54423.307 52352.263
SC 54432.505 52379.857
-2 Log L 54421.307 52346.263
包括的帰無仮説: BETA=0 の検定
検定 カイ 2 乗値 自由度 Pr > ChiSq
尤度比 2075.0443 2 <.0001
スコア 2108.2791 2 <.0001
Wald 2025.3779 2 <.0001
最尤推定値の分析
パラメータ 自由度 推定値 標準誤差 Wald
カイ 2 乗値		 Pr > ChiSq
Intercept 1 -3.7782 0.0638 3505.9533 <.0001
VehOdo 1 8.341E-6 8.526E-7 95.6991 <.0001
VehicleAge 1 0.2681 0.00677 1567.3289 <.0001
予測確率と観測データの応答との関連性
一致の割合 64.4 Somers の D 0.288
不一致の割合 35.6 ガンマ 0.288
タイの割合 0.0 Tau-a 0.062
組 574526832 c 0.644
パラメータ推定と Wald による信頼区間
パラメータ 推定値 95% 信頼限界
Intercept -3.7782 -3.9033 -3.6531
VehOdo 8.341E-6 6.67E-6 0.000010
VehicleAge 0.2681 0.2548 0.2814
オッズ比推定と Wald による信頼区間
効果 単位 推定値 95% 信頼限界
VehOdo 1.0000 1.000 1.000 1.000
VehicleAge 1.0000 1.307 1.290 1.325
プロット : オッズ比の 95% Wald 信頼限界

CLPARM= オプションと CLODDS= オプションは、パラメータ推定値と対応するオッズ比の信頼区間を出力します。

ポアソン回帰#

ポアソン回帰は、カウントデータに対して用いられます。このモデルは (単一の予測変数の場合は) \(Y|x\sim\text{Poisson}(\lambda(x))\)、ここで \(\lambda(x)=E[Y|x]\) であると仮定し、単一の予測変数の場合は

\[\log(E[Y|x])=\beta_0+\beta_1x.\]

となります。このとき、\(e^{\beta_1}\) は、x が 1 単位増加したときの率比を表します。このようなモデルを使用するには、PROC GENMOD を使用します。

例#

以下のプログラムは、曜日を予測変数として、オレンジのバス路線の乗客数のカウントに対してポアソン回帰モデルを使用します。

proc genmod data = circ;
  class day(ref='Friday');
  model orangeBoardings = day / dist = Poisson link = log;
run;
SAS 出力

SAS システム

GENMOD プロシジャ

モデルの情報
データセット WORK.CIRC
分布 Poisson
リンク関数 Log
従属変数 orangeBoardings
読み込んだオブザベーション数 1146
使用されたオブザベーション数 1079
欠損値の数 67
分類変数の水準の情報
分類 水準 値
day 7 Monday Saturday Sunday Thursday Tuesday Wednesday Friday
適合度評価の基準
基準 自由度 値 値/自由度
デビアンス 1072 465776.5310 434.4930
Scaled デビアンス 1072 465776.5310 434.4930
Pearson カイ 2 乗 1072 425148.2927 396.5936
Scaled Pearson カイ 2 乗 1072 425148.2927 396.5936
対数尤度   23002386.843  
完全対数尤度   -238139.4730  
AIC (小さいほどよい)   476292.9459  
AICC (小さいほどよい)   476293.0505  
BIC (小さいほどよい)   476327.8324  
アルゴリズムは収束しました。
最大尤度パラメータ推定値の分析
パラメータ   自由度 推定値 標準誤差 Wald 95% 信頼限界 Wald
カイ 2 乗		 Pr > ChiSq
Intercept   1 8.2279 0.0013 8.2254 8.2305 3.954E7 <.0001
day Monday 1 -0.1964 0.0019 -0.2002 -0.1926 10163.2 <.0001
day Saturday 1 -0.2691 0.0020 -0.2730 -0.2652 18119.2 <.0001
day Sunday 1 -0.6897 0.0023 -0.6942 -0.6852 90824.2 <.0001
day Thursday 1 -0.1523 0.0019 -0.1561 -0.1485 6213.44 <.0001
day Tuesday 1 -0.1719 0.0019 -0.1757 -0.1681 7860.04 <.0001
day Wednesday 1 -0.1396 0.0019 -0.1434 -0.1359 5260.55 <.0001
day Friday 0 0.0000 0.0000 0.0000 0.0000 . .
尺度   0 1.0000 0.0000 1.0000 1.0000    

Note:尺度パラメータは固定されています。

model ステートメントでは、dist=Poisson オプションを指定して、応答変数がポアソン分布に従うと仮定し、link=log オプションで log リンクを使用することを指定します。
ポアソン回帰モデルでオフセットを使用したい場合は、model ステートメントで OFFSET= オプションを使用できます。ただし、このオプションを使用する場合は、オフセット値の対数 (log()) を自分で設定する必要があります。

演習#

これらの演習では、乳児死亡率データセット (indicatordeadkids35.csv) と Kaggle 中古車オークションデータセット (kaggleCarAuction.csv) を使用します。以下のコードを修正して、このデータセットを読み込んでください。

filename cardata '/folders/myfolders/SAS_Notes/data/kaggleCarAuction.csv';

proc import datafile = cardata out = cars dbms = csv replace;
  getnames = yes;
  guessingrows = 1000;
run;

filename mortdat '/folders/myfolders/SAS_Notes/data/indicatordeadkids35.csv';

proc import datafile = mortdat out = mort dbms = csv replace;
  getnames = yes;
  guessingrows = 500;
run;

  1. 1980年、1990年、2000年、2010年の死亡率データ間の相関関係を計算します。結果を画面に表示するだけで十分です。次に、NOMMISS オプションを使用して計算します。(注意: 列名は数値ですが、SASの標準的な名前としては無効なので、コード内で 1980 年の変数を指すには ‘1980’n を使用します)

a. ミャンマー、中国、およびアメリカの死亡率データ間の相関関係を計算します。この相関行列をODS OUTPUT を使用してデータセット country_corに出力します。
   b. ミャンマーとアメリカの相関関係を相関行列から抽出します。

  1. 1990年と2000年の死亡率情報間に差があるかどうかを、対応のある t 検定とウィルコクソン符号順位検定で調べてください。ヒント: 1990年の情報列を取り出すには ‘1990’n を使用します。

  2. データセット carsを使用して、車両価格「VehBCost」 を従属変数とし、車両年齢「VehicleAge」とオンライン販売 「IsOnlineSale」の有無とその相互作用を予測変数とする線形回帰モデルを使用します。

  3. データセットcars内に車両価格が 10,000 ドルを超えていることを示す変数「expensive」を作成してください。カイ二乗検定により車が高価だったことと、損な取引だったことを示すラベル用変数「IsBadbuy」に関連があるかを評価してください。

  4. ロジスティック回帰を使用し、 応答変数を「IsBadbuy」予測変数を「expensive」「VehicleAge」としてオッズ比の信頼区間を算出してください。

Contents