Similarity Calculation 5 - Mutual Information Using SQL and SAS

I(X,Y) = H(X,Y) - H(X|Y) - H(Y|X) 
Intuitively, we can interpret the MI between X and Y as the reduction in uncertainty about X about observing Y, or by symmetry, the reduction in uncertainty about Y after observing X.

-- Calculate population counts/prob;

-- 4,856,590 33,853,878 38,710,468;
select sum(response), count(1)-sum(response), count(1)
from ANALYTICS_STG..LH_CAMPAIGN_DM_ADM1;


-- 1 Calculate positive marginal counts for x=1;
drop table ANALYTICS_STG..LH_CAMPAIGN_DM_ADM1_pos;
-- 1 rows affected;
create table ANALYTICS_STG..LH_CAMPAIGN_DM_ADM1_pos as
select
sum(UNIV_BBBUS)UNIV_BBBUS,
sum(UNIV_BABYUS)UNIV_BABYUS,
sum(UNIV_CTSUS)UNIV_CTSUS,
sum(UNIV_HRMUS)UNIV_HRMUS,
sum(UNIV_BBBCA)UNIV_BBBCA,
sum(UNIV_BBBMX)UNIV_BBBMX,
sum(UNIV_UNIV_BABYCA)UNIV_UNIV_BABYCA,
...
from ANALYTICS_STG..LH_CAMPAIGN_DM_ADM1
;

-- SAS
proc sql;
create table pos as
select * from connection to NETEZZA
(
select * from ANALYTICS_STG..LH_CAMPAIGN_DM_ADM1_pos
;
)
;
quit;

proc means data=pos noprint;
output out=meanout(drop=_type_ _freq_ where=(_stat_ in ('MIN')));
run;

proc transpose data=meanout out=meanout;
run;

* varlist;
proc sql;
select _name_
from meanout where col1>0.0;
quit;


-- 2 Calculate entropy for marginal attributes;
drop table ANALYTICS_STG..LH_CAMPAIGN_DM_ADM2_pos;
-- 1 rows affected;
create table ANALYTICS_STG..LH_CAMPAIGN_DM_ADM2_pos as
select 
-(UNIV_BBBUS/38710468)*ln(UNIV_BBBUS/38710468)-(1-UNIV_BBBUS/38710468)*ln(1-UNIV_BBBUS/38710468)UNIV_BBBUS,
-(UNIV_BABYUS/38710468)*ln(UNIV_BABYUS/38710468)-(1-UNIV_BABYUS/38710468)*ln(1-UNIV_BABYUS/38710468)UNIV_BABYUS,
-(UNIV_CTSUS/38710468)*ln(UNIV_CTSUS/38710468)-(1-UNIV_CTSUS/38710468)*ln(1-UNIV_CTSUS/38710468)UNIV_CTSUS,
-(UNIV_HRMUS/38710468)*ln(UNIV_HRMUS/38710468)-(1-UNIV_HRMUS/38710468)*ln(1-UNIV_HRMUS/38710468)UNIV_HRMUS,
-(UNIV_BBBCA/38710468)*ln(UNIV_BBBCA/38710468)-(1-UNIV_BBBCA/38710468)*ln(1-UNIV_BBBCA/38710468)UNIV_BBBCA,
-(UNIV_BBBMX/38710468)*ln(UNIV_BBBMX/38710468)-(1-UNIV_BBBMX/38710468)*ln(1-UNIV_BBBMX/38710468)UNIV_BBBMX,
...
from ANALYTICS_STG..LH_CAMPAIGN_DM_ADM1_pos a
;


drop table ANALYTICS_STG..LH_CAMPAIGN_DM_ADM3_pos;
-- 3 Caculate conditional prob for y=1;
create table ANALYTICS_STG..LH_CAMPAIGN_DM_ADM3_pos as
select 
sum(UNIV_BBBUS)UNIV_BBBUS,
sum(UNIV_BABYUS)UNIV_BABYUS,
sum(UNIV_CTSUS)UNIV_CTSUS,
sum(UNIV_HRMUS)UNIV_HRMUS,
sum(UNIV_BBBCA)UNIV_BBBCA,
sum(UNIV_BBBMX)UNIV_BBBMX,
sum(UNIV_UNIV_BABYCA)UNIV_UNIV_BABYCA,
...
from ANALYTICS_STG..LH_CAMPAIGN_DM_ADM1
where response = 1
;

-- SAS
proc sql;
create table pos as
select * from connection to NETEZZA
(
select * from ANALYTICS_STG..LH_CAMPAIGN_DM_ADM3_pos
;
)
;
quit;

proc means data=pos noprint;
output out=meanout(drop=_type_ _freq_ where=(_stat_ in ('MIN')));
run;

proc transpose data=meanout out=meanout;
run;

* varlist;
proc sql;
select _name_
from meanout where col1>1.0;
quit;


drop table ANALYTICS_STG..LH_CAMPAIGN_DM_ADM4_pos;
-- 4 Caculate conditional prob for y=1;
create table ANALYTICS_STG..LH_CAMPAIGN_DM_ADM4_pos as
select 
-(POSTAL_ZIP4_85205_7911/4856590)*ln(POSTAL_ZIP4_85205_7911/4856590)-(1-POSTAL_ZIP4_85205_7911/4856590)*ln(1-POSTAL_ZIP4_85205_7911/4856590)POSTAL_ZIP4_85205_7911,
-(POSTAL_ZIP4_90403_5704/4856590)*ln(POSTAL_ZIP4_90403_5704/4856590)-(1-POSTAL_ZIP4_90403_5704/4856590)*ln(1-POSTAL_ZIP4_90403_5704/4856590)POSTAL_ZIP4_90403_5704,
...
from ANALYTICS_STG..LH_CAMPAIGN_DM_ADM3_pos;


drop table ANALYTICS_STG..LH_CAMPAIGN_DM_ADM5_pos;
-- 5 Caculate conditional prob for y=0;
create table ANALYTICS_STG..LH_CAMPAIGN_DM_ADM5_pos as
select 
sum(UNIV_BBBUS)UNIV_BBBUS,
sum(UNIV_BABYUS)UNIV_BABYUS,
sum(UNIV_CTSUS)UNIV_CTSUS,
sum(UNIV_HRMUS)UNIV_HRMUS,
sum(UNIV_BBBCA)UNIV_BBBCA,
sum(UNIV_BBBMX)UNIV_BBBMX,
sum(UNIV_UNIV_BABYCA)UNIV_UNIV_BABYCA,
...
from ANALYTICS_STG..LH_CAMPAIGN_DM_ADM1
where response = 0
;

-- SAS
proc sql;
create table pos as
select * from connection to NETEZZA
(
select * from ANALYTICS_STG..LH_CAMPAIGN_DM_ADM5_pos
;
)
;
quit;

proc means data=pos noprint;
output out=meanout(drop=_type_ _freq_ where=(_stat_ in ('MIN')));
run;

proc transpose data=meanout out=meanout;
run;

proc sort data=meanout;
by col1;
run;

* varlist;
proc sql;
select _name_
from meanout where col1>1 and col1<33853878;
quit;

drop table ANALYTICS_STG..LH_CAMPAIGN_DM_ADM6_pos;
-- 6 Caculate conditional prob for y=0;
create table ANALYTICS_STG..LH_CAMPAIGN_DM_ADM6_pos as
select 
-(POSTAL_ZIP4_11004_1040/33853878)*ln(POSTAL_ZIP4_11004_1040/33853878)-(1-POSTAL_ZIP4_11004_1040/33853878)*ln(1-POSTAL_ZIP4_11004_1040/33853878)POSTAL_ZIP4_11004_1040,
-(POSTAL_ZIP4_11364_3015/33853878)*ln(POSTAL_ZIP4_11364_3015/33853878)-(1-POSTAL_ZIP4_11364_3015/33853878)*ln(1-POSTAL_ZIP4_11364_3015/33853878)POSTAL_ZIP4_11364_3015,
...
from ANALYTICS_STG..LH_CAMPAIGN_DM_ADM5_pos;



drop table ANALYTICS_STG..LH_CAMPAIGN_DM_ADM_Entropy;
-- H(X) - P(Y=1)H(X|Y=1) - P(Y=0)H(X|Y=0)
-- 38,710,468 rows affected
create table ANALYTICS_STG..LH_CAMPAIGN_DM_ADM_Entropy as
select 
a.POSTAL_ZIP4_85205_7911-(4856590.0/38710468.0)*b.POSTAL_ZIP4_85205_7911-(33853878.0/38710468.0)*c.POSTAL_ZIP4_85205_7911 POSTAL_ZIP4_85205_7911,
a.POSTAL_ZIP4_90403_5704-(4856590.0/38710468.0)*b.POSTAL_ZIP4_90403_5704-(33853878.0/38710468.0)*c.POSTAL_ZIP4_90403_5704 POSTAL_ZIP4_90403_5704,
...
from ANALYTICS_STG..LH_CAMPAIGN_DM_ADM2_pos a
cross join ANALYTICS_STG..LH_CAMPAIGN_DM_ADM4_pos b
cross join ANALYTICS_STG..LH_CAMPAIGN_DM_ADM6_pos c;
;


-- Apply Mutual Information
drop table ANALYTICS_STG..LH_CAMPAIGN_DM_ADM2_Entropy;
-- 38,710,468 rows affected
create table ANALYTICS_STG..LH_CAMPAIGN_DM_ADM2_Entropy as
select a.COUPON_BARCODE, a.RESPONSE,
a.POSTAL_ZIP4_85205_7911*b.POSTAL_ZIP4_85205_7911+
a.POSTAL_ZIP4_90403_5704*b.POSTAL_ZIP4_90403_5704+
a.POSTAL_ZIP4_92173_3150*b.POSTAL_ZIP4_92173_3150+
a.POSTAL_ZIP4_90631_1103*b.POSTAL_ZIP4_90631_1103+
...
from ANALYTICS_STG..LH_CAMPAIGN_DM_ADM1 a
cross join ANALYTICS_STG..LH_CAMPAIGN_DM_ADM_Entropy b
;

select a.DECILE, count(*), sum(response)
from
(
select coupon_barcode, response, score,
NTILE(10) OVER(ORDER BY score DESC NULLS LAST) as DECILE
from ANALYTICS_STG..LH_CAMPAIGN_DM_ADM2_Entropy
) a
group by 1
order by 1 
;

Similarity Calculation 4 - Naive Bayes Using SQL and SAS

-- Calculate population counts/prob;
-- 4,856,590 33,853,878 38,710,468;
select sum(response), count(1)-sum(response), count(1)
from ANALYTICS_STG..LH_CAMPAIGN_DM_ADM1;


-- Calculate counts/prob for positive;
drop table ANALYTICS_STG..LH_CAMPAIGN_DM_ADM1_pos;
-- 1 rows affected;
create table ANALYTICS_STG..LH_CAMPAIGN_DM_ADM1_pos as
select
sum(UNIV_BBBUS)UNIV_BBBUS,
sum(UNIV_BABYUS)UNIV_BABYUS,
sum(UNIV_CTSUS)UNIV_CTSUS,
sum(UNIV_HRMUS)UNIV_HRMUS,
sum(UNIV_BBBCA)UNIV_BBBCA,
sum(UNIV_BBBMX)UNIV_BBBMX,
sum(UNIV_UNIV_BABYCA)UNIV_UNIV_BABYCA,
sum(DM_RECENTOPT_OPT_BABYUS)DM_RECENTOPT_OPT_BABYUS,
sum(DM_RECENTOPT_OPT_BBBUS)DM_RECENTOPT_OPT_BBBUS,
sum(DM_RECENTOPT_OPT_CTSUS)DM_RECENTOPT_OPT_CTSUS,
...
from ANALYTICS_STG..LH_CAMPAIGN_DM_ADM1
where response = 1;


-- Caculate counts/prob for negative;
drop table ANALYTICS_STG..LH_CAMPAIGN_DM_ADM1_neg;
-- 1 rows affected;
create table ANALYTICS_STG..LH_CAMPAIGN_DM_ADM1_neg as
select
sum(UNIV_BBBUS)UNIV_BBBUS,
sum(UNIV_BABYUS)UNIV_BABYUS,
sum(UNIV_CTSUS)UNIV_CTSUS,
sum(UNIV_HRMUS)UNIV_HRMUS,
sum(UNIV_BBBCA)UNIV_BBBCA,
sum(UNIV_BBBMX)UNIV_BBBMX,
sum(UNIV_UNIV_BABYCA)UNIV_UNIV_BABYCA,
sum(DM_RECENTOPT_OPT_BABYUS)DM_RECENTOPT_OPT_BABYUS,
sum(DM_RECENTOPT_OPT_BBBUS)DM_RECENTOPT_OPT_BBBUS,
sum(DM_RECENTOPT_OPT_CTSUS)DM_RECENTOPT_OPT_CTSUS,
...
from ANALYTICS_STG..LH_CAMPAIGN_DM_ADM1
where response = 0;


-- Drop a probability of zero which causes divide by zero;
proc sql;
create table pos as
select * from connection to NETEZZA
(
select * from ANALYTICS_STG..LH_CAMPAIGN_DM_ADM1_pos
;
)
;
quit;

proc means data=pos noprint;
output out=meanout(drop=_type_ _freq_ where=(_stat_ in ('MIN')));
run;

proc transpose data=meanout out=meanout;
run;

* varlist;
proc sql noprint;
select _name_ into :varlist separated by ' '
from meanout where col1>0.0;
quit;
%put &varlist;


/*Calculate probabilities for negative*/
proc sql;
create table neg as
select * from connection to NETEZZA
(
select * from ANALYTICS_STG..LH_CAMPAIGN_DM_ADM1_neg
;
)
;
quit;

proc means data=neg noprint;
var &varlist;
output out=meanout(drop=_type_ _freq_ where=(_stat_ in ('MIN')));
run;

proc transpose data=meanout out=meanout;
run;

* varlist;
proc sql;
select _name_
from meanout where col1>0.0;
quit;


-- Caculate conditional prob ratio for positive;
drop table ANALYTICS_STG..LH_CAMPAIGN_DM_ADM2_pos;
-- 1 rows affected;
create table ANALYTICS_STG..LH_CAMPAIGN_DM_ADM2_pos as
select
(a.UNIV_BBBUS/4856590.0)/(b.UNIV_BBBUS/33853878.0) UNIV_BBBUS,
(a.UNIV_BABYUS/4856590.0)/(b.UNIV_BABYUS/33853878.0) UNIV_BABYUS,
(a.UNIV_CTSUS/4856590.0)/(b.UNIV_CTSUS/33853878.0) UNIV_CTSUS,
(a.UNIV_HRMUS/4856590.0)/(b.UNIV_HRMUS/33853878.0) UNIV_HRMUS,
(a.UNIV_BBBCA/4856590.0)/(b.UNIV_BBBCA/33853878.0) UNIV_BBBCA,
(a.UNIV_BBBMX/4856590.0)/(b.UNIV_BBBMX/33853878.0) UNIV_BBBMX,
(a.UNIV_UNIV_BABYCA/4856590.0)/(b.UNIV_UNIV_BABYCA/33853878.0) UNIV_UNIV_BABYCA,
(a.DM_RECENTOPT_OPT_BABYUS/4856590.0)/(b.DM_RECENTOPT_OPT_BABYUS/33853878.0) DM_RECENTOPT_OPT_BABYUS,
...
from ANALYTICS_STG..LH_CAMPAIGN_DM_ADM1_pos a
cross join ANALYTICS_STG..LH_CAMPAIGN_DM_ADM1_neg b;


drop table ANALYTICS_STG..LH_CAMPAIGN_DM_ADM2_neg;
-- Caculate conditional prob ratio for negative;
create table ANALYTICS_STG..LH_CAMPAIGN_DM_ADM2_neg as
select
(1-a.UNIV_BBBUS/4856590.0)/(1-b.UNIV_BBBUS/33853878.0) UNIV_BBBUS,
(1-a.UNIV_BABYUS/4856590.0)/(1-b.UNIV_BABYUS/33853878.0) UNIV_BABYUS,
(1-a.UNIV_CTSUS/4856590.0)/(1-b.UNIV_CTSUS/33853878.0) UNIV_CTSUS,
(1-a.UNIV_HRMUS/4856590.0)/(1-b.UNIV_HRMUS/33853878.0) UNIV_HRMUS,
(1-a.UNIV_BBBCA/4856590.0)/(1-b.UNIV_BBBCA/33853878.0) UNIV_BBBCA,
(1-a.UNIV_BBBMX/4856590.0)/(1-b.UNIV_BBBMX/33853878.0) UNIV_BBBMX,
(1-a.UNIV_UNIV_BABYCA/4856590.0)/(1-b.UNIV_UNIV_BABYCA/33853878.0) UNIV_UNIV_BABYCA,
(1-a.DM_RECENTOPT_OPT_BABYUS/4856590.0)/(1-b.DM_RECENTOPT_OPT_BABYUS/33853878.0) DM_RECENTOPT_OPT_BABYUS,
...
from ANALYTICS_STG..LH_CAMPAIGN_DM_ADM1_pos a
cross join ANALYTICS_STG..LH_CAMPAIGN_DM_ADM1_neg b;


create table ANALYTICS_STG..LH_CAMPAIGN_DM_ADM_NB as
select a.*,
a.UNIV_BBBUS*b.UNIV_BBBUS+(1-a.UNIV_BBBUS)*c.UNIV_BBBUS+
a.UNIV_BABYUS*b.UNIV_BABYUS+(1-a.UNIV_BABYUS)*c.UNIV_BABYUS+
a.UNIV_CTSUS*b.UNIV_CTSUS+(1-a.UNIV_CTSUS)*c.UNIV_CTSUS+
a.UNIV_HRMUS*b.UNIV_HRMUS+(1-a.UNIV_HRMUS)*c.UNIV_HRMUS+
a.UNIV_BBBCA*b.UNIV_BBBCA+(1-a.UNIV_BBBCA)*c.UNIV_BBBCA+
a.UNIV_BBBMX*b.UNIV_BBBMX+(1-a.UNIV_BBBMX)*c.UNIV_BBBMX+
a.UNIV_UNIV_BABYCA*b.UNIV_UNIV_BABYCA+(1-a.UNIV_UNIV_BABYCA)*c.UNIV_UNIV_BABYCA+
...
a.FST_STR_SHOP_NBR_651*b.FST_STR_SHOP_NBR_651+(1-a.FST_STR_SHOP_NBR_651)*c.FST_STR_SHOP_NBR_651 score
from ANALYTICS_STG..LH_CAMPAIGN_DM_ADM1 a
cross join ANALYTICS_STG..LH_CAMPAIGN_DM_ADM2_pos b
cross join ANALYTICS_STG..LH_CAMPAIGN_DM_ADM2_neg c
;

select a.DECILE, count(*), sum(response)
from
(
select coupon_barcode, response, score,
NTILE(10) OVER(ORDER BY score DESC NULLS LAST) as DECILE
from ANALYTICS_STG..LH_CAMPAIGN_DM_ADM_NB
) a
group by 1
order by 1
;

Similarity Calculation 2 - Cosine Similarity Using SQL and SAS

-- Calculate population counts/prob;
-- 4,856,590 33,853,878 38,710,468;
select sum(response), count(1)-sum(response), count(1)
from ANALYTICS_STG..LH_CAMPAIGN_DM_ADM1;

-- 1 Calculate positive marginal counts for x=1;
drop table ANALYTICS_STG..LH_CAMPAIGN_DM_ADM1_pos;
-- 1 rows affected;
create table ANALYTICS_STG..LH_CAMPAIGN_DM_ADM1_pos as
select
sum(UNIV_BBBUS)UNIV_BBBUS,
sum(UNIV_BABYUS)UNIV_BABYUS,
sum(UNIV_CTSUS)UNIV_CTSUS,
sum(UNIV_HRMUS)UNIV_HRMUS,
sum(UNIV_BBBCA)UNIV_BBBCA,
sum(UNIV_BBBMX)UNIV_BBBMX,
...

sum(FST_STR_SHOP_NBR_651)FST_STR_SHOP_NBR_651
from ANALYTICS_STG..LH_CAMPAIGN_DM_ADM1
;

-- SAS
/*Calculate probabilities for positive*/
proc sql;
create table pos as
select * from connection to NETEZZA
(
select * from ANALYTICS_STG..LH_CAMPAIGN_DM_ADM1_pos
;
)
;
quit;

proc means data=pos noprint;
output out=meanout(drop=_type_ _freq_ where=(_stat_ in ('MIN')));
run;

proc transpose data=meanout out=meanout;
run;

* varlist;
proc sql;
select _name_
from meanout where col1>1.0;
quit;
%put &varlist;

-- 2 Calculate positive marginal counts given x=1 and y=1;
drop table ANALYTICS_STG..LH_CAMPAIGN_DM_ADM2_pos;
-- 1 rows affected;
create table ANALYTICS_STG..LH_CAMPAIGN_DM_ADM2_pos as
select
sum(UNIV_BBBUS)UNIV_BBBUS,
sum(UNIV_BABYUS)UNIV_BABYUS,
sum(UNIV_CTSUS)UNIV_CTSUS,
sum(UNIV_HRMUS)UNIV_HRMUS,
sum(UNIV_BBBCA)UNIV_BBBCA,
sum(UNIV_BBBMX)UNIV_BBBMX,
...

sum(FST_STR_SHOP_NBR_651)FST_STR_SHOP_NBR_651
from ANALYTICS_STG..LH_CAMPAIGN_DM_ADM1
where response = 1
;


drop table ANALYTICS_STG..LH_CAMPAIGN_DM_ADM3_pos;
-- 1 rows affected
create table ANALYTICS_STG..LH_CAMPAIGN_DM_ADM3_pos as
select 
b.UNIV_BBBUS/(sqrt(a.UNIV_BBBUS)*sqrt(4856590.0) ) UNIV_BBBUS,
b.UNIV_BABYUS/(sqrt(a.UNIV_BABYUS)*sqrt(4856590.0) ) UNIV_BABYUS,
b.UNIV_CTSUS/(sqrt(a.UNIV_CTSUS)*sqrt(4856590.0) ) UNIV_CTSUS,
b.UNIV_HRMUS/(sqrt(a.UNIV_HRMUS)*sqrt(4856590.0) ) UNIV_HRMUS,
b.UNIV_BBBCA/(sqrt(a.UNIV_BBBCA)*sqrt(4856590.0) ) UNIV_BBBCA,
b.UNIV_BBBMX/(sqrt(a.UNIV_BBBMX)*sqrt(4856590.0) ) UNIV_BBBMX,
b.UNIV_UNIV_BABYCA/(sqrt(a.UNIV_UNIV_BABYCA)*sqrt(4856590.0) ) UNIV_UNIV_BABYCA,
...

b.FST_STR_SHOP_NBR_651/(sqrt(a.FST_STR_SHOP_NBR_651)*sqrt(4856590.0) ) FST_STR_SHOP_NBR_651
from ANALYTICS_STG..LH_CAMPAIGN_DM_ADM1_pos a
cross join ANALYTICS_STG..LH_CAMPAIGN_DM_ADM2_pos b
;


drop table ANALYTICS_STG..LH_CAMPAIGN_DM_ADM4_pos;
-- 38,710,468 rows affected;
create table ANALYTICS_STG..LH_CAMPAIGN_DM_ADM4_pos as
select a.COUPON_BARCODE, a.RESPONSE,
a.UNIV_BBBUS*b.UNIV_BBBUS+
a.UNIV_BABYUS*b.UNIV_BABYUS+
a.UNIV_CTSUS*b.UNIV_CTSUS+
a.UNIV_HRMUS*b.UNIV_HRMUS+
a.UNIV_BBBCA*b.UNIV_BBBCA+
a.UNIV_BBBMX*b.UNIV_BBBMX+
...

a.FST_STR_SHOP_NBR_651*b.FST_STR_SHOP_NBR_651 score
from ANALYTICS_STG..LH_CAMPAIGN_DM_ADM1 a
cross join ANALYTICS_STG..LH_CAMPAIGN_DM_ADM3_pos b
;

select a.DECILE, count(*), sum(response)
from
(
select coupon_barcode, response, score,
NTILE(10) OVER(ORDER BY score DESC NULLS LAST) as DECILE
from ANALYTICS_STG..LH_CAMPAIGN_DM_ADM4_pos
) a
group by 1
order by 1 
;

PROC LOGISTIC to handle Randomized Complete Block Design (RCBD)

The main idea of blocking is to ensure that the random sample you draw from the population of interest does not, by random chance, have some undesirable properties. 

*Use risk scores (riskScr) and response groups (respGrp) to create blocks;
data blocks;
set s.dem0303;
if riskScr < 650 then riskGrp=1;
else if riskScr < 750 then riskGrp=2;
else riskGrp =3;
if respScr < 260 then respGrp=1;
else if respScr < 280 then respGrp=2;
else respGrp =3;
block = compress(riskGrp||respGrp);
random=ranuni(8675309);
run;

proc sort data=blocks;
by block;
run;

*randomly assigned units to treatments within each block;
proc rank data=blocks out=blocks groups=8;
by block;
var random;
ranks treatment;
run;

*Generate design matrix (factors: interestRate sticker priceGraphic);
data blocks;
set blocks;
array factors(3) interestRate sticker priceGraphic;
do i = 1 to dim(factors);
factors(i)=substr(put(treatment,binary3.),i,1);
end;
run;

*Analysis after collecting the data;
*Block | InterestRate | Sticker | PriceGraphic;
ods select globalTests type3;
title "Block | InterestRate | Sticker | PriceGraphic";
proc logistic data=s.dem0304 namelen=40;
class block interestRate sticker priceGraphic;
model response(event="1")=
block
interestRate|sticker|priceGraphic @1
interestRate|sticker|priceGraphic @2
interestRate|sticker|priceGraphic;
run;

Note:
It is possible to get the same parameter estimates in a slightly easier to read format. The above code uses the '@' notation in the MODEL statement to arrange the terms in the model by their polynomial order. A|B|C|D in a MODEL statement yields every term from the main effects up to the three-factor interaction. 

SAS Tricks: Padding zeros & Transpose Cross product sales

## Padding zeros
data expiredate;

do j=2010 to 2011;
if _n_=1 then exp_date=mdy(1,1,j);
i=1;
format exp_date monyy.;
exp_wk=year(exp_date)||put(i,z2.);
output;

do i=2 to int(365/7);
exp_date=exp_date+7;
format exp_date monyy.;
exp_wk=year(exp_date)||put(i,z2.);
output;
end;

end;

goptions ftext='Arial' ctext=BLACK htext=1.5 cells dev=activex;
ods html file="C:\expdate.html";
ods listing close;
options nodate nonumber;
proc print data=expiredate;
run;
ods html close;
ods listing;

## Transpose Cross product sales

data out3.multiple;

set out1.multiple;
length multiple_purchses $50 purchse_type $20;
array num(10) num_1110 num_1140 num_1160 num_1210 num_1610 num_1720 num_1750 num_1810 num_1820 num_1860;
multiple_purchses=put(num(1),3.);
do i=2 to 10;
multiple_purchses=compress(multiple_purchses||'-'||put(num(i),3.));
end;

purchse_type=put(min(num(1),1),3.);
do i=2 to 10;
purchse_type=compress(purchse_type||'-'||put(min(num(i),1),3.));
end;
run;

How SAS handles missing values

It is important to understand how SAS procedures handle missing data if you have missing data. To know how a procedure handles missing data, you should consult the SAS manual. Here is a brief overview of how some common SAS procedures handle missing data. 

• proc means
  For each variable, the number of non-missing values are used
• proc freq
  By default, missing values are excluded and percentages are based on the number of non-missing values. If you use the missing option on the tables statement, the percentages are based on the total number of observations (non-missing and missing) and the percentage of missing values are reported in the table.
• proc corr
  By default, correlations are computed based on the number of pairs with non-missing data (pairwise deletion of missing data). The nomiss option can be used on the proc corr statement to request that correlations be computed only for observations that have non-missing data for all variables on the var statement (listwise deletion of missing data).
• proc reg
  If any of the variables on the model or var statement are missing, they are excluded from the analysis (i.e., listwise deletion of missing data)
• proc factor
  Missing values are deleted listwise, i.e., observations with missing values on any of the variables in the analysis are omitted from the analysis.
• proc glm
  The handling of missing values in proc glm can be complex to explain. If you have an analysis with just one variable on the left side of the model statement (just one outcome or dependent variable), observations are eliminated if any of the variables on the model statement are missing. Likewise, if you are performing a repeated measures ANOVA or a MANOVA, then observations are eliminated if any of the variables in the model statement are missing. For other situations, see the SAS/STAT manual about proc glm.

*Get missing values freqencies of numeric variables:
proc freq data=raw;
table variables /missing;
run;

*Get missing values levels of categorical variables:
data class;
if 0 then set sashelp.class;
do i=1 to 10;output;
end;
stop;
run;

proc format;
value allmiss ._-.z=. other=1;
value $allmiss ' '=' ' other='1';
run;

ods select nlevels;
ods output nlevels=nlevels(keep=TableVar NNonMissLevels where=(NNonMissLevels=0));
proc freq levels;
format _character_ $allmiss. _numeric_ allmiss.;
run;
ods output close; 

Hypothesis Test in R and SAS

### R

The t.test( ) function produces a variety of t-tests.
Unlike most statistical packages, the default assumes unequal variance and applies the Welch df modification.

# one sample t-test
t.test(y,mu=3) # Ho: mu=3

# one sample median test
wilcox.test(y, mu=3)

# one sample binomial test
prop.test(sum(y), length(y), p=0.5)

# Chi-square goodness of fit
X <- matrix(c(172, 7, 6, 15), 2, 2)
chisq.test(table(X), p=c(10,10,10,70)/100)

# McNemar's Chi-squared test with continuity correction
X <- matrix(c(172, 7, 6, 15), 2, 2)
mcnemar.test(X)

# Wilcoxon-Mann-Whitney rank sum test
wilcox.test(y ~ x)

# Chi-square test
chisq.test(table(y, x))

# Fisher's exact test
fisher.test(table(y, x))

# Correlation
cor(y, x)
cor.test(y, x)

# Kruskal Wallis rank sum test
kruskal.test(y, x)

# Non-parametric correlation
cor.test(y, x, method = "spearman")

## Multual Information
mi.plugin(table(all$y, all$is_verified))

## Jaccard Similarity
1-dist(rbind(all$y, all$is_verified), method= "binary")

# Simple linear regression
lm(y ~ x)

# Multiple regression
lm(y ~ x1 + x2 + x3 + x4 + x5)

# independent 2-group t-test
t.test(y~x) # where y is numeric and x is a binary factor

# independent 2-group t-test
t.test(y1,y2) # where y1 and y2 are numeric

t.test(y1,y2, alternative="greater", var.equal=T)

# paired t-test
t.test(y1,y2,paired=TRUE) # where y1 & y2 are numeric

# Wilcoxon signed rank sum test
wilcox.test(y1, y2, paired = TRUE)

# One-way ANOVA
summary(aov(y ~ x))

# Factorial ANOVA
anova(lm(y ~ x1 * x2, data = hsb2))

# One-way repeated measures ANOVA
require(car)
require(foreign)
model <- lm(y ~ a + s, data = kirk)
analysis <- Anova(model, idata = kirk, idesign = ~s)
print(analysis)

# Repeated measures logistic regression
require(lme4)
glmer(y ~ x + (1 | id), data = exercise, family = binomial)

# Friedman test
friedman.test(cbind(x1, x2, x3))

# Simple logistic regression
glm(y ~ x, family = binomial)

# Factorial logistic regression
summary(glm(y ~ x1 * x2, data = hsb2, family = binomial))

# Multiple logistic regression
glm(y ~ x1 + x2, family = binomial)

# Analysis of covariance
summary(aov(y ~ x1 + x2))

# One-way MANOVA
summary(manova(cbind(x1, x2, x3) ~ x4))

# Multivariate multiple regression
require(car)
M1 <- lm(cbind(x1, x2) ~ x3 + x4 + x5 + x6, data = hsb2)
summary(Anova(M1))

# Canonical correlation
require(CCA)
cc(cbind(x1, x2), cbind(x3, x4))

# Factor analysis
require(psych)
fa(r = cor(model.matrix(~x1 + x2 + x3 + x4 + x5 - 1, data = hsb2)), rotate = "none", fm = "pa", 5)

# Principle component
princomp(~x1 + x2 + x3 + x4 + x5, data = hsb2)


Examples

1 Smoker vs NonSmoker

nonsmokers = c(18,22,21,17,20,17,23,20,22,21)
mean(nonsmokers)
sd(nonsmokers)
sqrt(var(nonsmokers)/length(nonsmokers))
smokers = c(16,20,14,21,20,18,13,15,17,21)
mean(smokers)
sd(smokers)
sqrt(var(smokers)/length(smokers))

plot(density(nonsmokers))            # output not shown
plot(density(smokers))               # output not shown
boxplot(nonsmokers,smokers,ylab="Scores on Digit Span Task",
names=c("nonsmokers","smokers"),
main="Digit Span Performance by\n Smoking Status")

2 Formula

mean.diff = mean(smokers) - mean(nonsmokers)
df = length(smokers) + length(nonsmokers) - 2
pooled.var = (sd(smokers)^2 * (length(smokers)-1) + sd(nonsmokers)^2 * (length(nonsmokers)-1)) / df
se.diff = sqrt(pooled.var/(length(smokers) + pooled.var/length(nonsmokers))
t.obt = mean.diff / se.diff
t.obt
p.value = 2*pt(t.obt,df=df)
p.value

3 odds ratios

odds.ratios <- function(clicks1, clicks2, nonclicks1, nonclicks2) {

#nonblank over blank

or = (clicks2/nonclicks2)/(clicks1/nonclicks1);

or.lower <- exp(log(or)-1.96*sqrt(1/clicks1 + 1/clicks2 + 1/nonclicks1 + 1/nonclicks2));

or.upper <- exp(log(or)+1.96*sqrt(1/clicks1 + 1/clicks2 + 1/nonclicks1 + 1/nonclicks2));

return(c(or, or.lower, or.upper))

}

#odds ratios

clicks1 <- as.numeric(data[2,3])

clicks2 <- as.numeric(data[3,3])

nonclicks1 <- as.numeric(data[2,2])-as.numeric(data[2,3])

nonclicks2 <- as.numeric(data[3,2])-as.numeric(data[3,3])

odds.ratios(clicks1, clicks2, nonclicks1, nonclicks2)

# Fisher's exact test

fisher.test(matrix(c(nonclicks1, clicks1, nonclicks2, clicks2), ncol = 2, byrow = T))

#### SAS

proc sql;

CONNECT TO NETEZZA (USER=*** PASSWORD=*** SERVER="***"DATABASE="ANALYTICS_STG");

create table data as

select * from connection to NETEZZA

(

select CAMPAIGN_ID,

SRCSYS_CELL_ID,

min(SRCSYS_CELL_NAME) SRCSYS_CELL_NAME,

SRCSYS_AUDIENCE_ID,

min(SRCSYS_AUDIENCE_NAME) SRCSYS_AUDIENCE_NAME,

count(distinct EMAIL_ADDRESS_ID) sent,

sum(case when upper(FEEDBACK_EVENT_CD)='OPEN' then 1 else 0 end) open,

sum(case when upper(FEEDBACK_EVENT_CD)='CLICK' then 1 else 0 end) click,

sum(case when upper(FEEDBACK_EVENT_CD)='CONVERSION' then 1 else 0 end) conv,

sum(case when upper(FEEDBACK_EVENT_CD)='UNSUB' then 1 else 0 end) unsub

from ANALYTICS_STG..LH_persado_email_and_resp

group by 1,2,4

order by 1,2,4;

)

;

quit;

data data1;

set data;

select(SRCSYS_AUDIENCE_ID);

when (841791043) control = 1;
...

when (834298683) control = 1;

when (834299793) control = 0;

when (834299803) control = 0;

when (834021323) control = 0;

when (834021333) control = 0;

when (834021373) control = 1;

otherwise control=.;

end;

if campaign_id = 9910 then delete;

open_rate=open/sent;

run;

/** One sample Binomial test;

proc freq data=data1;

by SRCSYS_AUDIENCE_ID;

table control / binomial(p=0.5);

exact binomial;

run;

**/



** Each Individual Campaigns;

** Two sample Binomial test;



proc ttest data=data1;

by campaign_id;

class control;

var open_rate;

run;

** Noparameteric Wilcoxon-Mann-Whitney test;

proc npar1way data=data1;

by campaign_id;

class control;

var open_rate;

run;

data data5;

set data1;

open_flag=1;

counts=open; output;

open_flag=0;

counts=sent-open; output;

run;

** Chisq Test;

proc freq data = data5;

by campaign_id;

weight counts;

tables control*open_flag / chisq;

run;



** Fisher's exact test;

proc freq data = data5;

by campaign_id;

weight counts;

tables control*open_flag / fisher;

run;



** McNemar test;

proc freq data=data5;

by campaign_id;

tables control*open_flag

exact mcnem;

weight counts;

run;




### All Campaigns

** One-way ANOVA;

proc glm data = data1;

class control campaign_id;

model open_rate = control campaign_id;

means control;

run;

quit;







** Paired t-test;

proc sql;

create table data2 as

select CAMPAIGN_ID, control,

sum(sent) as sent,

sum(open) as open,

sum(open)/sum(sent) as open_rate

from data1

group by CAMPAIGN_ID, control

;

run;

quit;

proc transpose data=data2 prefix=control out=data3;

by campaign_id;

id control;

var open_rate;

run;

proc ttest data = data3;

paired control0*control1;

run;







** Wilcoxon signed rank sum test;

data data4;

set data3;

diff = control0 - control1;

run;

proc univariate data=data4;

var diff;

run;

** McNemar test;

proc freq data=data5;

tables control*open_flag;

exact mcnem;

weight counts;

run;

proc freq data=data5;

by campaign_id;

tables control*open_flag;

exact mcnem;

weight counts;

run;



** logistic regression;

proc genmod data= data5 descending;

class campaign_id control;

freq counts;

model open_flag = campaign_id control/ dist=binomial link=logit;

run;

proc logistic data= data5 descending;

class campaign_id control;

freq counts;

model open_flag = campaign_id control/ scale = none aggregate;

output out=plot p=p l=l u=u;

run;



proc sgplot data=plot;

band x=control

lower=l

upper=u /group=campaign_id

transparency=.5;

series x=control

y=p /group=campaign_id;

run;

INFORMAT AND FORMAT

INFORMAT gives SAS special instructions for reading a variable.

FORMAT gives SAS special instructions for writing a variable.

Character to Numeric: newvar=input(oldvar, informat);
The informat must be the type you are converting to-numeric;

Numeric to Character: newvar=put(oldvar, format);
The format must be the type you are converting from-numeric. 

PROC CLUSTER, PROC FASTCLUS to CLUSTERING

The following SAS codes is to use Proc cluster to generate random seeds (centroids/initial points), and then use Proc fastclus to create the clusters.

%let demo=var1 var2 ;

title2 'Hierarchical Solution (WARD''S)';
proc cluster data=out1.training method=ward k=10 trim=0.1 outtree=tree noprint;
var &demo;
copy keys;
run;

proc tree data=tree nclusters=5 dock=5 out=out1.results noprint;
copy &demo;
run;

proc freq data=out1.results ;
table cluster;
run;

/* generate the centroids of the hierarchical clusters */
title1 'Cluster Centroids';
proc means data=out1.results;
class cluster;
var &demo;
output mean= out=centroids(where=(_type_ = 1));
run;

title1 'Score Development Data against the Centroids';
proc fastclus data=out1.training seed=centroids maxclusters=5 least=2 out=results noprint;
var &demo;
run;

title1 'USS (5 clusters)';
proc means data=results uss;
var distance;
run;

proc freq data=results;
table cluster;
run; 

%missingPattern to Studying Missing Data Patterns

The macro is designed to look at missing data in four ways: the proportion of subjects with each pattern of missing data, the number and percentage of missing data for each individual variable, the concordance of missingness in any pair of variables, and possible unit nonresponse.

The SAS macro is %missingPattern:

%missingPattern(datain=, varlist=, exclude=, missPattern1=, dataout1=, missPattern2=, dataout2=, missPattern3=, dataout3=, missPattern4=, dataout4=)

Reference: T. Schwartz, Q. Chen, and N. Duan, (2011). Studying missing data patterns using a SAS macro, SAS Global Forum 2011 proceedings.

Example:
data tmp;
set out1.training_motorcycle1;
if vn0451^=. then vn0451_log=log(1+vn0451);
else vn0451_log=.;
if vn0467^=. then vn0467_log=log(1+vn0467);
else vn0467_log=.;
if vn0479^=. then vn0479_log=log(1+vn0479);
else vn0479_log=.;
if vn0712^=. then vn0712_log=log(1+vn0712);
else vn0712_log=.;
if vn0717^=. then vn0717_log=log(1+vn0717);
else vn0717_log=.;
if vn0722^=. then vn0722_log=log(1+vn0722);
else vn0722_log=.;
run;

%let varlist=vn0451_log vn0467_log vn0479_log vn0712_log vn0717_log vn0722_log;

%missingPattern(datain=tmp, varlist=&varlist, exclude='FAULT', missPattern1='TRUE',
dataout1=result1,missPattern2='TRUE', dataout2=result2,
missPattern3='TRUE', dataout3=result3, missPattern4='TRUE', dataout4=result4);

ods html path = "/u/lhuang/" (url = none)
body = "temp.html"
style = Default;
Ods listing close;

proc print data=result1;
proc print data=result2;
proc print data=result3;
*proc print data=result4;
run;

* DOWNLOAD HTML FILES;
proc download infile="/u/lhuang/temp.html"
outfile="&LocalLocation/missing pattern.html";
run; 

PROC CLUSTER, PROC FASTCLUS to run Variable Selection for Clustering Analysis

Variable Standardization

Standardization refers to data transformation that involves correction of variables using either means or standard deviations or both, depending on the data and context of analysis. SAS PROC STDIZE allows users to standardize a set of variables using several criteria, including the following methods:

RANGE standardization is most helpful when variables are measured on different scales. Variables with large values and a wide range of variation have significant effects on the final similarity measure. Hence, it is essential to make sure that each variable is evenly constituted in the distance measurement by means of data standardization.

MEAN/MEADIAN offers centering refers to variables being adjusted using only the mean or median across the variables. The variable mean or median is subtracted from its original score. Standard deviation is not adjusted in this process.

STD standardization refers to correction of variables using the variable mean and standard deviation. The variable mean is subtracted from its original score and then divided by the standard deviation.

L (p) refers to adjusting the scores using Minkowski distance. For example, if p=2, then Euclidean distance is applied.

proc stdize data=dataset method=range out=outdataset outstat=stats;

var &varlist;

run;

Standardization to unit variability can be seen as a special case of weighting, where the weights are the reciprocals of variability. In some cases, however, defining weights as inversely proportional to some measure of total variable can actually dilute the difference between clusters. This is major reason why weights based on sample range are usually more effective when clustering than weights based on the standard deviation.

Clustering Methodology

Types of Models

SAS offers different clustering algorithms including k-means clustering, non-parametric clustering, and hierarchical clustering.

1. k-means clustering is, perhaps, the most popular partitive clustering algorithm. One reason for its popularity is that the time required to reach convergence is proportion to the number of observations being clustered, which means it can be used to cluster larger data sets. However, k-means clustering is inappropriate for small data sets (<100 cases). The solution becomes sensitive to the order of the observations and the variables, which is known as order effect.

title 'K-Means Clustering using Adaptive Training';

proc fastclus data=dataset maxclusters=5 maxiter=100 least=2 drift replace=full distance out=cluster;

var &varlist;

run;

Option MAXCLUSTERS=specifies the maximum number of clusters allowed. Instead of using MAXCLUSTERS=, option RADIUS= is also able to establishe the minimum distance criterion for selecting new seeds. No observation is considered as a new seed unless its minimum distance to previous seeds exceeds the value given by the RADIUS= option. The default value is 0.

Option MAXITER= specifies the number of iterations to facilitate convergence. By default, MAXITER=1. When the value of the MAXITER= option is greater than 0, each observation is assigned to the nearest seed, and the seeds are recomputed as the means of the clusters.

Option LEAST= specifies distance metrics. By default, LEAST=2, which is Euclidean distance. LEAST=1 implements city block distance, and LEAST=MAX implements maximum absolute deviation. In other words, PROC FASTCLUS requires numeric data, and be sensitive to extreme values, so standardization in data preparation is important.

Option DRIFT is specified the closest seed moves as each case is assigned to it.

Option REPLACE= specifies how seed replacement is performed. Option REPLACE=FULL requests default seed replacement. Option REPLACE=PART requests seed replacement only when the distance between the observation and the closest seed is greater than the minimum distance between seeds. Option REPLACE=NONE suppresses seed replacement. Option REPLACE=RANDOM selects a simple pseudo-random sample of complete observations as initial cluster seeds.

Option DISTANCE computes distances between the cluster means.

2. Nonparametric methods can detect clusters of unequal size and dispersion, or clusters with irregular shapes. If the covariance matrices are unequal or radically non-normal, nonparametric density estimation is often the best approach. And nonparametric methods are less sensitive than most clustering techniques to changes in scale.

title 'Nonparametric Clustering';

proc modeclus data=dataset method=1 r=5.75330 join out=nonpar_results;

var &varlist;

run;

Option R= specifies the radius of the sphere of support for uniform-kernel density estimation and the neighborhood for clustering. Option METHOD= specifies what clustering method to use. For most purposes, METHOD=1 is recommended.

3. Hierarchical methods form the backbone of cluster analysis. The popularity of hierarchical method is partly due to the fact that they are not subject to the order effect. And also, some hierarchical methods can even recover irregular clusters directly. Unfortunately, hierarchical methods often require processing times on the order of the square. This limits their use to small and mid-sized data sets.

title2'Hierarchical Solution (Average)';
proc cluster data=dataset method=average outtree=tree simple rmsstd rsquare;
var &varlist;
copy customer_key;
run;

Option SIMPLE displays simple, descriptive statistics.

The METHOD= specification determines the clustering method used by the procedure. The example used average method. Other popular methods include single, centriod, twostage, and ward.

Option RMSSTD displays the pooled standard deviation of all the variables of each cluster. Since the objective of cluster analysis is to form homogeneous groups, the RMSSTD of a cluster should be as small as possible.

Option RSQUARE displays the R-squared and semi-partial R-squared to evaluate cluster solution. R-squared measures the extent to which groups or clusters are different from each other (so, when you have just one cluster R-squared value is, intuitively, zero). Thus, the R-squared value should be high. semi-partial R-squared is the loss of homogeneity due to combining two groups or clusters to form a new group or cluster. Thus, the semi-partial R-squared value should be small to imply that we are merging two homogeneous groups.

proc tree data=tree nclusters=5 dock=5 out=results2 noprint;
copy &varlist;
run;

PROC TREE procedure produces a tree diagram, also known as a dendrogram or phenogram, using a data set created by PROC CLUSTER. PROC CLUSTER creates output data sets that contain the results of hierarchical clustering as a tree structure. PROC TREE uses the output data set to produce a diagram of the tree structure.

Option NCLUSTERS= specifies the number of clusters desired in the OUT= data set.

The COPY statement specifies one or more character or numeric variables to be copied to the OUT= data set.

Sample Size Calculation using PROC POWER for proportions under CRD

Q: How to ensure you get the correct sample size?

A: Leverage power analysis to determine sample size.

One of the pivotal aspects of planning an experiment is the calculation of the sample size. It is naturally neither practical nor feasible to study the whole population in any study. Hence, a set of users is selected from the population, which is less in number (size) but adequately represents the population from which it is drawn so that true inferences about the population can be made from the results obtained. This set of individuals is known as the “sample" int DOE (design of experiment).

It is a basic statistical principle with which we define the sample size before we start a study so as to avoid bias in interpreting results. If we include very few users in a study, the results cannot be generalized to the population as this sample will not represent the size of the target population. Further, the study then may not be able to detect the difference between test groups, making the study unethical. On the other hand, if we study more subjects than required, we put more users than needed, also making the study wasting precious resources, including the researchers’ time.

Generally, the sample size for any study depends on the:

•  Acceptable level of significance
•  Power of the study
•  Expected effect size
•  Expected variance

Level of Significance

This is the pre-defined threshold value to reject the null hypothesis. Usually p<0.05 will be taken as statistically significant to accept that the result is observed due to chance. To put in different words, it is possible to accept the detection of a difference 5 out of 100 times when actually no difference exists (i.e., get a false positive result). It is denoted by letter α(alpha) as Type 1 error.

Power

Exactly conversely, type II of error indicate failing to detect a difference when actually there is a difference (i.e., false negative result). The false negative rate is the proportion of positive instances that were erroneously reported as negative and is referred to in statistics by the letter β(beta). The power of the study then is equal to (1 –β) and is the probability of failing to detect a difference when actually there is a difference. The power of a study increases as the chances of committing a Type II error decrease. Usually most studies accept a power of 80%. This means that we are accepting that one in five times (that is 20%) we will miss a real difference. 

Expected effect size

The effect size is the minimum deviation from the null hypothesis that the hypothesis test set up to detect. Suppose you are comparing different responses between treatment group and control group, and the measuring metrics are group means denoted as μ1 and  μ2,  the expected effect size equals to abs( μ1- μ2). Generally speaking, the smaller the difference you want to detect, the larger the required sample size.

Expected variance 

The variance describes the variability/noise among measurement units. As variance gets larger, it gets harder to detect a significant difference, so the bigger sample size requires. The pilot study or similar study could help to determine the expected variance.

Reference

1. Design and Analysis of Experiments. by Montgomery, D.C.X (Hoboken, New Jersey: John Wiley & Sons, Inc., 2000)

* Power with One Sample; 

proc power;
oneSampleFreq
/* hypotheses of interest */
nullProportion=0.1
proportion=0.2
sides=1
/* decision rule */
alpha=0.3
/* sample size */
nTotal=10 20 25 30 40 50 100
power=.
;
run;

*normal approximation;
proc power;
oneSampleFreq
/* hypotheses of interest */
nullProportion=0.1
proportion=0.2
sides=1
/* decision rule */
alpha=0.3
/* solve for sample size */
nTotal=.
power=.80
/* different distributional assumptions */
test=z
method=normal
;
run;

* Power with Two Samples;
** SAS
proc power;
twoSampleFreq
/* hypotheses of interest */
refProportion=0.01
proportionDiff=0.001/*first factor difference*/ .0025/*second factor difference*/
sides=1 2
/* decision rule */
alpha=0.05
/* sample size */
nTotal=.
power=.6 .99
;
plot y=power
yopts=(crossref=YES ref=.8)
vary(color by proportionDiff,
symbol by sides);
run;

proc power;
twoSampleFreq
/* hypotheses of interest */
refProportion=0.1
proportiondiff=0.1
sides=1
/* decision rule */
alpha=0.05
/* sample size */
nTotal=.
power=.6 .8 .99
/* how balanced the test is */
groupWeights=(1 1) (10 15) (1 2) (1 3) (1 10)
;
plot;
run;

** R
install.packages("pwr")
require(pwr)

delta <- 20
sigma <- 60
d <- delta/sigma
pwr.t.test(d=d, sig.level=.05, power = .80, type = 'two.sample')

PROC POWER, PROC GLMPOWER for proportions under CRD

* Balance with Two Factor; 

proc power;
twosamplefreq
refproportion= 0.010
proportiondiff=0.001 0.0025
ntotal=.
power=.8;
run;

data work.a;
input intro $1-4
goto $6-9
responseRate;
variance= responseRate*(1-responseRate);
datalines;
LOW LOW .0135
LOW HIGH .0125
HIGH LOW .011
HIGH HIGH .010
;
run;

proc glmpower data=work.a;
class intro goto;
model responseRate=intro|goto;
power
power=.8
ntotal=.
stddev=%sysfunc(sqrt(.01*.99))
%sysfunc(sqrt(.011*.989));
run;

* Three Two-Level Factors;
* 8 trts;
data work.dem0302;
do a=-1 to 1 by 2;
do b=-1 to 1 by 2;
do c=-1 to 1 by 2;
respRate=.01+.00050*sum(a,b,c);
output;
end;
end;
end;
run;

proc print data=work.dem0302;
run;

proc glmpower data=work.dem0302;
class a b c;
model respRate=a|b|c;
power
power=.8
ntotal=.
stddev=%sysfunc(sqrt(.01*.99));
run;

*Analysis Examples:
data work.ex0302;
informat responseRate percent6.4;
input interestRate $1-4
sticker $6-8
priceGraphic $10-14
responseRate;
datalines;
LOW YES SMALL 1.00%
LOW YES LARGE 1.20%
LOW NO SMALL 0.85%
LOW NO LARGE 1.00%
HIGH YES SMALL 0.80%
HIGH YES LARGE 1.00%
HIGH NO SMALL 0.60%
HIGH NO LARGE 0.75%
;
run;

*N=4,973,040;
proc glmpower data=work.ex0302;
class interestRate Sticker priceGraphic;
model responseRate=interestRate|Sticker|priceGraphic;
power
power=.8
ntotal=.
stddev=%sysfunc(sqrt(.01*.99));
run;

proc logistic data=s.two3;
class interestRate Sticker priceGraphic;
model orders/volume=
interestRate|Sticker|priceGraphic @1
interestRate|Sticker|priceGraphic @2
interestRate|Sticker|priceGraphic;
output out=work.plot p=p l=l u=u;
run; 

proc QLIM, proc LOGISTIC, PROC GENMOD, PROC GLIMMIX to fit Population-averaged model vs. Subject-specific model

The GLM consists of three elements:

1. A probability distribution from the exponential family.

2. A linear predictor η = Xβ .

3. A link function g such that E(Y) = μ = g^-1(η).

For GLM, there are two different models: population-averaged model without mix effects, subject-specific model with mix effects.

The population-averaged approach is focused on modeling the mean response across the population of units at each time point as a function of time. Thus, the model described how the averages across the population of responses at different time points are related over time.

The subject-specific approach is focused on modeling the individual unit trajectories rather than the mean across all units. We did this by the introduction of random effects, e.g., the random coefficient model that says each unit has its own intercept and slope.

The MLE of beta is the solution of the normal equation, (sum_i x_i inv(sigma) (y_i-x_i*beta)=0). It is consistent and asymptotically normal with estimated asymptotic variance matrix H. There are three types of variances can be computed.
First, Hessian matrix H=-der^2.(l(beta))/der(beta)der(beta')
Second, asymptotic variance matrix of outer product B=der.(l(beta))/der(beta). der.(l(beta))/der(beta').
Third, asymptotic robust-sandwich estimate H^(-1)BH^(-1).

proc QLIM compute all three variances.
proc Logistic compute Hessian variance.
proc GENMOD and proc GLIMMIX compute Hessian variance and robust-sandwich variances. 

Proc QLIM to fit Tobit Models

Standard Tobit Model:

yi=yi^star*I(yi^star>0), where yi^star=xi*beta+ei, ei~N(0,sigma^2).
The inverse Mills ratio, also called the selection hazard, is used to take account of a possible selection bias.
IMratio=f(xi*beta/sigma)/Pi(xi*beta/sigma), f is pdf, Pi is cdf of normal distribution. As the probability of censoring increases, the ratio approaches infinity. As the probability of censoring decreases, the ratio approaches 0.

The likelihood function of Tobit models consists of two parts: the first part is the likelihood function for OLS regression and the second part is the likelihood function for a Probit model. If there was no censoring, then the tobit model would produce unbiased OLS estimates. But if there was censoring, the Tobit model would not only predict the non-censored response values, but also the estimated probability that the observation is censored.

Probit Model:
Pi^(-1)(pi)=xi*beta.
The probit link function transforms a probability to a standard normal z-score at which the left-tailed probability equals the posterior probability. Pi^(-1)(0.5)=0, Pi^(-1)(1.96)=.0975.

Proc QLIM ( Qualitative and Limited dependent variable Model) analyzes univariate and multivariate limited dependent variables where dependent variables are observed in a limited range of values.
Option ENDOGENOUS CENSORED specifies the dependent variable is censored.
Option HETERO specifies the ways of heteroscedasticity of the residuals.
Option ENDOGENOUS TRUNCATED specifies the dependent variable is truncated. In a truncated distribution, the values of the predictor variables are known only when the response variable is observed. This differs from censoring where the predictor variables are known even when the response variable is unknown.

* Tobit Model in PROC QLIM;
** ENDOGENOUS CENSORED (lower bound=0);
proc qlim data=pva;
class gender;
format gender $gender.;
model donation = gender age income_group wealth_rating pep_star months_since_last_gift median_home_value recent_avg_gift_amt recent_card_response_prop card_prom_12 recent_response_count;
endogenous donation ~ censored(lb=0);
output out=selection conditional expected marginal predicted xbeta mills;
title 'Tobit Model in PROC QLIM';
run;

* Tobit Model Corrected for Heteroscedasticity;
**HETERO specified variance function Var(ei);
proc qlim data=pva;
class gender;
format gender $gender.;
model donation= gender age income_group wealth_rating
pep_star months_since_last_gift
median_home_value recent_avg_gift_amt
recent_card_response_prop card_prom_12
recent_response_count;
endogenous donation ~ censored(lb=0);
hetero donation~pep_star months_since_last_gift
recent_avg_gift_amt card_prom_12
recent_response_count;
output out=selection xbeta;
title "Tobit Model Corrected for Heteroscedasticity";
run;

* Tobit Model in PROC QLIM;
** ENDOGENOUS CENSORED (lower bound=500 upper bound=1500);
proc qlim data=censored_aids;
model basecd4 = age cigarettes drug partners depression;
endogenous basecd4 ~ censored(lb=500 ub=1500);
title 'Censored Distribution in PROC QLIM';
run;

* Tobit Model in PROC QLIM;
** ENDOGENOUS CENSORED (lower bound=500 upper bound=1500);
proc qlim data=censored_aids;
model basecd4 = age cigarettes drug partners depression;
endogenous basecd4 ~ censored(lb=500 ub=1500);
hetero basecd4 ~ cigarettes;
output out=selection xbeta;
title 'Censored Distribution in PROC QLIM';
run;

* Tobit Model in PROC QLIM;
**Truncated Regression Model;
proc qlim data=housing;
model median_value = crime_rate large_lots nonretail charles_river nitric_oxide
average_number_rooms percent_lower_status
distance_Boston teacher_pupil_ratio access_highway;
endogenous median_value ~ truncated (ub=25);
output out=selection conditional marginal predicted mills xbeta errstd;
title "Truncated Regression Model for the Boston Housing Data";
run;

* Tobit Model in PROC QLIM;
*Truncated Regression Model Accounting for Heteroscedasticity;
proc qlim data=housing;
model median_value = crime_rate large_lots nonretail charles_river nitric_oxide
average_number_rooms percent_lower_status
distance_Boston teacher_pupil_ratio access_highway;
endogenous median_value ~ truncated (ub=25);
hetero median_value~ percent_lower_status distance_Boston;
output out=selection xbeta;
title "Truncated Regression Model for the Boston Housing Data Accounting for Heteroscedasticity";
run;

Sample selection Models;
In this model, the response variable is only observed when some selection criterion is met. The selection equation has a latent response variable which includes predictor variable, coefficients and an error term.

* Tobit Model in PROC QLIM;
**Sample Selection Model ;
proc qlim data=mroz;
model labor_force = age education experience kidslt6 income / discrete;
model wage = experience experience_sq education marg_tax_rate / select(labor_force=1);
hetero wage ~ experience education;
title "Sample Selection Model of Married Female Wages in Labor Force";
run; 

PROC DISCRIM to fit Linear Discriminant Analysis (LDA) , Quadratic Discriminant Analysis (QDA)

LDA

LDA attempts to express one dependent variable as a linear combination of other features or measurements.

LDA explicitly attempts to model the difference between the classes of data.

LDA for two classes

PROC DISCRIM DATA=iris; 

CLASS Species;
Var X1-XK;

RUN;

PROC DISCRIM DATA=Train TESTDATA=Test TESTOUT=Pred;
CLASS Species;
Var X1_XK;
RUN;

• Multivariate normal distribution assumptions holds for the response variables. This means that each of the dependent variables is normally distributed within groups, that any linear combination of the dependent variables is normally distributed, and that all subsets of the variables must be multivariate normal. 
• Each group must have a sufficiently large number of cases.
• Different classification methods may be used depending on whether the variance-covariance matrices are equal (or very similar) across groups.

Difference btw LDA and logistic regression

• LDA operates by maximizing the log-likelihood based on an assumption of normality and homogeneity 
• Logistic regression makes no assumption about Pr(X), and estimates the parameters of Pr(G|x) by maximizing the conditional likelihood
• Intuitively, it would seem that if the distribution of x is indeed multivariate normal, then we will be able to estimate our coefficients more efficiently by making use of that information by using LDA.
• On the other hand, logistic regression would presumably be more robust if LDA’s distributional assumptions are violated

QDA

QDA is assumed that the measurements from each class are normally distributed.When the normality assumption is true, the best possible test for the hypothesis that a given measurement is from a given class is the likelihood ratio test.

Proc LOGISTIC, PROC GENMOD, PROC GLIMMIX to fit Binary Logit, Cumulative Logit and Cumulative Probit Models

Binary Logit Models

Suppose there are J-levels of the outcome Yi, 1..J. The cumulative probabilities of Yi, pj(xi)=Pr(Yi<=j|xi), reflect the ordering, with p1(xi)<=p2(xi)<=...<=pJ(xi)=1.

Proc logistic, genmod and glimmix with the option link=cumlogit will fit the cumulative logit model s.t. log(pj(xi)/(1-pj(xi))=aj+xi.beta. The parameter beta describe the effect of a covariate on the log odds of response in the category j or below.

proc logistic data=test;
class vn0476_M(ref='female') vn0435_34(ref='rent') vn0455_1(ref='unmarried')/param=ref;
model product_pref=&numeric &norminal/link=cumlogit;
format vn0476_M gender. vn0435_34 rent. vn0455_1 marital.;
run;

Cumulative Probit Models

Proc logistic, genmod and glimmix with the option link=cumprobit will fit the cumulative probit model s.t. Pi^(-1)=ai+xi.beta.

proc logistic data=test;
class vn0476_M(ref='female') vn0435_34(ref='rent') vn0455_1(ref='unmarried')/param=ref;
model product_pref=&numeric &norminal/link=cumprobit;
format vn0476_M gender. vn0435_34 rent. vn0455_1 marital.;
run;

Cumulative logit and cumulative probit models assume the effect of a covariate is identical for all J-1 cumulative logits, which is called the proportional odds property. 

PROC QLIM to fit Ordered Logit Models and Ordered Probit Models

Suppose that uj, 0..J such that -inf=u0<u1<...<uJ=inf. The observed outcome is a categorization of a latent variable Yi^star=xi.beta+ei s.t. Yi=j iff uj-1<Yi^star<=uj. The probability of response is pj(xi)=Pr(Yi=j|xi)=F(uj-xi.beta)-F((uj-1)-xi.beta).

In the ordered logit model s.t. log(pj(xi)/(1-pj(xi))=uj-xi.beta. The parameter beta describe the effect of a covariate on the log odds of response in the category j, or marginal effect of the covariate E(Yi^star|xi).
In the cumulative probit model s.t. Pi^(-1)=uj-xi.beta.

*Ordered Logit;
proc logistic data=test desc;
class vn0476_M(ref='female') vn0435_34(ref='rent') vn0455_1(ref='unmarried')/param=ref;
model product_pref=&numeric &norminal;
format vn0476_M gender. vn0435_34 rent. vn0455_1 marital.;
run;

Proc qlim fits the ordered logit and ordered probit models. It uses the latent variable formulation. By default, an intercept is included in beta and the first threshold parameter u1=0. The model option limit1=varying overrides the default.

proc qlim data=test covest=qml;
class vn0476_M vn0435_34 vn0455_1;
endogenous product_pref ~ discrete (dist=logistic order=formatted);
*hetero product_pref ~&numeric;
model product_pref=&numeric &norminal/limit1=varying;
format vn0476_M gender. vn0435_34 rent. vn0455_1 marital.;
run;

proc qlim fit the equivalent homoscedastic ordered logit model. However, the signs for the covariates for intercepts are reversed.

Ordered logit and ordered probit models assume the effect of a covariate is identical for all J-1 cumulative logits, which is called the proportional odds property. 

PROC LOGISTIC and PROC PHREG to fit Conditional Logistic Regression

The usual MLE sometimes is not appropriate. For example, there may be insufficient sample size for logistic regression, particularly if the data are highly stratified and there are a small number of subjects in each stratum. Highly stratified data often come from a design with cluster sampling, that is designs with tow or more observations for each primary sampling unit or cluster. The appropriate form of logistic regression for these types of data is called conditional logistic regression.

You can fit a model based on conditional probabilities that condition away the random cluster effects, which results in a model that contains substantially fewer parameters. The random cluster effects, e.g., centers in clinical trials, or dealers in auto risk analysis. Those are nuisance parameters.

Since the likelihood conditioned on the discordant pairs, the concordant pairs are non-informative and thus can be ignored. Considering Pr(yi1=0, yi2=1) and Pr(yi1=1, yi2=0), the conditional likelihood for the entire data will reduce. The explanatory variables are the differences in values of the explanatory variables for the treatment and control i.
data trial_test;
drop center1 i_sex1 age1 initial1 improve1 trtsex1 trtinit1
trtage1 isexage1 isexint1 iageint1;
retain center1 i_sex1 age1 initial1 improve1 trtsex1 trtinit1
trtage1 isexage1 isexint1 iageint1 0;
input center treat $ sex $ age improve initial @@;
/* compute model terms for each observation */
i_sex=(sex='m'); i_trt=(treat='t');
trtsex=i_sex*i_trt; trtinit=i_trt*initial;
trtage=i_trt*age; isexage=i_sex*age;
isexinit=i_sex*initial;iageinit=age*initial;
/* compute differences for paired observation*/
if (center=center1) then do;
pair=10*improve + improve1;
i_sex=i_sex1-i_sex;
age=age1-age;
initial=initial1-initial;
trtsex=trtsex1-trtsex;
trtinit=trtinit1-trtinit;
trtage=trtage1-trtage;
isexage=isexage1-isexage;
isexint=isexint1-isexinit;
iageinit=iageint1-iageinit;
if (pair=10 or pair=1) then do;
/* output discordant pair observations */
improve=(pair=1); output trial_test; end;
end;
else do;
center1=center; age1=age;
initial1=initial; i_sex1=i_sex; improve1=improve;
trtsex1=trtsex; trtinit1=trtinit; trtage1=trtage;
isexage1=isexage; isexint1=isexinit; iageint1=iageinit;
end;
cards;
1 t f 27 0 1 1 p f 32 0 2
2 t f 41 1 3 2 p f 47 0 1
3 t m 19 1 4 3 p m 31 0 4
4 t m 55 1 1 4 p m 24 1 3
5 t f 51 1 4 5 p f 44 0 2
6 t m 23 0 1 6 p f 44 1 3
;

proc logistic data=trial descending;
model improve = initial age i_sex isexage isexinit iageinit trtsex trtinit trtage /
selection=forward include=3 details;
run;

And also, you can use proc phreg for conditional logistic regression so that you can operate directly on the actual observations; you don't have to create difference observations.
data trial2;
drop center1 i_sex1 age1 initial1 improve1 trtsex1 trtinit1
trtage1 isexage1 isexint1 iageint1;
retain center1 i_sex1 age1 initial1 improve1 trtsex1 trtinit1
trtage1 isexage1 isexint1 iageint1 0;
input center treat $ sex $ age improve initial @@;
/* compute model terms for each observation */
i_sex=(sex='m'); i_trt=(treat='t');
trtsex=i_sex*i_trt; trtinit=i_trt*initial;
trtage=i_trt*age; isexage=i_sex*age;
isexinit=i_sex*initial;iageinit=age*initial;
cards;
1 t f 27 0 1 1 p f 32 0 2
2 t f 41 1 3 2 p f 47 0 1
3 t m 19 1 4 3 p m 31 0 4
4 t m 55 1 1 4 p m 24 1 3
5 t f 51 1 4 5 p f 44 0 2
6 t m 23 0 1 6 p f 44 1 3
;

proc phreg data=trial2 nosummary;
strata center;
model improve = initial age i_sex i_trt trtsex trtinit trtage isexage isexinit iageinit/ ties=discrete details;
run;

proc logistic data=trial2;
class sex treat / param=ref;
strata center;
model improve = initial age i_sex i_trt trtsex trtinit trtage isexage isexinit iageinit;
run;