SAS_codes_multivariate_3(Lab2&hw3)

TEMP·2021년 10월 9일

SAS

목록 보기

9/13

tinv(0.90, 14)는 자유도 14인 t분포에서 -inf 부터 해서 0.9가 되는 값을 반환한다.
따라서 유의확률 5%에서 양측검정을 하고 싶으면 tinv(0.975, d.f)를 구해야 한다.

평균 검정(일변량)

data

%let path = D:\Chrome;
proc import datafile = "&path\sweat.csv" dbms = csv replace out = sweat;
getnames =  yes;
run;

proc print data = sweat;
run;

proc means data = sweat maxdec = 2;
var sweatrate;
run;

모평균 검정

모분산을 안다면 Z를 사용할 텐데 모르니까 모분산(모표준편차) 대신 표분분산(표본표준편차)으로 바꾸고 t 분포 사용하면 된다.
https://documentation.sas.com/doc/en/pgmsascdc/9.4_3.5/procstat/procstat_univariate_syntax01.htm
위 doc 참고함.

proc univariate data = sweat mu0=4;
var sweatrate;
run;
해당 변수에 대하여 모평균이 4인지 검정하는 방법이다.

일변수 식이야 어렵지 않으니 직접 t-value를 유도해보고 같은지도 확인해보자.
위의 proc means에서 나온 표본평균과 표본표분편차을 이용하자.
SAS가 이래서 싫다.
data t_value;
x_bar = 4.64;
std = 1.70;
n = 20;
mu_0 = 4;
t_value = ( sqrt(n) *(  x_bar - mu_0 ) )/ (std );
run;
해보니 결과가 같다. 잘했다.
proc ttest data = sweat h0 = 4;
var sweatrate;
run;
이렇게 해도 된다.
당연히 결과는 같은데 출력해주는게 쫌 다르다.
univariate는 기술통계량을 출력하고 test값을 추가적으로 출력해주는 느낌이고
ttest는 진짜 test procedure인거 같다. QQ plot하고 CL도 출력해준다.
찾아보니까 mu0는 말 그대로 모평균을 의미하고 h0는 귀무가설을 의미하는거 같다.

평균 검정(다변량)

data

proc means data = sweat maxdec = 2;
var sweatrate sodium potassium;
run;
또는
data meanQ1;
set sweat_corr;
if _TYPE = "MEAN";
keep x1-x3;
run;
각 3개의 변수에 대한 평균
proc corr data = sweat cov outp = sweat_corr;
var sweatrate sodium potassium;
run;
3개의 변수에 대한 공분산과 상관계수
data sweat_cov_mat;
set sweat_corr;
if _type_ = "COV";
drop _character_;
run;
공분산 dataset

모평균 검정

마찬가지로 모공분산행렬은 모른다. (근데 아는 경우가 있나? 평균을 모르는데?)
$\mu$ = (4, 50, 10)인지 검정 (동시에 만족하는지 이다.)

직접계산해보기

proc iml;
sample_mean = {4.64 45.40 9.965};
mu0 = {4 50 10};
n = 20;
p = 3;
use sweat_cov_mat;
read all into cov_mat;
mu_ = sample_mean - mu0;
mu_t = mu_`;
inv_cov = inv(cov_mat);
T2 = n * mu_* inv_cov * mu_t ;
q1 = ( (n-1) * p / (n - p) )* finv(.95, p, n-p);
q2 = ( (n-1) * p / (n - p) )* finv(.90, p, n-p);
print cov_mat, inv_cov, mu_ , T2, mu_t, q1, q2 ;
run;
proc mean을 이용하여 표본평균을 계산한다.
귀무가설을 설정한다.
표본의 개수와 변수의 개수를 입력한다.
공분산 dataset을 matrix로 만든다.
호텔링의 T_square 값을 계산한다.
다음 finv는 F 분포에서 해당하는 확률과 자유도에 대한 값을 반환한다.
이를 통하여 95%와 90%에 대한 기각역을 계산한다.
이때 순서대로 $\alpha$ 는 0.05와 0.1이다.

proc glm이용하여 계산하기

data sweat_h;
set sweat;
nx1 = sweatrate - 4;
nx2 = sodium - 50;
nx3 = potassium - 10;
run;

proc glm data = sweat_h;
model nx1 nx2 nx3 = / nouni ;
manova h = intercept;
run;

proc iml; 
T = 19 * (1-0.66112774)/0.66112774; 
print T;
run;

산점도

proc G3D data = sweat;
scatter sodium * sweatrate = potassium;
run;

두개의 모집단 모평균이 같은지 검정

직접계산

proc corr data = man cov outp = out_man; 
run;
proc corr data = woman cov outp = out_woman; 
run;

공분산 dataset을 저장한다.

data out_man_woman ;
set out_man out_woman; 
x1 = round(x1, 0.001);
x2 = round(x2, 0.001);
x3 = round(x3, 0.001);
x4 = round(x4, 0.001);
run;

한번에 계산해야 하므로 둘을 같이 저장하되 유효숫자를 정하여 반올림한다.

PROC IML;
p = 4;
USE out_man_woman;
read all into A;
S_man = A[1:4, 1:4];
S_woman = A[12:15, 1:4];
M_man = A[5, 1:4];
M_woman = A[16, 1:4];
N_man = A[7,3];
N_woman = A[18,3];
S_p = ( ( N_man-1 )*S_man + ( N_woman-1 )*S_woman )/( N_man + N_woman -2 );
T2 = ( N_man*N_woman / ( N_man + N_woman  ) )*( M_man - M_woman )*inv(S_p)*( M_man` - M_woman` );
constant_F =  ( ( N_man + N_woman -2 )*p)/( N_man + N_woman - p - 1 )*finv(.95, p,  N_man + N_woman - p-1);
print S_p, T2, constant_F, F ;
run;

HW

위 내용이긴 한데 연습하면서 code가 간결해짐.

1.

$H_0$ : $\mu$ = $(4, 48, 11)$
(a) 검정통계량 $T^2$ 을 구하라.
(b) 위에서 구한 $T^2$ 값을 이용하여 가설검정을 실시하라.

proc corr data = Q1 cov outp = co_Q1; run;
공분산과 상관계수 행렬을 구하고
data mean_Q1;
set co_Q1;
if _TYPE_ = "MEAN";
keep x1-x3;
run;
위의 dataset은 표본평균도 구해주므로 따로 구하지 않고 여기서 추출한다.
data cov_Q1;
set co_Q1;
if _TYPE_ = "COV";
keep x1-x3;
run;
공분산 dataset을 추출한다.
proc iml;
h0 = {4 48 11};
n = 20;
p = 3;
use mean_Q1;
read all into samplemeans;
use cov_Q1;
read all into S;
T2 = n*( samplemeans - h0 )*inv(S)*( samplemeans` - h0` );
q5 = ( (n-1) * p / (n - p) )* finv(.95, p, n-p);
q10 = ( (n-1) * p / (n - p) )* finv(.90, p, n-p);
print samplemeans S T2 q5, q10;
귀무가설, 전체 데이터 개수, 변수개수를 적는다.
위에서 만들어 놓은 표본평균 dataset을 matrix(row vector)로 바꾼다.
공분산 행렬을 만든다.
$T^2$ 를 계산한다.
다음 해당 유의확률별로 기각역의 임계값을 계산한다.
(원래는 $T^2$ 분포의 값을 계산해야 하는데 상수를 곱해서 $F$ 분포를 따름이 알려져 있음.)
$F$ 분포는 대칭이 아니니까 예를들어 유의확률 5%이면 $-\infty$ 부터하여 넓이가 0.95가 되는 값을 계산해주면 된다.

2.

a) 남자와 여자의 표본 평균벡터를 구하라.
(b) 남자와 여자의 표본 공분산행렬을 구하라.
(c) 남자와 여자의 표본 공통 공분산행렬을 구하라.
(d) 표본 공통 공분산 행렬을 이용하여 남자와 여자의 평균벡터가 같은지 검정하라.
(e) 각 변수별로 평균이 같은지 일변량 검정을 실시하라.

proc corr data = male cov outp = cov_male; run;
proc corr data = female cov outp = cov_female; run;
data cov; set cov_male cov_female; proc print; run;

proc iml;
p = 4;
USE cov;
read all into A;
S_man = A[1:4, 1:4];
S_woman = A[12:15, 1:4];
M_man = A[5, 1:4];
M_woman = A[16, 1:4];
N_man = A[7,3];
N_woman = A[18,3];
S_p = ( ( N_man-1 )*S_man + ( N_woman-1 )*S_woman )/( N_man + N_woman -2 );
T2 = ( N_man*N_woman / ( N_man + N_woman  ) )*( M_man - M_woman )*inv(S_p)*( M_man` - M_woman` );
constant_F5 =  ( ( N_man + N_woman -2 )*p)/( N_man + N_woman - p - 1 )*finv(.95, p,  N_man + N_woman - p-1);
constant_F10 =  ( ( N_man + N_woman -2 )*p)/( N_man + N_woman - p - 1 )*finv(.9, p,  N_man + N_woman - p-1);
print S_p, T2, constant_F5, constant_F10 ;
quit;

proc iml;
n = 28;
USE cov;
read all into A;
mean_vector_man = A[5, 1:4];
mean_vector_woman = A[16, 1:4];
var_male_y1 = A[1,1];
var_female_y1 = A[12,1];
var_male_y2 = A[2,2];
var_female_y2 = A[13,2];
var_male_y3 = A[3,3];
var_female_y3 = A[14,3];
var_male_y4 = A[4,4];
var_female_y4 = A[15,4];
S_p_y1 = ( ( n-1 )*var_male_y1 + ( n-1 )*var_female_y1 ) / ( n+n-2 );
S_p_y2 = ( ( n-1 )*var_male_y2 + ( n-1 )*var_female_y2 ) / ( n+n-2 );
S_p_y3 = ( ( n-1 )*var_male_y3 + ( n-1 )*var_female_y3 ) / ( n+n-2 );
S_p_y4 = ( ( n-1 )*var_male_y4 + ( n-1 )*var_female_y4 ) / ( n+n-2 );
t_y1 = ( mean_vector_man[1] - mean_vector_woman[1] ) / ( S_p_y1*sqrt( 1/n + 1/n ) );
t_y2 = ( mean_vector_man[2] - mean_vector_woman[2] ) / ( S_p_y2*sqrt( 1/n + 1/n ) );
t_y3 = ( mean_vector_man[3] - mean_vector_woman[3] ) / ( S_p_y3*sqrt( 1/n + 1/n ) );
t_y4 = ( mean_vector_man[4] - mean_vector_woman[4] ) / ( S_p_y4*sqrt( 1/n + 1/n ) );
t_value_5 = tinv(0.975, n+n-2);
t_value_10 = tinv(0.95, n+n-2);
print S_p_y1, S_p_y2, S_p_y3, S_p_y4, t_y1, t_y2, t_y3, t_y4, t_value_5, t_value_10;
quit;

3.

data Q3;
input y1 y2 x1 x2;
cards;
148 20 137 15
159 24 164 25
144 19 224 27
103 18 208 33
121 17 178 24
89 11 128 20
119 17 154 18
123 13 158 16
76 16 102 21
217 29 214 25
148 22 209 24
151 21 151 16
83 7 123 13
135 20 161 22
178 15 175 23
;
run;

data D;
set Q3;
d1 = x1 - y1;
d2 = x2- y2;
keep d1 d2;
run;

proc corr data = D cov outp = cov_D; run;

proc iml;
n = 15;
p = 2;
use cov_D;
read all into A;
samplemean = A[3,];
S = A[1:2,1:2];
T2 = n*samplemean*inv(S)*(samplemean`);
constant_F_value5 = finv(.95, p, n-p)*( n - 1 )*p / ( n - p );
constant_F_value10 = finv(.90, p, n-p)*( n - 1 )*p / ( n - p );
print T2, constant_F_value5, constant_F_value10;

proc iml;
n = 15;
use cov_D;
read all into A;
samplemean = A[3,];
S = A[1:2,1:2];
t_1 = samplemean[1] / sqrt( S[1,1]/n );
t_2 = samplemean[2] / sqrt( S[2,2]/n );
t_value_5 = tinv(0.975, n-1);
t_value_10 = tinv(0.95, n-1);
print t_1, t_2, t_value_5, t_value_10; quit;