통컨(통계컨설팅)

[SAS]

DATA iris;INFILE '/home/joinos0/sas_class/iris.csv' DLM=',' FIRSTOBS=2;
LENGTH Species $15;
INPUT sepal_length sepal_width petal_length petal_width species $;
PROC PRINT;
RUN;

DATA iris1;SET iris;
IF _N_ <= 10;
PROC CLUSTER DATA=iris1 METHOD=SINGLE;                     
PROC TREE HORIZONTAL SPACES=2;
RUN;

[R]

head(iris)
dist(iris[1:10,1:4])                 # 편의상 7개 관측치 간의 거리를 구한다...
round(dist(iris[1:10,1:4]),digits=4) # 소수7자리 -> 소수 4 자리
dist01 <- round(dist(iris[1:10,1:4]),digits=4);dist01 # 거리를 객체 dis01 저장 
model_hc <- hclust(dist01, method="ave")
plot(model_hc)
plot(model_hc,hang=-1)

[Python]

from sklearn.datasets import load_iris
iris = load_iris()
# 계층적 군집 model
from scipy.cluster.hierarchy import linkage, dendrogram
clusters = linkage(y=iris['data'], method='complete', metric='euclidean')
clusters
clusters.shape # (149, 4)
import matplotlib.pyplot as plt
plt.figure( figsize = (25, 10) )
dendrogram(clusters, leaf_rotation=90, leaf_font_size=12,)
# leaf_rotation=90 : 글자 각도
# leaf_font_size=20 : 글자 사이즈
plt.show() 

[SAS]

DATA iris;INFILE '/home/joinos0/sas_class/iris.csv' DLM=',' FIRSTOBS=2;
LENGTH Species $15;
INPUT sepal_length sepal_width petal_length petal_width species $;
DATA iris1;SET iris;
IF species ^= "Iris-setosa  ";  * 붓꽃 중에서 setosa를 제외, versicolor, virginica 두 종류만 추출
PROC CANDISC DATA=iris1;      
CLASS species;                          
VAR sepal_length sepal_width petal_length petal_width;                      
RUN;

SAS는 분석결과물이 상당히 많습니다. 이중에서 R 결과와 비교할 만한 부분을 추출하였습니다.

R 결과와 동일하게 나오는 것을 볼 수 있습니다.

[R]

# R 이용한 다변량분석 
library(MASS)
iris1 <- iris   #  iris 데이터를 iris1 으로 복사
head(iris1)    # iris1 데이터 구성 살펴보기
str(iris1)

# (1) setosa 를 제외하기
levels(iris1$Species)
iris1 <- subset(iris1,Species=="versicolor" | Species=="virginica")
levels(iris1$Species)
iris1$Species <- factor(iris1$Species)
model_lda <- lda(Species~.,data=iris1)  # lda 를 이용하여 판별분석
model_lda

[Python]

. 보완예정

from sklearn.discriminant_analysis import LinearDiscriminantAnalysis
from sklearn.preprocessing import StandardScaler
from sklearn.datasets import load_iris
import pandas as pd
import matplotlib.pyplot as plt
%matplotlib inline

# 붓꽃 데이터 로드
iris = load_iris()

# 데이터 정규 스케일링
iris_scaled = StandardScaler().fit_transform(iris.data)

# 2개의 클래스로 구분하기 위한 LDA 생성
lda = LinearDiscriminantAnalysis(n_components=2)

# fit()호출 시 target값 입력 
lda.fit(iris_scaled, iris.target)
iris_lda = lda.transform(iris_scaled)

[SAS]

SAS에서 요인분석을 하려면 PROC FACTOR을 사용합니다.

DATA a1;INPUT x1 x2 ;CARDS;
 4  15
 6  16   
 7  11
 8  10  
 9   6
 11  8
 12 10
 13 14   
;
PROC FACTOR ;VAR x1 x2;RUN;      /* 가장 기본적인 요인분석 실행*/          
PROC FACTOR MEHTOD=PRINIPAL;VAR x1 x2;RUN; /* 디폴트로 principal */

PROC FACTOR NFACTORS=2 ;VAR x1 x2;RUN;  /* 요인 갯수 지정  */
PROC FACTOR MINEIGEN=0 ;VAR x1 x2;RUN;  /* 고유값의 하한값 지정  */
PROC FACTOR NFACTORS=2 EIGENVECTORS;VAR x1-x2; RUN; /* 고유벡터를 구함 */
PROC FACTOR NFACTORS=2 EIGENVECTORS SCORES;VAR x1-x2; RUN; /* 요인점수를 구함 */
RUN; 
PROC FACTOR MEHTOD=ML;VAR x1 x2;RUN; /* 요인분석 초기치 지정  */
PROC FACTOR NFACTORS=2 ROTATE=VARIMAX;VAR x1-x2;RUN; /* 요인의 회전 */      
PROC FACTOR NFACTORS=2 OUT=out1 OUTSTAT=out2; VAR x1-x2;RUN;/* 요인분석 결과를 데이터로 저장 */

PROC PRINT DATA=out1;RUN;
PROC PRINT DATA=out2;RUN;

[R]

R에서 요인분석을 하려면 factoranal 을 사용합니다. 변수가 2개인 경우는 에러가 발생합니다.

패키지 factanal 이용
x1 <- c( 4,  6,  7,  8, 9, 11, 12, 13);
x2 <- c(15, 16, 11, 10, 6,  8, 10, 14);
(dataf01 <- data.frame(x1,x2))
factanal(dataf01,factors=2) 
factanal(dataf01,factors=2,rotation="none") 

패키지 psych 를 이용

install.packages("psych")
# package ‘tmvnsim’ successfully unpacked and MD5 sums checked
# package ‘mnormt’ successfully unpacked and MD5 sums checked
# package ‘psych’ successfully unpacked and MD5 sums checked
library(psych)
principal(dataf01,2,scores=TRUE)

[Python]

import pandas as pd
df = pd.read_csv("d:/sas_class/iris.csv")  # iris 데이터 불러오기
df
df.head()
X=df.drop(['species'],axis=1)
X.head()
# pip install factor_analyzer : factror_analyzer 가 설치되지 않은 경우, 설치함
from factor_analyzer import FactorAnalyzer
fa=FactorAnalyzer(method="principal", n_factors=3,rotation="varimax").fit(X) # varimax로 회전
pd.DataFrame(fa.loadings_,index=X.columns)
ev, v=fa.get_eigenvalues()
ev

fa.transform(X)

(1-1) SAS 프로그램 결과(상관행렬 이용)

PROC PRINCOMP ;VAR x1 x2;RUN; 의 결과 (default로 상관행렬을 이용함)

아무런 옵션이 없이 그냥 PROC PRINCOMP를 이용하면 상관행렬을 이용합니다.

주성분1 Prin1 = 0.707107 X x1 - 0.707107 X x2

주성분2 Prin2 = 0.707107 X x1 + 0.707107 X x2 로 구할 수 있습니다.

주성분1이 설명하는 변량은 1.40148이고

주성분2가 설명하는 변량은 0.59852입니다.

다음은 주성분갯수를 선택하는데 도움이 되는 Scree 플롯인데 변수가 2개이라서 별로 의미가 없어 보이지만

변수가 여러 개인 경우에는 주성분의 갯수를 결정하는 데 도움이 됩니다.

(1-2) R 프로그램 결과(상관행렬 이용)

prin_cor <-princomp(dataf01,cor=T) # 상관계수 이용, prin_cor 객체 생성

summary(prin_cor)  # 객체 통계량을 보면 주성분의 표준편차와 설명력을 출력

prin_cor$sdev ^2   #3 표준폍차를 제곱해 보면 주성분의 설명력 출력

prin_cor$loadings  # loadings 를 출력해 보면 고유벡터를 확인해 볼 수 있습니다. 

주성분1 Prin1 = 0.707107 X x1 - 0.707107 X x2

주성분2 Prin2 = 0.707107 X x1 + 0.707107 X x2 로 구할 수 있습니다.

주성분1이 설명하는 변량은 1.40148이고

주성분2가 설명하는 변량은 0.59852입니다.

* 참고 R 을 이용한 고유치와 고유벡터의 결과

상관행렬을 구한 다음 eigen( ) 함수를 이용하면 고유벡터와 고유치를 확인할 수 있습니다.

> eigen(cor(dataf01))

SAS와 비교하면 주성분계수의 부호가 음과 양으로 바뀌어 있고

주성분의 고유치(eigen value, 설명하는 변량)는 동일한 것을 볼 수 있습니다.

(2-1) SAS 프로그램 결과(공분산행렬 이용)

PROC PRINCOMP COV;VAR x1 x2;RUN; 의 결과 (COV 옵션에 의해) 공분산행렬을 이용함

공분산을 사용하려면 PROC PRINCOMP에 cov 옵션을 이용합니다.

주성분1 Prin1 = -0.598741 X x1 + 0.800942 X x2

주성분2 Prin2 =  0.800942 X x1 + 0.598741 X x2 로 구할 수 있습니다.

주성분1이 설명하는 변량은 15.47145이고

주성분2가 설명하는 변량은 6.38569입니다.

(2-2) R 프로그램 결과(공분산행렬 이용)

* 참고 R 의 결과 - SAS와 같은 결과

x1 <- c( 4, 6, 7, 8, 9, 11, 12, 13);

x2 <- c(15,16,11,10, 6, 8, 10, 14);

(dataf01 <- data.frame(x1,x2));

eigen(cov(dataf01))

성분1 Prin1 = -0.598741 X x1 + 0.800942 X x2

주성분2 Prin2 =  0.800942 X x1 - 0.598741 X x2 로 구할 수 있습니다.

주성분1이 설명하는 변량은 15.47145이고

주성분2가 설명하는 변량은 6.38569입니다.

AS와 비교하면 주성분계수의 부호가 음과 양으로 바뀌어 있고

주성분의 고유치(eigen value, 설명하는 변량)는 동일한 것을 볼 수 있습니다.

주성분분석에서 상관행렬(또는 공분산행렬)이 있을 때에는 covmat 옵션을 사용하면 됩니다.

?princomp

?prcomp 를 이용하여 도움말을 보면 covmat 옵션은 princomp 에만 가능한 것을 확인했습니다,

> x1 <- c( 4, 6, 7, 8, 9, 11, 12, 13);

> x2 <- c(15,16,11,10, 6, 8, 10, 14);

> (dataf01 <- data.frame(x1,x2));

> mat_cov <- cov(dataf01)

> princomp(covmat=mat_cov)

> x1 <- c( 4, 6, 7, 8, 9, 11, 12, 13);

> x2 <- c(15,16,11,10, 6, 8, 10, 14);

> (dataf01 <- data.frame(x1,x2));

> mat_cor <- cor(dataf01)

> princomp(covmat=mat_cor)

x1 <- c( 4, 6, 7, 8, 9, 11, 12, 13);

x2 <- c(15,16,11,10, 6, 8, 10, 14);

(dataf01 <- data.frame(x1,x2));

out <-princomp(dataf01,cor=F,scores=T)

names(out)

out$score

out$loadings

summary(out)

(1) dataf01 <- data.frame(x1,x2)  변수 x1, x2를 데이터프레임으로 묶습니다.

(2) princomp(dataf01,cor=F,scores=T)의 객체 out를 만듭니다.

주성분분석으로 princomp 함수를 사용합니다.
첫 인자 dataf01 은 분석에 사용되는 데이터이고

cor=F 상관계수행렬을 사용하지 않고 공분산행렬을 사용하겠다는 의미

scores=T 는 주성분점수를 계산하겠다는 의미

names(out) 를 보면 score 라는 변수가 있는 것을 확인할 수 있습니다.

(3) 주성분점수를 출력합니다. 객체 out의 score 객체에 저장되어 있습니다.

(4) 주성분점수를 구하는 계수를 구합니다.

comp1= 0.599  X x1 - 0.801  X x2

comp2 = 0.801 X x1 + 0.599 X x2

(5) 주성분 1과 2가 설명하는 설명력을 출력합니다.

summary( )함수를 사용하면 주성분들이 각 변수의 변량을 얼마나 설명하는지 비율을 구할 수 있습니다.

제1주성분이 70.78%를 설명하고, 제2주성분이 29.21%를 설명하는 것을 보이고 있습니다.

*주성분분석에 대해서 궁금하신 부분이 있으시면 댓글에 글을 남겨 주시면

제가 아는 범위내에서 설명드리도록 하겠습니다.