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;
沒有留言:
張貼留言