Dear @Rick_SAS Thanks for your prompt reply. The complete macro is as below: %macro gwnbr(data=,y=,x=,lat=,long=,h=,grid=,latg=,longg=,gwr=,method=,alphag=,offset=,geocod=,out=);
proc iml;
use &data;
read all var {&y} into y;
read all var {&x} into x;
read all var{&long &lat} into COORD;
n=nrow(y);
%if &offset= %then %do; offset=j(n,1,0); %end;
%else %do; read all var {&offset} into offset; %end;
%if &grid= %then %do;
read all var{&long &lat} into POINTS;
read all var{&geocod} into geocod_;
%end;
close &data;
%if &grid^= %then %do;
use &grid;
read all var{&longg &latg} into POINTS;
close &grid;
geocod_=nrow(points,1,0);
%end;
x=j(n,1,1)||x;
yhat=j(n,1,0);
h=&h;
gwr="&gwr"; *global,local, poisson;
method="&method"; *fixed, adaptive1, adaptiven;
m=nrow(POINTS);
bii=j(ncol(x)*m,2,0); alphaii= j(m,2,0);
xcoord=j(ncol(x)*m,1,0); ycoord=j(ncol(x)*m,1,0);
&geocod= j(ncol(x)*m,1,0);
sebi=j(ncol(x)*m,1,0); sealphai= j(m,1,0);
S=j(n,n,0);
yp=y-sum(y)/n;
probai=j(m,1,0); probbi=j(m,1,0);
yhat=j(m,1,0);
res= j(m,1,0);
if gwr^="poisson" then do;
ym=sum(y)/nrow(y);
u=(y+ym)/2;
n=log(u);
par=1; ddpar=1; j=0; aux2=0;
do while (abs(ddpar)>0.00001);
aux1=0;
dpar=1;
parold=par;
do while (abs(dpar)>0.001);
aux1=aux1+1;
if par<0 then par=0.00001;
par=choose(par<1E-10,1E-10,par);
g=sum(digamma(par+y)-digamma(par)+log(par)+1-log(par+u)-(par+y)/(par+u));
hess=sum(trigamma(par+y)-trigamma(par)+1/par-2/(par+u)+(y+par)/((par+u)#(par+u)));
hess=choose(abs(hess)<1E-23,sign(hess)*1E-23,hess); *CONFERIR!!!;
hess=choose(hess=0,1E-23,hess);
par0=par;
par=par0-inv(hess)*g;
if aux1>50 & par>1E5 then do;
dpar= 0.0001;
aux2=aux2+1;
if aux2=1 then par=2 ;
else if aux2=2 then par=1E5;
else if aux2=3 then par=0.0001;
end;
else dpar=par-par0;
end;
a=1/par; dev=0; ddev=1; i=0;
do while (abs(ddev)>0.00001);
i=i+1;
w=(u/(1+a*u))+(y-u)#(a*u/(1+2*a*u+a*a*u#u));
w=choose(w<=0,1E-5,w);
z=n+(y-u)/(w#(1+a*u)) - offset;
b=inv((x#w)`*x)*(x#w)`*z;
n=x*b + offset;
n=choose(n>1E2,1E2,n);
u=exp(n);
olddev=dev;
tt=y/u;
tt=choose(tt=0,1E-10,tt);
dev=2*sum(y#log(tt)-(y+1/a)#log((1+a*y)/(1+a*u)));
ddev=dev-olddev;
end;
if aux2>4 then ddpar=1E-9;
else ddpar=par-parold;
end;
%if &alphag= %then %do; alphag=a;%end;
%else %if &alphag=0 %then %do; alphag=1e-8;%end;
%else %do; alphag=&alphag;%end;
bg=b;
parg=par;
end;
if gwr="global" then print alphag aux2;
n=nrow(y);
aux2=0;
do i=1 to m;
d=j(1,3,0);
do j=1 to n;
if abs(COORD[,1])<180 then do;
dif=abs(POINTS[i,1]-COORD[j,1]);
raio=arcos(-1)/180;
ang=sin(POINTS[i,2]*raio)*sin(COORD[j,2]*raio)+cos(POINTS[i,2]*raio)*cos(COORD[j,2]*raio)*cos(dif*raio);
if round(ang,0.000000001)=1 then arco=0;
else arco=arcos(ang);
d1=arco*6371 /*Earth's Radius = 6371 (approximately)*/;
end;
else d1=sqrt((POINTS[i,1]-COORD[j,1])**2+(POINTS[i,2]-COORD[j,2])**2);
d[1]=i; d[2]=j; d[3]=d1;
if j=1 then dist=d; *cleaning dist where i value changes;
else dist=dist//d;
end;
w=j(n,1,0);
if method= "fixed" then do;
do jj=1 to n;
w[jj]=exp(-0.5*(dist[jj,3]/h)**2);
end;
end;
else if method= "adaptiven" then do;
do jj=1 to n;
if dist[jj,3]<=h then w[jj]=(1-(dist[jj,3]/h)**2)**2;
else w[jj]= 0;
end;
end;
else if method= "adaptive1" then do;
w=j(n,2,0);
call sort(dist,{3});
dist=dist||(1:n)`;
hn=dist[h,3]; *bandwith for the point i;
do jj=1 to n;
if dist[jj,4]<=h then w[jj,1]=(1-(dist[jj,3]/hn)**2)**2;
else w[jj,1]=0;
w[jj,2]=dist[jj,2];
end;
call sort(w,{2});
end;
wi=w[,1];
ym=sum(y)/nrow(y);
uj=(y+ym)/2;
nj=log(uj);
ddpar=1; jj=0; count=0; aux2=0;
if i=1 | aux2=5 | count=4 then par=1; else par=alphaii[i-1,2];
do while (abs(ddpar)>0.000001);
dpar=1;
if ddpar=1 then parold=1.8139;
else parold=par;
aux1=0;
if gwr="global" | gwr="poisson" then do;
dpar=0.00001;
if gwr= "global" then par=1/alphag;
end;
/* computing alpha=1/par, where par=theta=r */
do while (abs(dpar)>0.001);
aux1=aux1+1;
if gwr="local" then do;
par=choose(par<1E-10,1E-10,par);
g=sum((digamma(par+y)-digamma(par)+log(par)+1-log(par+uj)-(par+y)/(par+uj))#w[,1]);
hess=sum((trigamma(par+y)-trigamma(par)+1/par-2/(par+uj)+(y+par)/((par+uj)#(par+uj)))#w[,1]);
end;
par0=par;
hess=choose(abs(hess)<1E-23,sign(hess)*1E-23,hess);
hess=choose(hess=0,1E-23,hess);
par=par0-inv(hess)*g;
if par<=0 then do;
count=count+1;
if count=1 then par=0.000001;
else if count=2 then par=0.0001;
else par=1/alphag;
end;
if aux1>100 & par>1E5 then do; *MAXINTA;
dpar= 0.0001;
if aux2=0 then par=1/alphag + 0.0011;
if aux2=1 then par=2 ;
else if aux2=2 then par=1E5;
else if aux2=3 then par=0.0001;
aux2=aux2+1;
end;
else do;
dpar=par-par0;
if par<1E-3 then dpar=dpar*100;
end;
end;
if gwr= "poisson" then alpha=0;
else alpha=1/par;
dev=0; ddev=1; cont=0;
/* computing beta */
do while (abs(ddev)>0.000001);
cont=cont+1;
Ai=(uj/(1+alpha*uj))+(y-uj)#(alpha*uj/(1+2*alpha*uj+alpha*alpha*uj#uj));
Ai=choose(Ai<=0,1E-5,Ai);
zj=nj+(y-uj)/(Ai#(1+alpha*uj))-offset;
if det(x`*(wi#Ai#x))=0 then bi=j(ncol(x),1,0);
else bi=inv(x`*(wi#Ai#x))*x`*(wi#Ai#zj);
nj=x*bi + offset;
nj=choose(nj>1E2,1E2,nj);
uj=exp(nj);
olddev=dev;
uj=choose(uj<1E-150,1E-150,uj);
tt=y/uj;
tt=choose(tt=0,1E-10,tt);
if gwr= "poisson" then dev=2*sum(y#log(tt)-(y-uj));
else dev=2*sum(y#log(tt)-(y+1/alpha)#log((1+alpha*y)/(1+alpha*uj)));
if cont>100 then ddev= 0.0000001; *MAXINTB;
else ddev=dev-olddev;
end;
jj=jj+1;
*print jj bi;
if gwr="global" | gwr="poisson" | aux2>4 | count>3 | jj>200 then ddpar=1E-9;
else do;
ddpar=par-parold;
if par<1E-3 then ddpar=ddpar*100;
end;
/* print j aux1 cont aux2 count parold par ddpar;*/
end;
if aux2>4 then probai[i]=1;
if count>3 then probai[i]=2;
Ai2=(uj/(1+alpha*uj))+(y-uj)#(alpha*uj/(1+2*alpha*uj+alpha*alpha*uj#uj));
if Ai2[><,]<1E-5 then do;
probbi[i]=1;
Ai2=choose(Ai2<1E-5,1E-5,Ai2);
end;
Ai=Ai2;
%if &grid= | &grid=&data %then %do;
if det(x`*(wi#Ai#x))=0 then S[i,]=j(1,n,0);
else S[i,]= x[i,]*inv(x`*(wi#Ai#x))*(x#wi#Ai)`;
%end;
C=inv(x`*(wi#Ai#x));
varb= C;
seb=sqrt(vecdiag(varb));
if gwr^="poisson" then do;
ser=sqrt(1/abs(hess));
r=1/alpha;
sealpha=ser/(r**2);
sealphai[i,1]=sealpha;
alphaii[i,1]=i;
alphaii[i,2]= alpha;
end;
m1=(i-1)*ncol(x)+1;
m2=m1+(ncol(x)-1);
sebi[m1:m2,1]=seb;
bii[m1:m2,1]=i;
bii[m1:m2,2]=bi;
xcoord[m1:m2,1]= POINTS[i,1];
ycoord[m1:m2,1]= POINTS[i,2];
&geocod[m1:m2,1]= geocod_[i,1];
%if &grid= | &grid=&data %then %do;
yhat[i]=uj[i];
%end;
end;
tstat= bii[,2]/sebi;
probtstat=2*(1-probnorm(abs(tstat)));
if gwr^="poisson" then do;
atstat= alphaii[,2]/sealphai;
aprobtstat=2*(1-probnorm(abs(atstat))); *check for normality;
end;
else do;
atstat=j(n,1,0);
aprobtstat=j(n,1,1);
end;
b=bii[,2];
alphai=alphaii[,2];
_id_= bii[,1];
_ida_=alphaii[,1];
_beta_=shape(bii[,1:2],n);
i=do(2,ncol(_beta_),2);
_beta_=_beta_[,i];
call qntl(qntl,_beta_);
qntl=qntl//(qntl[3,]-qntl[1,]);
descriptb=_beta_[:,]//_beta_[><,]//_beta_[<>,];
print qntl[label="Quantiles of GWNBR Parameter Estimates"
rowname={"P25", "P50", "P75", "IQR"} colname={'Intercept' &x}],,
descriptb[label="Descriptive Statistics" rowname={"Mean", "Min", "Max"}
colname={'Intercept' &x}];
_stdbeta_=shape(sebi,n);
call qntl(qntls,_stdbeta_);
qntls=qntls//(qntls[3,]-qntls[1,]);
descripts=_stdbeta_[:,]//_stdbeta_[><,]//_stdbeta_[<>,];
print qntls[label="Quantiles of GWNBR Standard Errors"
rowname={"P25", "P50", "P75", "IQR"} colname={'Intercept' &x}],,
descripts[label="Descriptive Statistics of Standard Errors" rowname={"Mean", "Min", "Max"}
colname={'Intercept' &x}];
%if &grid= | &grid=&data %then %do;
yhat=choose(yhat<1E-150,1E-150,yhat);
tt=y/yhat;
tt=choose(tt=0,1E-10,tt);
if gwr= "poisson" then do;
dev=2*sum(y#log(tt)-(y-yhat));
tt2=y/y[:];tt2=choose(tt2=0,1E-10,tt2);
devnull=2*sum(y#log(tt2)-(y-y[:]));
pctdev=1-dev/devnull;
end;
else do;
dev=2*sum(y#log(tt)-(y+1/alphai)#log((1+alphai#y)/(1+alphai#yhat)));
tt2=y/y[:];tt2=choose(tt2=0,1E-10,tt2);
devnull=2*sum(y#log(tt2)-(y+1/alphai)#log((1+alphai#y)/(1+alphai#y[:])));
pctdev=1-dev/devnull;
end;
if gwr^="poisson" then do;
a2=y+1/alphai; b2=1/alphai;
algamma=j(n,1,0); blgamma=j(n,1,0);
do i=1 to nrow(y);
algamma[i]=lgamma(a2[i]);
blgamma[i]=lgamma(b2[i]);
end;
end;
c2=y+1;
clgamma=j(n,1,0);
do i=1 to nrow(y);
clgamma[i]=lgamma(c2[i]);
end;
if gwr^="poisson" then do;
ll=sum(y#log(alphai#yhat)-(y+1/alphai)#log(1+alphai#yhat)+ algamma - blgamma - clgamma );
if gwr="global" & alphai^=1/parg then npar=trace(S);
else npar=trace(S)+1;
tt=y/(alphai#yhat);tt=choose(tt=0,1E-10,tt);
ll1=sum(y#log(tt)-y+(y+1/alphai)#log(1+alphai#yhat)-algamma+blgamma);
tt2=y/y[:];tt2=choose(tt2=0,1E-10,tt2);
llnull=sum(y#log(tt2));
pctll=1-ll1/llnull;
end;
else do;
ll=sum(-yhat+y#log(yhat)-clgamma);
npar=trace(S);
pctll=pctdev;
end;
adjpctdev=1-((nrow(y)-1)/(nrow(y)-npar))*(1-pctdev);
adjpctll=1-((nrow(y)-1)/(nrow(y)-npar))*(1-pctll);
resord=y-yhat;
sigma2= (resord`*resord)/(n-npar);
sii=vecdiag(S);
res=resord/sqrt(sigma2#(1-sii));
res=unique(_id_)`||COORD[,1]||COORD[,2]||y||yhat||res||resord;
/*AIC= 2*npar + dev;*/
AIC= 2*npar - 2*ll;
AICC= AIC +(2*npar*(npar+1))/(n-npar-1);
BIC= npar*log(n) - 2*ll ;
_malpha_=0.05*(ncol(x)/npar);
_t_critical_=abs(tinv(_malpha_/2,n-npar));
print _malpha_[label="alpha-level=0.05"] _t_critical_[format=comma6.2 label="t-Critical"] npar;
print gwr method ll dev pctdev adjpctdev pctll adjpctll npar aic aicc bic;
create _res_ from res[colname={"_id_" "xcoord" "ycoord" "yobs" "yhat" "res" "resraw"}];
append from res;
stat=ll|| dev|| pctdev || adjpctdev|| pctll || adjpctll || npar|| aic|| aicc|| bic;
create _stat_ from stat[colname={"l1" "dev" "pctdev" "adjpctdev" "pctll" "adjpctll" "npar" "aic" "aicc" "bic"}];
append from stat;
%end;
%else %do; print gwr method; %end;
create _beta_ var{_id_ &geocod xcoord ycoord b sebi tstat probtstat}; * _beta_ has beta vector for each point i;
append;
xcoord=COORD[,1];ycoord=COORD[,2];
&geocod=unique(&geocod)`;
sig_alpha=j(n,1,"not significant at 90%");
v1=npar;
do i=1 to n;
if aprobtstat[i]<0.01*(ncol(x)/v1) then sig_alpha[i]="significant at 95%";
else if aprobtstat[i]<0.1*(ncol(x)/v1) then sig_alpha[i]="significant at 90%";
else sig_alpha[i]="not significant at 90%";
end;
create _alpha_ var{_ida_ &geocod xcoord ycoord alphai sealphai atstat aprobtstat sig_alpha probai probbi}; * _alpha_ has alpha vector for each point i;
append;
_tstat_=_beta_/_stdbeta_;
_probt_=2*(1-probnorm(abs(_tstat_)));
_bistdt_=geocod_||COORD||_beta_||_stdbeta_||_tstat_||_probt_;
_colname1_={"Intercept" &x};
_label_=repeat("std_",ncol(x))//repeat("tstat_",ncol(x))//repeat("probt_",ncol(x));
_colname_={"&geocod" "x" "y"}||_colname1_||concat(_label_,repeat(_colname1_`,3))`;
call change(_colname_, "_ ", "_");
call change(_colname_, "_ ", "_");
create _parameters_ from _bistdt_[colname=_colname_];
append from _bistdt_;
close _parameters_;
_sig_=j(n,ncol(x),"not significant at 90%");
v1=npar;
do i=1 to n;
do j=1 to ncol(x);
if _probt_[i,j]<0.01*(ncol(x)/v1) then _sig_[i,j]="significant at 99%";
else if _probt_[i,j]<0.05*(ncol(x)/v1) then _sig_[i,j]="significant at 95%";
else if _probt_[i,j]<0.1*(ncol(x)/v1) then _sig_[i,j]="significant at 90%";
else _sig_[i,j]="not significant at 90%";
end;
end;
_colname1_={"Intercept" &x};
_label_=repeat("sig_",ncol(x));
_colname_=concat(_label_,repeat(_colname1_`,1))`;
create _sig_parameters2_ from _sig_[colname=_colname_];
append from _sig_;
/*
%let nvar=0;
%do %while(%scan(%str(&x),&nvar+1)~=);
%let nvar=%eval(&nvar+1);
%end;
use _beta_;
read all into b;
close _beta_;
n=nrow(b);
npar=&nvar+1;
%do i=0 %to &nvar;
b&i=j(1,8,0);
nome={"_id_" "&geocod" "xcoord" "ycoord" "b" "sebi" "tstat" "probtstat"};
create &out._b&i from b&i[colname=nome];
do i=1 to (n/npar);
b&i[1,]=b[(i-1)*npar+&i+1,];
append from b&i;
end;
%end;
*/
quit;
%mend gwnbr; The Do loop functions is to generate weighting matrix for each data points. Since I have 227 data points, the first weighting matrix contains weights of point 1 to point 1, point 1 to point 2......, point 1 to point 227, then the next matrix contains weights of point 2 to point 1, point 2 to point 2,......, point 2 to point 227, and the last matrix will be the weights of point 227 to point 1, point 227 to point 2,......point 227 to point 227. These weights are based on the distances between data points, and the SAS/IML program use the coordinates of each data points to calculate the distances. I was wondering after I calculate the weights like the way you suggested, how can I modify the code and let the program use the correct weighting matrix in each iterations like the original program did(the first should be the weights between point 1 and other points, and the second one are the weights between point 2 and other points, the last one are the weights between point 227 and other points). I'd only like to conduct GWNBR under gwr=global and method=fixed as the snippet in my first post. Attached are the gwnbr macro, my data to conduct gwnbr and the distance table for each data points. Thank you for your time and help.
... View more
