BookmarkSubscribeRSS Feed
Henry1101
Calcite | Level 5

Hello, I'm using the GWNBR macro that developed by Alan. The distances between the data points are calculated by the coordinates of each data point and the weights are calculated by the distances in the original code below:

 

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;

For some reasons, I'd like to use the distances that created by myself to calculate the weights and conduct GWNBR. It seems that I should replace the step of creating the distances matrix from the coordinates with the reading of my distance table, but I'm not sure how to modify the code above. Is there a way to do this? 

Attached is the distances table I created, which contains 227 distances matrices(since I have 227 data points).

 

Thank you all in advance,
Henry

 

5 REPLIES 5
Rick_SAS
SAS Super FREQ

I am not familiar with the macro, and you've only provided a snippet of a SAS/IML program. If you want to skip the first portion, which computes dist from POINTS and COORDS, then you can just calculate the weights (for method="fixed") as

 

w = exp(-0.5*(d/h)##2);

 

Here d is the third column in your Excel file. You'll have to define the bandwidth h somehow. It must be defined prior to the snippet that you've posted.  I'm not sure how the DO I = 1 to m loop gets used.

Henry1101
Calcite | Level 5

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.

 



 

Rick_SAS
SAS Super FREQ

> how can I modify the code and let the program use the correct weighting matrix in each iterations like the original program did

I'd only like to conduct GWNBR under gwr=global and method=fixed 

 

Unfortunately, I do not have the time to study the paper and help you decipher a 400-line macro. I suggest

1. Start with a simple data set that has 5 points.

2. Extract the program in the macro and just work with the SAS/IML program and some %LET statements

3. Get rid of the program sections that are not related to gwr=global and method=fixed 

 

Good luck.

Henry1101
Calcite | Level 5

Dewa @Rick_SAS 

Thank you for your response and help. I will follow your suggestions and try to figure it out.

 

I'd like to ask a question about SAS/IML function. After reading the excel file "revised_dist" in SAS, I'd like to generate the matrix according to variable "idi" in the data set. In other words, I want to  DO idi = 1 to 227, and generate the matrix that contains the rows for the value of "idi" equals to 1, and the second matrix contains the rows of "idi" equals to 2, then the last matrix contains the rows of "idi" equals to 227. Are there any statements or functions can do this? Thank you for the help.

idi1.PNGimage.pngimage.png

Rick_SAS
SAS Super FREQ

The following statements read the data into a 227 x 227 matrix:

 

proc iml;
use revised_dist;
read all var "d1";
close;
D = shape( d1, sqrt(nrow(d1)) );  /* make into nxn matrix */
print (dimension(D));  /* should print  227 x 227 */

If you want to access the values where id1=1, use

r = D[1, ]; /* extract the first row */

if you want the 227th row, use 

r = D[227, ]; /* extract the last row */

 

BTW, you might want to contact the macro author. He is a nice guy and might be able to help you modify the macro to your use-case.

SAS Innovate 2025: Save the Date

 SAS Innovate 2025 is scheduled for May 6-9 in Orlando, FL. Sign up to be first to learn about the agenda and registration!

Save the date!

How to Concatenate Values

Learn how use the CAT functions in SAS to join values from multiple variables into a single value.

Find more tutorials on the SAS Users YouTube channel.

SAS Training: Just a Click Away

 Ready to level-up your skills? Choose your own adventure.

Browse our catalog!

Discussion stats
  • 5 replies
  • 1839 views
  • 1 like
  • 2 in conversation