2011-09-24

SAS main macro for GWGLM

The following program should work well with SAS version 9.2 TS2M2 or higher. Users
are advised to check their SAS version prior to the use of the program.


/*The MACRO for conducting GWGLM*/
%macro gwglm(data=,out=,h=,kernel=gaussian,distrb=normal,response=,vars=,offset=);
%global nvar;
%let nvar=0;
%do %while(%scan(%str(&vars),&nvar+1)~=);
    %let nvar=%eval(&nvar+1);
    %let var&nvar=%scan(%str(&vars),&nvar);
%end;

%global nobs;
proc contents data=&data out=temp noprint;
data _null_; set temp; call symput ('nobs',compress(nobs)); run;

%do j=1 %to &nobs;
data use1; set &data;
   idnum=symget('j');
   if _n_=idnum then do;
      call symput ('xstart', X);
      call symput ('ystart', Y);
   end;
run;
data use1; set use1;
   d=sqrt((X-&xstart)**2+(Y-&ystart)**2);
proc sort; by d;
run;

data aa; set use1;
   if _n_=&h;
   call symput ('b', d);
run;
data use1; set use1;
%if &kernel=gaussian %then %do;
       w=exp(-0.5*((d/&b)**2)); /*gaussian kernel*/
%end;
%else %if &kernel=bisquare %then %do; /*bisquare nearest neighbor*/
   if _n_ <=&h then do;
       w=(1-(d/&b)**2)**2;
   end;
   else do; w=0; end;
%end;
%else %if &kernel=exponential %then %do;  /*exponential kernel*/
     w=exp(-d/&b);
%end;
run;

ods noresults;
ods select none;
proc GENMOD data=use1 %if &distrb=binomial %then %do; DESCENDING %end;;
    model &response=&vars
          / expected scoring=100
          %if &distrb=normal %then %do;
                dist=normal;
          %end;
          %if &distrb=poisson %then %do;
                 OFFSET=&offset dist=poisson;
          %end;
          %if &distrb=binomial %then %do;
                 dist=binomial;
          %end;
    output out=p&j p=muhat xbeta=xbeta leverage=hh;
    ods output ParameterEstimates=est&j;
    weight w;
run;

data est&j; set est&j; idnum=&j;
data a0; set est&j;
     if Parameter='Intercept';
        rename Estimate=Est_int StdErr=str_int LowerWaldCL=lcl_int UpperWaldCL=ucl_int ChiSq=chisq_int ProbChiSq=p_int;
     drop Parameter DF;
     attrib _all_ label=" ";
run;

%do i=1 %to &nvar;
data a&i; set est&j;
    if Parameter="&&var&i.";
    rename Estimate=Est_&&var&i StdErr=str_&&var&i LowerWaldCL=lcl_&&var&i UpperWaldCL=ucl_&&var&i ChiSq=chisq_&&var&i ProbChiSq=p_&&var&i;
    drop Parameter DF;
    attrib _all_ label=" ";
run;
%end;
data tt&j;
   merge %do i=0 %to &nvar; a&i %end;;
   by idnum;
run;

data pp&j; set p&j; if _n_=1; run;
%end;

ods results;
ods select all;

/*Five number summary of the estimates*/
data finalbeta; set %do j=1 %to &nobs; tt&j %end;; run;
proc means data=finalbeta QMETHOD=os min q1 median q3 max;
   var Est_int %do i=1 %to &nvar; Est_&&var&i. %end;;
run;

data finalp; set %do j=1 %to &nobs; pp&j %end;;
   %if &distrb=normal %then %do;
       deviancei=(&response-muhat)**2;
   %end;
   %if &distrb=poisson %then %do;
          %if &response=0 %then %do ;
                   deviancei=-2*(-muhat);
          %end;
          %if &response>0 %then %do ;
                   deviancei=2*(&response*log(&response/muhat)-(&response-muhat));
          %end;

    %end;
   %if &distrb=binomial %then %do;
       deviancei=(-2)*(&response*log(muhat) + (1-&response)*log(1-muhat));
   %end;
   llri=(-1)*deviancei/2;
   aici=deviancei+2*hh;
   bici=deviancei+hh*log(&nobs);
run;

**Model fitting Statistics;
proc sql;
    create table dd as
    select sum(llri) as LLR,
           sum(deviancei) as D label='Deviance',
           sum(hh) as K label='Effective number of parameters',
           &nobs- calculated K as EDF label='Effective degrees of freedom',
           %if &distrb=normal %then %do;
           sqrt(calculated D/(&nobs- calculated K)) as sigma,
           sqrt(calculated D/&nobs) as sigma_ML,
           2*(&nobs)*log(sqrt(calculated D/&nobs))+(&nobs)*log(2*CONSTANT('PI'))+&nobs+ calculated K as AIC,
           calculated AIC - &nobs - calculated K +(&nobs)*((&nobs+ calculated K)/(&nobs-2- calculated K)) as AICC,
           %end;
           %else %do;
           sum(aici) as AIC,
           calculated D+2* calculated K*(&nobs/(&nobs- calculated K-1)) as AICC,
           %end;
           sum(bici) as BIC          
    from finalp;
quit;
proc print data=dd label;

data cc; set &data; idnum+1; keep idnum X Y; run;
data &out; merge cc finalbeta; by idnum; run;
%mend gwglm;

沒有留言:

張貼留言