2012-02-10

Program update note

According to the feedback provided by some users, there exist a copule of bugs in the GWGLM program. We have fixed them; see the blue texts for the previous post created on Sep. 24, 2011.
We appreciate any comments from you and hope to better our proposed SAS macros for public use.

2011-09-25

Link to download all datasets used in CMPB paper

The GWGLM macro and related sas datasets used in our paper accepted by Computer Methods and Programs in Biomedicine (this paper is available at: http://dx.doi.org/10.1016/j.cmpb.2011.10.006) can be download via the website. Please go to the link , click the files you desired and then select "download".

Suggested citation when using the macro in your research: 
V.Y.-J. Chen, T.-C. Yang, SAS macro programs for geographically weighted generalized linear modeling with spatial point data: Applications to health research, Comput. Methods Programs Biomed. (2011), doi:10.1016/j.cmpb.2011.10.006

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;