Study Notes2 Predictive Modeling with Logistic Regression: collapse levels of contingency table

/* In the contingency table if there are many levels for a variable, then at some levels */

/* of this variable there may be no events(or all events) happening. This is called quasi */

/* completion problem(all separation). If this happens, the MLE will not exist since the */

/* regression coefficient will be infinit if want to maximize the likelihood */

/* A good solution is to collapse(cluster) some levels for the given categorical variable */

/*Variable Cluster_code has many levels which cause 0 events at some levels of cluster_code*/

/*This is called Quasi-completion problem*/

proc means data=pva1 nway missing;

class cluster_code;

var target_b;

output out=level mean=prop;

run;

proc print data=level width=min;

run;

ods output clusterhistory=cluster;

/* method=ward will use Greenacre's method to collapse levels of contingency tables */

/* the levels are clusterd based on the redunction in the chisq test of association */

/* levels with similar marginal response rate will be merged */

/* levels with least chisq value decreased will be merged, this is from the */

/* Semipartial R-Square value of the proc cluster output */

proc cluster data=level method=ward outtree=fortree plots=(dendrogram(vertical height=rsq));

freq _freq_;

var prop;

id cluster_code;

run;

proc freq data=pva1;

tables cluster_code*target_b / chisq;

output out=chi(keep=_pchi_) chisq;

run;

/* _pchi_ = 112.009 */

/* # of levels final is the one that giving least log(p_value) */

data cutoff;

if _n_ = 1 then set chi;

set cluster;

chisquare=_pchi_*rsquared;

degfree=numberofclusters-1;

logpvalue=logsdf('CHISQ',chisquare,degfree);

run;

proc sgplot data=cutoff;

scatter y=logpvalue x=numberofclusters / markerattrs=(color=blue symbol=circlefilled);

xaxis label="Number of Clusters";

yaxis label="Log of p-value";

title "Plot of the Log of p-value by Number of Clusters";

run;

/* number of clusters = 5 having the minimum logpvalue */

proc sql;

select numberofclusters into :ncl

from cutoff

having logpvalue=min(logpvalue);

quit;

ods html close;

proc tree data=fortree nclusters=&ncl out=clus;

id cluster_code;

run;

ods html;

proc sort data=clus;

by clusname;

run;

data clus_fmt;

retain fmtname "clus_fmt";

set clus;

start=cluster_code;

end=cluster_code;

label=clusname;

run;

/* at last use proc format to format the variable levels into clusters */

/* The formated fmt_cluster_code can be used as a categorical variable */

proc format cntlin=clus_fmt;

select clus_fmt;

run;

data pva1;

set pva1;

fmt_cluster_code=put(cluster_code,clus_fmt.);

run;

proc freq data=pva1;

table fmt_cluster_code / missing list;

run;

Study Notes1 Predictive Modeling with Logistic Regression: impute missing data with proc stdize for continuous variable

/*This is the study notes from SAS online training. It includes why use Logistic Regression, how to clean the data */

/*(impute missing value, cluster rare events levels, variable clustering, variable screening), how to build */

/*logistic regression and how to measure the performance of the model */

libname mydata "D:\SkyDrive\sas_temp";

proc datasets lib=mydata;

contents data=_all_;

run;

*** How to impute the missing data by proc stdize ***;

/*data preparation*/

data pva(drop=control_number);

set mydata.pva_raw_data;

run;

/*use proc means with nmiss to check how many obs are missing for each variable*/

proc means data=pva nmiss min max median;

var donor_age income_group wealth_rating;

run;

/*use array to set indicator for missing obs, if missing then indicatd as 1*/

data pva;

set pva;

array a_mi{*} mi_donor_age mi_income_group mi_wealth_rating;

array a_var{*} donor_age income_group wealth_rating;

do i=1 to dim(a_mi);

a_mi(i)=(a_var(i)=.);

end;

run;

/*group data into 3 groups by recent_response_prop from low to high: first 1/3 grp_resp is 0, next 1/3 of data is 1 */

/*the last 1/3 is 2. In the same way for grp_amt. So totally there are 9 groups of data considering grp_resp and grp_ant */

/*obs numbers in grp_resp level or grp_amt level are similar, but at grp_resp*grp_amt level is different */

proc rank data=pva out=pva groups=3;

var recent_response_prop recent_avg_gift_amt;

ranks grp_resp grp_amt;

run;

proc freq data=pva;

table grp_resp grp_amt / missing list;

table grp_resp*grp_amt / missing list;

run;

/*sort the data by grp_resp and grp_amt*/

proc sort data=pva;

by grp_resp grp_amt;

run;

/*impute the missing data in each group formed by grp_resp*grp_amt, by the median value of non-missing data in that group*/

/*after the imputation, the data will have imputed value as well as the missing value indicator */

proc stdize data=pva method=median reponly out=pva1;

by grp_resp grp_amt;

var donor_age income_group wealth_rating;

run;

/*check the imputed value in each group*/

proc means data=pva median;

class grp_resp grp_amt;

var donor_age income_group wealth_rating;

run;

/*there are some other ways to impute, like cluster imputation using proc fastclus */

/*or using EM, MCMC, Regression to impute in proc mi(see UCLA ATS) */

/*proc mi can also use logistic regression to impute the categorical variables */

/* http://www.ats.ucla.edu/stat/sas/seminars/missing_data/part1.htm */

/* http://www.ats.ucla.edu/stat/sas/seminars/missing_data/part2.htm */

using prxchange in SAS to replace non-alpha and non-numeric characters to blanks

For example, we want to replace the special characters like "# $ ! & * ) ^" in the following words into blanks.

san francisco, @California
Oregon, && U^S
Google, *Mountain View

Function prxchange can be used:

data a;
input old & $100.;
cards;
san francisco, @California
Oregon, && U^S
Google, *Mountain View
;
run;


data a2;
set a;
new=prxchange('s/[^a-zA-Z0-9]/ /i', -1, old);
run;

proc print data=a2;
run;

Number of obs for multi-join in proc sql

It's necessary to remember that proc sql first do the Carditian product and then select obs based on the conditions.

For example, the following code: a.a 1 has 3 replicates, 2 has 2 replicates, b.a is unique, c.a 2 has 3 replicates and 3 has 3 replicates, so totally 1 will have 3*1*1=3 replicates, 2 has 2*1*3=6 replicates, 3 will have 1*1*3=3 replicates. if no other conditions, then totally there will be 3+6+3=12 observations.

data a;
input a b;
cards;
1 2
1 3
1 5
2 3
2 7
3 5
;
run;

data b;
input a c;
cards;
1 2
2 3
3 5
;
run;

data c;
input a d;
cards;
1 1
2 1
2 2
2 3
3 1
3 2
3 3
;
run;


proc sql;
create table test as
select a.*, b.c, c.d
from a,b,c
where a.a=b.a and b.a=c.a;
quit;

proc print data=test;
run;

The output is:

Obs    a    b    c    d

 1    1    2    2    1
 2    1    3    2    1
 3    1    5    2    1
 4    2    3    3    1
 5    2    3    3    3
 6    2    3    3    2
 7    2    7    3    1
 8    2    7    3    3
 9    2    7    3    2
10    3    5    5    1
11    3    5    5    3
12    3    5    5    2

HAVING and CALCULATED in PROC SQL

HAVING: The distinction between HAVING and the WHERE clause is that HAVING conditions can reference summary statistics and are evaluated after aggregations are performed. Thus they take effect “downstream,” on the output side of the process. So we can use WHERE statement in the data sets to replace this HAVING, as is done below.

CALCULATED: if you want to use the calculated result during the SQL select clause, you need to add CALCULATED otherwise there will be error like: ERROR: The following columns were not found in the contributing tables: price_range

In the following the first one is faster than the second, although they generate the same results.

proc sql;
       create table outdata.event_trends1 as
       select event_id, active_max_price-active_min_price as price_range, max(calculated price_range) as max_range, datepart(trend_time) as trend_date
       from fansnapd.event_trends
       group by event_id
       having calculated price_range=max_range
       order by event_id;
quit;

proc sql;
       create table outdata.event_trends(where=(price_range=max_range)) as
       select event_id, active_max_price-active_min_price as price_range, max(calculated price_range) as max_range, datepart(trend_time) as trend_date
       from fansnapd.event_trends
       group by event_id
       order by event_id;
quit;

Variables order in the KEEP RENAME WHERE INDEX clause in the data sets

Variables order in the KEEP RENAME WHERE INDEX clause in the data sets:

1) KEEP = variables name exactly the same as original data set

2) RENAME change the original variable name to the new one

3) WHERE INDEX should use the new name

In the following clicks table has variable event_id guid is_cpa.

proc sort data= clicks (keep=event_id guid is_cpa rename=(event_id=id is_cpa=cpa guid=fsid)) where=(cpa=1) out=test1(index=(fsid));

by id;

run;

SAS study notes: index in sas: when to use index, how to create index, how to modify index / datasets / variables without losing index

Index

When to use index:

1: yield faster access to small subsets of obs for WHERE processing

2: return obs in sorted order for BY processing

3: perform table lookup operations

4: join obs

5: modify obs

Simple index: (index=(myindex)) or (index=(lastname firstname))

Composite index: (index=(comp_index=(lastname firstname)))

Pcoc datasets library=libref;

Modify sas_data_set;

Index delete index_name;

Index create index_specification;

Quit;

Proc sql;

Create index index_name on table_name;

Drop index index_name from table_name;

Quit;

Task	Effect
Add observation(s) to data set	Value/identifier pairs are added to index(es).
Delete observation(s) from data set	Value/identifier pairs are deleted from index(es).
Update observation(s) in data set	Value/identifier pairs are updated in index(es).
Delete data set	The index file is deleted.
Rebuild data set with DATA step	The index file is deleted.
Sort the data in place with the FORCE option in PROC SORT	The index file is deleted.

Use PROC DATASETS to copy, rename(change) datasets, rename variables

DATA SAS-data-file-name (INDEX=

(index-specification-1</UNIQUE><...index-specification-n</UNIQUE>>));

SET SAS-data-set-name ;

RUN;

PROC DATASETS LIBRARY=libref ;

MODIFY SAS-data-set-name;

INDEX DELETE index-name ;

INDEX CREATE index-specification;

QUIT;

PROC SQL;

CREATE INDEX index-name

ON table-name(column-name-1<...,column-name-n>);

DROP INDEX index-name FROM table-name ;

QUIT;

PROC CONTENTS DATA=SAS-data-set-name;

RUN;

PROC DATASETS <LIBRARY=libref> ;

CONTENTS DATA=SAS-data-set-name;

QUIT;

PROC DATASETS LIBRARY=old-libref ;

COPY OUT=new-libref;

SELECT SAS-data-set-name;

QUIT;

PROC COPY OUT=new-libref IN=old-libref ;

SELECT SAS-data-set-name(s);

RUN;

QUIT;

PROC DATASETS LIBRARY=libref ;

CHANGE old-data-set-name = new-data-set-name;

QUIT;

PROC DATASETS LIBRARY=libref ;

MODIFY SAS-data-set-name;

RENAME old-var-name-1 = new-var-name-1 ;

<...old-var-name-n = new-var-name-n> ;

QUIT;

R vs SAS 1: R aggregate v.s. SAS proc summary

In SAS, it's convenient to calculate mean/sum alike statistics over different subset of the original data using proc summary.

In R we can get the similar result using function "aggregate", or use "tapply" for simple condition.

Example:

library(stats)
aggregate(cbind(ncases, ncontrols) ~ alcgp + tobgp, data = esoph, sum)-> data1
aggregate(cbind(ncases, ncontrols) ~ alcgp , data = esoph, sum)-> data2
merge(data1, data2, by.x="alcgp", by.y="alcgp")

gives us:



We can get this from SAS:

data a;
     infile "./esoph.txt" firstobs=2;
     input agegp $ alcgp $ tobgp $ ncases ncontrols;
run;

proc print data=a;
run;


proc summary data=a ;
     class alcgp tobgp;
     var ncases ncontrols;
     output out=temp(drop=_freq_) sum=;
run;

proc print data=temp;
run;

proc sort data=temp;
     by alcgp;
run;

data final(drop=_type_);
     merge temp(where=(_type_=3)) temp(where=(_type_=2) rename=(ncases=tot_ncases ncontrols=tot_ncontrols));
     by alcgp;
run;

proc print data=final;
run;

The output is:

proc surveyselect method=pps: repeatedly sample if sample size is greater than ops maximum available sample size

*** when use pps in proc surveyselect, there is the limitation of maximum sample size (say 1000) every time. If your wanted sample size (say 3000) ***;

*** is greater than this max available sample size, it will show error message. One way to solve this is to sample several times (3 times) with ***;

*** each time to sample the max available sample size. That is, sample 3 times without replacement with each time 1000 samples ***;

%macro ppssample(dataname, m_tot_sample_size, sizevar);
proc sql;
 select floor(sum(&sizevar)/max(&sizevar)) into: m_max_sample_size from &dataname;
quit;

data _null_;
       pps_loop=ceil(&m_tot_sample_size/&m_max_sample_size);
       pps_each_size=ceil(&m_tot_sample_size/pps_loop);
       call symput('pps_loop',pps_loop);
       call symput('pps_each_size', pps_each_size);
run;

%put m_max_sample_size=&m_max_sample_size;

data first_sample;
 set &dataname;
run;

%if &pps_loop>1 %then %do;

 %do i=1 %to &pps_loop;

  proc surveyselect data=first_sample  method=pps sampsize=&pps_each_size
   out=pps_sample_&i. seed=-1;
   size &sizevar;
  run;
  title "to compare: first_sample";
  proc contents data=first_sample;
  run;
  title "to coampare: pps_sample_&i.";
  proc contents data=pps_sample_&i.;
  run;

  proc sql;
   create table first_sample as
   select * from first_sample
   EXCEPT CORR
   select * from pps_sample_&i.(drop=samplingweight SelectionProb);
  quit;

                      title "contents of first_sample for loop &i";
                      proc contents data=first_sample;
                      run;

  data pps_sample;
   %if &i=1 %then %do;
    set pps_sample_1;
   %end;
   %else %do;
    set pps_sample pps_sample_&i.;
   %end;
  run;
              %end;

      %end;

      %else %do;

              proc surveyselect data=first_sample method=pps sampsize=&m_tot_sample_size out=pps_sample seed=-1;
     size &sizevar;
  run;

      %end;

proc surveyselect data=pps_sample method=srs n=&m_tot_sample_size out=final;
run;

%mend ppssample;

%ppssample(prep,5000,revenue);

function fileexist to determine a file exist or not

in linux, the filename is case sensitive and remember to add surfix .

data a;
a=fileexist("/data02/temp/temp_hsong/seo/seo_expr/seo_visitnew_raw_incr.sas7bdat");
run;

proc print data=a;
run;


Obs a

1 1

IS MISSING / IS NULL to select missing rows

In sas, if want to select those missing rows, you can use VAR=. (for numeric var) or VAR='' OR VAR=' ' (for char var).

A convenient way is to use VAR IS MISSING   OR   VAR IS NULL.

Using Multiple SET Statements

When you use multiple SET statements,

• processing stops when SAS encounters the end-of-file (EOF) marker on either data set (even if there is more data in the other data set) 
•  
• the variables in the program data vector (PDV) are not reinitialized when a second SET statement is executed. 

This is useful when want to combine the summary info with the details rows of the data set.

suppose sasuser.summary has one var named sum, with one obs valued as 100; sasuser.monthsum has two vars: month sale:
month sale
1      28
2      32
3      40
we want to combine them together to calculate the percentage for each month.

data a; 

  if _N_=1 then set sasuser.summary;

  set sasuser.monthsum;

  pct=sale/sum;

run;

the automatic variable _N_ keeps track of how many times the DATA step has begun to execute. The following DATA step uses _N_ to keep SAS from reaching the EOF marker for Sasuser.Summary after the first iteration of the step. Since the variables in the PDV will not be reinitialized on each iteration, the first value of Summary.Cargosum will be retained in the PDV for each observation that is read from Sasuser.Monthsum. 

zz: SAS与R语言的数据加载与转化

第一：R加载/调用SAS
--在SAS中生成传送文件
LIBNAME SAS_R xport 'C:\sea.xpt';
DATA SAS_R.sea;
SET custdet1;
RUN;
--在R中读入
library(foreign)
library(Hmisc)
sea<-sasxport.get("c:/sea.xpt")
head(mydata)

第二：SAS调用R的数据
Library(foreign)
write.foreign(sea,"c:/sea.txt","c:/sea.sas",package="SAS")
在C盘中会生成二个文件：
SAS程序代码
一个是TXT文件
然后在SAS中直接加载程序就行。

zz: 如何在R中调用matlab

install.packages("R.matlab")
library(R.matlab)
path <- system.file("mat-files", package="R.matlab")
mat <- readMat(file.path(path, "structLooped.mat"))
s <- mat$s
fields <- dimnames(s)[[1]]
cat("Field names: ", paste(fields, collapse=", "), "\n", sep="");

print(s)

check difference of PROC SUMMARY <> PROC MEANS (not finished)

***  if VAR is missing, sas will drop that records when calculating mean. To include that records, we can impute that missing var by 0 ***;

data test;
  input a $ 1-3 b 4-5 ;
  cards;
  a 1
  a
  a 3
  b
  b 2
  b 1
  c 1
  c
    1
    4
  ;
run;

proc print data=test noobs;
  title "print of original data";
run;

proc summary data=test nway missing;
  var b;
  class a;
  output out=sum1(drop=_type_ _freq_) mean=;
run;

proc print data=sum1;
  title "Summary with MISSING opinion";
run;

proc summary data=test nway;
  var b;
  class a;
  output out=sum2(drop=_type_ _freq_) mean=;
run;

proc print data=sum2;
  title "Summary without MISSING opinion";
run;

proc means data=test nway missing;
  var b;
  class a;
  output out=means1(drop=_type_ _freq_) mean=;
run;

proc print data=means1;
  title "Means with MISSING opinion";
run;

proc summary data=test nway;
  var b;
  class a;
  output out=means2(drop=_type_ _freq_) mean=;
run;

proc print data=means2;
  title "Means without MISSING opinion";
run;

^LSummary with MISSING opinion

Obs    a     b

 1          2.5
 2     a    2.0
 3     b    1.5
 4     c    1.0


^LSummary without MISSING opinion

Obs    a     b

 1     a    2.0
 2     b    1.5
 3     c    1.0


^LMeans with MISSING opinion

Obs    a     b

 1          2.5
 2     a    2.0
 3     b    1.5
 4     c    1.0


^LMeans without MISSING opinion

Obs    a     b

 1     a    2.0
 2     b    1.5
 3     c    1.0

zz: proc means and proc mixed for paired t test

/*********************************************************************
FILENAME: MIXVMEAN.SAS
SUBJECT HEADING: STAT
INITIALS:  KBW
DATE:  7/26/96
PROGRAM:  SAS
VERSION:  6.11
PLATFORM:  WINDOWS 3.11 TS040
TITLE:  USING PROC MIXED FOR A PAIRED T-TEST COMPARED TO PROC MEANS

DESCRIPTION:  THIS PROGRAM DOES A PAIRED T-TEST, FOR THE SAME DATA,
              FIRST USING PROC MEANS, AND THEN USING PROC MIXED.
              NOTE THAT THE DATA STRUCTURE FOR THESE TWO METHODS
              IS DIFFERENT.  THE DATA STEPS AT THE BEGINNING SET
              UP THE TWO DIFFERENT DATA STRUCTURES.
              NOTE THAT THE T-TEST, THE DEGREES OF FREEDOM AND
              THE P-VALUES ARE THE SAME FOR BOTH METHODS.
**********************************************************************/


data test;
  input y1 y2;
  diff=y1-y2;
  Pair+1;
  datalines;
  13 15
  12 14
  17 17.2
  14 18
  11 12
  5  4.1
  7 9.3
     ;
run;

proc print data=test;
  title 'printout of original data set';
run;

proc means n mean t prt;
  var diff;
  title 'paired t-test using proc means';
run;

data test2(keep=y pair group);
  set test;
     y=y1;
     group=1;
     output;
     y=y2;
     group=0;
     output;
run;

proc print data=test2;
  title 'rearranged data for proc mixed';
run;

proc mixed;
    class pair group;
    model y=group;
    random pair;
    lsmeans group / pdiff;
    title 'paired t-test using proc mixed';
run; 

from SAS: Determining the Number of Variables and Observations in a Data Set

%macro obsnvars(ds);
   %global dset nvars nobs;
   %let dset=&ds;
   %let dsid = %sysfunc(open(&dset));
   %if &dsid %then
      %do;
         %let nobs =%sysfunc(attrn(&dsid,NOBS));
         %let nvars=%sysfunc(attrn(&dsid,NVARS));
         %let rc = %sysfunc(close(&dsid));
         %put &dset has &nvars  variable(s) and &nobs observation(s).;
      %end;
   %else
      %put Open for data set &dset failed - %sysfunc(sysmsg());
%mend obsnvars;

%obsnvars(sasuser.houses)

******************************************************************************;

%macro nobs(Dsn= /*Data set name */);
  if exist("&Dsn") then do;
    Dsid = open("&Dsn","i");
    Nobs = attrn(Dsid,"Nlobs");
  end;
  else Nobs=.;
  rc = close(Dsid);
%mend nobs;