黃新智隨手筆記: [教學示範]「使用工具：SAS」僅有單一變項，如何對此資料的數據分布加以歸類？

2019年10月8日星期二

[教學示範]「使用工具：SAS」僅有單一變項，如何對此資料的數據分布加以歸類？

先帶入2個專有名詞。
聚類分析：用在事先不知道類別的情況下，完全按照反映對象特徵的數據將對象進行分類。
判別分析：是基於已知分類的樣本建立判別函數，再對未知類別的大量個體歸屬於哪個類別進行判別。
所以聚類分析是一種沒有類別信息可參考的情況下，利用距離或相似性係數將一個集合劃分成若干個子集的過程，它是種監督式的學習過程。相關前述粉紅色內容，通常具有多個變項供軟體綜合判定，但也可能只有一個變項可供軟體判定，故接下來，本文章將聚焦在後者，如果你只有一個連續性的變項（即下例中的x），該如何分類？

＝＝＝以上是版大集中腦力寫出來的段落，然為節省體內能量，以下段落寫法會較隨性＝＝＝

版大先以2個迴圈產生2批虛擬資料，其描述性統計值設定如下：
data test;
grp=1; /*平均數0 標準差0.5，共有100筆觀察值，批數grp=1*/
do i=1 to 100;
x=0+0.5*rannor(1234);
output;
end;
grp=2; /*平均數3 標準差1，共有200筆觀察值,批數grp=2*/
do i=1 to 200;
x=3+1*rannor(1234);
output;
end;
run;

... ...

data test_2 test_3;
set test;
drop i;
if grp=1 then output test_2;
else output test_3;
run;

data test_3;
set test_3;
rename grp=grp_2 x=x_2;
run;

data test_4;
merge test_2 test_3;
run;

proc sgplot data=test_4;
density x / type=normal lineattrs=(color=red) legendlabel='grp=1';
density x_2 / type=normal lineattrs=(color=blue) legendlabel='grp=2';
histogram x / transparency=0.75 fillattrs=(color=red);
histogram x_2 / transparency=0.75 fillattrs=(color=blue);
keylegend / location=outside position=bottom;
xaxis label="Normal Curves" labelattrs=(size=15) valueattrs=(size=12);
yaxis grid;
run;

版大的虛擬資料設定的較誇張，這二批資料分的比較開，如果把這兩批資料合併成一個直方圖，將如下圖。

好！！！！從現在開始，假設你只知道這是一個連續變項x的直方圖，畫出的結果如下：

proc sgplot data=test;
density x / type=normal lineattrs=(color=red) legendlabel='grp=all';
histogram x / transparency=0.75 fillattrs=(color=green);
keylegend / location=outside position=bottom;
xaxis label="Normal Curves" labelattrs=(size=15) valueattrs=(size=12);
yaxis grid;
run;

請問它可能是由幾個群組(類別)的資料所混出？
各位！這個時候～～

Proc fmm data=test gconv=0;
model x=/ kmin=2 kmax=5; /*跟SAS說等一下你給我的分群數建議，要介於2～5組間，若膽敢超過，就打斷你的腿！*/
output out=abc_output class; /*分析後的結果，輸出至abc_output資料集，class是要求列出分群的組別，值得一提的是，class也可用category或group取代(SAS官方教學文件有提，見本文下方proc fmm進階設定)。另外，若你想要每筆觀察值對應的所有統計值，可以改成output out=abc_output / allstats */
run;

以下只列出重要統計資訊：

上圖：SAS認為上述的x可能是由3個群組所貢獻。當等於3時，AIC最小，所以此分群最適合。AIC常用於模型比較，在時間序列分析如ARIMA(差分整合移動平均自迴歸模型)中，AIC也常被拿出來說教XD

上圖：這3個群組，其平均數與變異數，還有其在混合模型中貢獻的佔比如下：

第1個群組：(0.07047, 0.3695, 0.3541)

第2個群組: （6.7314, 0, 0.0033）

第3個群組：(3.1484, 0.8506,0.6426)

分布和估計密度圖如下：

此圖先不解釋！

咦，版大不是用2個群組混進來嗎？SAS怎會推薦3組？來，我們換用另一支程式畫圖，如下：

ods graphics / imagefmt=png width=1200px height=900px ;
proc univariate data=test noprint;
var x;
HISTOGRAM / NORMAL kernel(c= 0.8 ) midpoints=-2 to 6.5 by 0.4 barlabel=count vaxislabel="count";
inset n="資料數" median="中位數" mean="平均數" std="標準差" / POSITION = ne ;
run;
ods graphics off;

終於抓到你了，就是最右邊那支bar, 就是上面講的第2個群組。
它只出現1次，我們可以忽略它(SAS是死的，人才是活的)，所以程式重改再跑。

ods graphics / imagefmt=png width=1200px height=900px ;
proc fmm data=test gconv=0;
model x=/ kmin=2 kmax=2; /*控制在=2*/
output out=abc_output class;
run;
ods graphics off;

結果如下：

再來複習一下：

原混入的群組1為
平均數0 標準差0.5，共有100筆觀察值，批數grp=1
0.5的平方等於0.25 ，已很接近0.3492這個變異數，而0也很接近0.05777這個平均數
此群組在混合模型中佔了0.3492，在實際值中，它是佔了100/300，很接近的值

原混入的群組2為
平均數3 標準差1，共有200筆觀察值,批數grp=2
1的平方還是1，很接近0.9428這個變異數，而3也很接近3.1505這個平均數
此群組在混合模型中佔了0.6508，在實際值中，它是佔了200/300，很接近的值

畫在圖上的視覺就是這樣。

紅線群組，其平均數為0.06，變異數為0.36

綠線群組，其平均數為3.15，變異數為0.94

至於0.3492旁，有個-0.6226的數字，這個道理與Logistic regression中的Logit(p)轉機率一樣，須透過一個公式。但公式不一樣！。這裡的公式是：

匯出的資料集abc_output如下(侷部擷圖)：
GRP是真實的分群，CLASS是SAS推估建議的分群

意義：
當你多了class變項(在實際操作面上，常常不存在grp這樣的變項)後，就能再跑些原來無法跑的統計(例如one-way ANOVA，或Chi-square test)，也許就因為多了這一步而讓你發現新資訊。

最後，以交叉表做結。
proc freq data=abc_output;
tables grp*class/ chisq;
run;