Skip to content

Instantly share code, notes, and snippets.

@cjdinger

cjdinger/surveygenmod.sas

Last active August 8, 2022 21:32
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save cjdinger/022cbec1b827a2265d4aab92c6e85f1c to your computer and use it in GitHub Desktop.
Save cjdinger/022cbec1b827a2265d4aab92c6e85f1c to your computer and use it in GitHub Desktop.
Download ZIP
Raw
surveygenmod.sas
/*
%SURVEYGENMOD macro copied from https://support.sas.com/resources/papers/proceedings17/0268-2017.pdf
%SURVEYGENMOD Macro: An Alternative to Deal with Complex Survey Design for the GENMOD Procedure
Alan Ricardo da Silva, Universidade de Brasilia, Dep. de Estatística, Brazil
Original PDF was not conducive for copy/paste, so this version attempts to correct
for many end-of-line wrapping problems. Might not have caught them all!
*/
%macro surveygenmod(data=, y=, x=, offset=, weight=,
strata=, cluster=, domain=,fpc=, dist=normal,
link=identity, xzip=, intercept=y, scale=deviance,
vadjust=n, alpha=0.05, delta=1, fi=0.5,
maxiter=100, eps=0.000001);
/********************************************************* *
The distributions available in %surveygenmod macro are:
* NORMAL, GAMMA, IG (INVERSE GAUSSIAN), POISSON, NEGBIN (NEGATIVE BINOMIAL)
* BINOMIAL (AND MULTINOMIAL), ZIP AND ZINB.
********************************************************** *
* The link functions available in %surveyglm macro are:
* IDENTITY, INV (INVERSE), INV2 (INVERSE SQUARED), LOG, LOGIT (GENERALIZED LOGIT).
**********************************************************/
%if &data= %then
%do;
%put ERROR: The database must be informed;
%end;
%else %if &y= %then
%do;
%put ERROR: The response variable must be informed;
%end;
%else %if &x= %then
%do;
%put ERROR: The covariates must be informed;
%end;
%else %if (%upcase(&dist)=NORMAL | %upcase(&dist)=GAMMA | %upcase(&dist)=IG) & &scale= %then
%do;
%put ERROR: The scale estimation method (DEVIANCE or PEARSON method) must be informed;
%end;
%else
%do;
/*Sorting data*/
%if %upcase(&strata) NE & %upcase(&cluster)= %then
%do;
proc sort data=&data;
by &strata descending &y;
run;
%end;
%if %upcase(&cluster) NE & %upcase(&strata)= %then
%do;
proc sort data=&data;
by &cluster descending &y;
run;
%end;
%if %upcase(&cluster) NE & %upcase(&strata)= NE %then
%do;
proc sort data=&data;
by &strata &cluster descending
&y;
run;
%end;
/*Default link functions*/
%if %upcase(&link)= & %upcase(&dist)=POISSON %then
%do;
%let link=log;
%end;
%else;
%if %upcase(&link)=%upcase(&link) &
%upcase(&dist)=POISSON %then
%do;
%let link=&link;
%end;
%if %upcase(&link)= & %upcase(&dist)=BINOMIAL %then
%do;
%let link=logit;
%end;
%else;
%if %upcase(&link)=%upcase(&link) &
%upcase(&dist)=BINOMIAL %then
%do;
%let link=&link;
%end;
%if %upcase(&link)= & %upcase(&dist)=GAMMA %then
%do;
%let link=inv;
%end;
%else;
%if %upcase(&link)=%upcase(&link) &
%upcase(&dist)=GAMMA %then
%do;
%let link=&link;
%end;
%if %upcase(&link)= & %upcase(&dist)=IG %then
%do;
%let link=inv2;
%end;
%else;
%if %upcase(&link)=%upcase(&link) &
%upcase(&dist)=IG %then
%do;
%let link=&link;
%end;
%if %upcase(&link)= & %upcase(&dist)=NORMAL %then
%do;
%let link=identity;
%end;
%else;
%if %upcase(&link)=%upcase(&link) &
%upcase(&dist)=NORMAL %then
%do;
%let link=&link;
%end;
proc iml;
MaxIter=&maxiter;
eps=&eps;
*Input data;
use &data;
read all var{&y} into
y[colname=depname];
read all var{&x} into
x0[colname=varname];
%if %upcase(&xzip) NE %then
%do;
read all var{&xzip} into
G0[colname=varnamezip];
G0=j(nrow(y),1,1)||G0;
%end;
%else
%do;
G0=j(nrow(y),1,1);
%end;
n=nrow(y);
nt=nrow(y);
%if %upcase(&weight) NE %then
%do;
read all var{&weight} into
w[colname=weightname];
%end;
%else
%do;
w=repeat(1, n, 1);
%end;
%if %upcase(&strata) NE %then
%do;
read all var{&strata} into h1;
%end;
%if %upcase(&fpc) NE %then
%do;
read all var{&fpc} into fh;
%end;
%else
%do;
fh=repeat(0, n, 1);
%end;
%if %upcase(&cluster) NE %then
%do;
read all var{&cluster} into i;
%end;
%if %upcase(&domain) NE %then
%do;
read all var{&domain} into domain;
domainclass=unique(domain);
vvd=repeat(w,1,ncol(domainclass));
do di=1 to ncol(domainclass);
vvd[loc(domain^=domainclass[di]),di]=0;
lld=loc(domain=domainclass[di]);
nd=ncol(lld);
end;
%end;
%if %upcase(&offset) NE %then
%do;
read all var{&offset} into offset;
%end;
%else
%do;
offset=repeat(0, n, 1);
%end;
close &data;
X=x0;
%if %upcase(&intercept)=Y %then
%do;
X=repeat(1,nrow(y),1)||x0;
%end;
p=ncol(X)*ncol(y);
weight=w;
%if %upcase(&domain) NE %then
%do;
w2=w;
y2=y;
x2=x;
G2=G0;
offset2=offset;
w=w[lld];
weight=w;
y=y[lld];
x=x[lld,];
G0=G0[lld,];
offset=offset[lld,];
n=nrow(y);
%end;
tab=y||x||w;
%if %upcase(&dist)=BINOMIAL %then
%do;
tab=y||x||w;
tabb=y||x||w;
call sort(tabb,{1});
freq=j(1,2,1);
do j=2 to nrow(tabb);
if tabb[j,1]=tabb[j-1,1] then
do;
freq[nrow(freq),1]=tabb[j-1,1];
freq[nrow(freq),2]=freq[nrow(freq),2]+1;
end;
else freq=freq//j(1,2,1);
end;
freq=(1:nrow(freq))`||freq;
mi=1;
y1=y;
yd=j(nrow(y),nrow(freq)-1,0);
do f=2 to nrow(freq);
do ff=1 to nrow(y);
if y[ff]=freq[f,2] then
yd[ff,f-1]=1;
else yd[ff,f-1]=0;
end;
end;
y=yd;
%end;
ni=ncol(unique(i));
nh=ncol(unique(h1));
sw=w[+];
smy=sum(y#w);
my=smy/sw;
title3 "The SURVEYGENMOD Macro";
%if %upcase(&dist) NE BINOMIAL %then
%do;
report=n//sw//my//smy;
reportname={"Number of Observations" "Sum of Weights" "Weighted Mean of &y." "Weighted Sum of &y."};
%if %upcase(&domain) NE %then
%do;
report=nt//nd//sw//my//smy;
reportname={"Number of Observations" "Number of Observations in Domain" "Sum of Weights" "Weighted Mean of &y." "Weighted Sum of &y."};
%end;
print report[label="Data Summary"
rowname=reportname];
%end;
*GLM routine;
start GLM;
* Get initial values;
%if %upcase(&link)=INV %then
%do;
yc = y + .000000000001;
y_trans = 1/yc;
%end;
%if %upcase(&link)=LOG %then
%do;
yc = (y + y[:])/2;
y_trans=log(yc);
%end;
%if %upcase(&link)=INV2 %then
%do;
yc = y + .000000000001;
y_trans = sqrt(1/yc);
%end;
%if %upcase(&link)=IDENTITY %then
%do;
yc = y;
y_trans=yc;
%end;
%if %upcase(&link)=LOGIT %then
%do;
yc=j(nrow(y), ncol(y),0);
do f=1 to ncol(y);
yc[,f]=y[+,f]/n;
end;
y_trans = log(yc/(1-yc[,+]));
%end;
%if %upcase(&link)=INV2 %then
%do;
b=repeat(.000000000001,p,1);
%end;
%else
%do;
b=j(ncol(x), ncol(y),0);
do ii=1 to ncol(y);
b[,ii]=inv(X`*(X#weight))*X`*(y_trans[,ii]#weight);
end;
%end;
phi=repeat(1, n, 1);
eta=j(nrow(y), ncol(y),0);
mu=j(nrow(y),ncol(y),0);
si=j(ncol(y)*ncol(x), nrow(y),0);
we=j(nrow(y), ncol(y),0);
wo=j(nrow(y), ncol(y),0);
H=j(ncol(y)*ncol(x), ncol(y)*ncol(x),0);
*IRLS algorithm;
k=1;
do iter=1 to MaxIter until(sum(abs((boldb)/(
oldb+0.0000001)))<eps);
*oldb=shapecol(b, ncol(x),ncol(y));
oldb=T(shape(b`, ncol(y),ncol(x)));
do ii=1 to ncol(y);
eta[,ii]=X*oldb[,ii];
end;
eta=eta;
/*Link functions*/
%if %upcase(&link)=INV %then
%do;
eta=choose(eta=0,0.0000000000001,eta);
mu=1/eta;
/*Firt derivative - link function*/
gderiv=-1/mu##2;
/*Second derivative - link
function*/
gderiv2=2/(mu##3);
%end;
%if %upcase(&link)=LOG %then
%do;
mu=exp(eta+offset);
/*Firt derivative - link function*/
gderiv=1/mu;
/*Second derivative - link
function*/
gderiv2=-1/(mu##2);
%end;
%if %upcase(&link)=IDENTITY %then
%do;
mu=eta;
/*Firt derivative - link function*/
gderiv=1;
/*Second derivative - link
function*/
gderiv2=0;
%end;
%if %upcase(&link)=INV2 %then
%do;
eta=choose(eta<0,0.0000000000001,eta);
mu=sqrt(1/eta);
/*Firt derivative - link function*/
gderiv=-2/(mu##3);
/*Second derivative - link
function*/
gderiv2=6/(mu##4);
%end;
%if %upcase(&link)=LOGIT %then
%do;
do ii=1 to ncol(y);
eta[,ii]=exp(X*oldb[,ii]);
end;
mu=mi#(eta/(1+eta[,+]));
/*Firt derivative - link function*/
gderiv=mi/(mu#(mi-mu));
/*Second derivative - link
function*/
gderiv2=(mi#(2#mu-mi))/((mu#(mimu))##2);
%end;
%if %upcase(&dist)=GAMMA %then
%do;
/*Variance function*/
vmu=mu##2;
/*Firt derivative - variance
function*/
vmuderiv=2#mu;
%end;
%if %upcase(&dist)=NORMAL %then
%do;
/*Variance function*/
vmu=1;
/*Firt derivative - variance
function*/
vmuderiv=0;
%end;
%if %upcase(&dist)=POISSON %then
%do;
/*Variance function*/
vmu=mu;
/*Firt derivative - variance
function*/
vmuderiv=1;
%end;
%if %upcase(&dist)=IG %then
%do;
/*Variance function*/
vmu=mu##3;
/*Firt derivative - variance
function*/
vmuderiv=3#(mu##2);
%end;
%if %upcase(&dist)=BINOMIAL %then
%do;
/*Variance function*/
vmu=(1/mi)#mu#(1-mu);
/*Firt derivative - variance
function*/
vmuderiv=(1/mi)#(1-2#mu);
%end;
%if %upcase(&dist)=NEGBIN %then
%do;
/*Variance function*/
vmu=mu+k#(mu##2);
/*First derivative - variance
function*/
vmuderiv=1+2#k#mu;
aux1=0;
aux2=0;
par=1;
dpar=1;
parold=par;
do while (abs(dpar)>eps);
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+mu)-(par+y)/(par+mu));
hess=sum(trigamma(par+y)-
trigamma(par)+1/par-
2/(par+mu)+(y+par)/((par+mu)#(par+mu)));
hess=choose(hess=0,1E-23,hess);
par0=par;
par=par0-inv(hess)*g;
dpar=par-par0;
end;
k=1/par;
ser=sqrt(-1/hess);
sek=ser/(par**2);
*print k par sek aux1 aux2;
%end;
/*Gradient vector*/
do ii=1 to ncol(y);
m1=(ii-1)*ncol(x)+1;
m2=ii*ncol(x);
si[m1:m2,] =
((weight#(y[,ii]-
mu[,ii]))/(vmu[,ii]#gderiv[,ii]#phi))`#X`;
end;
s = si[,+];
/*Weights*/
do ii=1 to ncol(y);
we[,ii]=weight/(phi#vmu[,ii]#(gderiv[,ii]##2))
;
wo[,ii]=we[,ii]+weight#(y[,ii]-
mu[,ii])#((vmu[,ii]#gderiv2[,ii]+vmuderiv[,ii]#gderiv[,ii])/((vmu[,ii]##2)#(gderiv[,ii]##3)#phi));
end;
/*Hessian matrix*/
do ii=1 to ncol(y);
m1=(ii-1)*ncol(x)+1;
m2=ii*ncol(x);
H[m1:m2, m1:m2]=-(X#wo[,ii])`*X;
end;
*oldb=shapecol(oldb,0,1);
oldb=T(shape(oldb`,1,0));
b=oldb-inv(H)*s;
end;
finish;
%if %upcase(&dist) ne ZIP & %upcase(&dist) ne ZINB %then
%do;
run GLM;
%end;
%else
%do;
%if %upcase(&dist)=ZIP %then
%do;
pos0=loc(y=0);
pos1=loc(y>0);
yc = (y[pos1] + y[pos1][:])/2;
y_trans=log(yc);
beta=inv((X[pos1,]#weight[pos1])`*X[pos1,])*(X[pos1,]#weight[pos1])`*y_trans;
*print beta;
*print pos0 pos1;
mi=1;
lambda0=((ncol(pos0)-sum(exp(-exp(X*beta))))/nrow(y));
lambda0=log(lambda0/(1-lambda0));
%if &xzip ne %then
%do;
lambda=lambda0//j(ncol(G0)-1,1,0);
%end;
%else
%do;
lambda=lambda0;
%end;
*print lambda;
zk=1/(1+exp(G0*lambda+exp(X*beta+offset)));
zk1=zk;
zk=choose(y>0,0,zk);
*print y zk;
dllike=1;
llike=0;
j=1;
do while (abs(dllike)>eps);
oldllike=llike;
dev=0;
ddev=1;
eta=X*beta;
mu=exp(eta+offset);
do while (abs(ddev)>eps);
gderiv=1/mu;
gderiv2=-1/(mu##2);
vmu=mu;
vmuderiv=1;
z=eta+(y-mu)#gderiv;
w=((1-zk)#weight)/(gderiv#gderiv#vmu);
beta=inv((X#w)`*X)*(X#w)`*z;
eta=X*beta;
mu=exp(eta+offset);
olddev=dev;
dev=sum((1-zk)#(y#eta-mu));
ddev=dev-olddev;
*print beta dev olddev ddev;
end;
H=-(X#w)`*X;
dev=0;
ddev=1;
eta=G0*lambda;
pi=exp(eta)/(1+exp(eta));
do while (abs(ddev)>eps);
gderiv=mi/(pi#(mi-pi));
gderiv2=(mi#(2#pi-mi))/((pi#(mi-pi))##2);
vmu=(1/mi)#pi#(1-pi);
vmuderiv=(1/mi)#(1-2#pi);
z=eta+(zk-pi)#gderiv;
w=weight/(gderiv#gderiv#vmu);
lambda=inv((G0#w)`*G0)*(G0#w)`*z;
eta=G0*lambda;
pi=exp(eta)/(1+exp(eta));
olddev=dev;
dev=sum(zk#eta)-sum(log(1+exp(eta)));
ddev=dev-olddev;
*print lambda dev olddev ddev;
end;
zk=1/(1+exp(-G0*lambda-exp(X*beta+offset)));
zk=choose(y>0,0,zk);
llike=sum(log(exp(G0[pos0,]*lambda)+exp(-
exp(X[pos0,]*beta))))+sum(y[pos1]#(X[pos1,]*beta)-
exp(X[pos1,]*beta))
-sum(log(1+exp(G0*lambda)))-sum(log(fact(y[pos1])));
dllike=llike-oldllike;
*print j beta lambda llike dllike;
j=j+1;
end;
b=beta;
phi=repeat(1, n, 1);
gderiv=1/mu;
gderiv2=-1/(mu##2);
vmu=mu+mu##2#zk/(1-zk);
vmuderiv=1+2#mu#zk/(1-zk);
we=weight/(phi#vmu#(gderiv##2));
wo=we+weight#(ymu)#((
vmu#gderiv2+vmuderiv#gderiv)/((vmu##2)#(gderiv##3)
#phi));
%end;
%if %upcase(&dist)=ZINB %then
%do;
pos0=loc(y=0);
pos1=loc(y>0);
yc = (y[pos1] + y[pos1][:])/2;
y_trans=log(yc);
beta=inv(X[pos1,]`*X[pos1,])*X[pos1,]`*y_trans;
*print
beta;
*print pos0 pos1;
mi=1;
lambda0=((ncol(pos0)-sum(exp(-exp(X*beta))))/nrow(y));
lambda0=log(lambda0/(1-lambda0));
%if &xzip ne %then
%do;
lambda=lambda0//j(ncol(G0)-1,1,0);
%end;
%else
%do;
lambda=lambda0;
%end;
*print lambda;
zk=1/(1+exp(G0*lambda+exp(X*beta+offset)));
zk1=zk;
zk=choose(y>0,0,zk);
*print y zk;
par=1;
dllike=1;
llike=0;
j=1;
do while (abs(dllike)>eps);
oldllike=llike;
olddllike=dllike;
dev=0;
ddev=1;
eta=X*beta;
mu=exp(eta+offset);
aux1=1;
do while (abs(ddev)>eps);
gderiv=1/mu;
gderiv2=-1/(mu##2);
dpar=1;
parold=par;
aux2=1;
do while (abs(dpar)>eps);
if par<0 then
par=0.00001;
gf=sum(digamma(par+(1-zk)#y)-
digamma(par)+log(par)+1-log(par+(1-zk)#mu)-(par+(1-
zk)#y)/(par+(1-zk)#mu));
hess=sum(trigamma(par+(1-zk)#y)-
trigamma(par)+1/par-2/(par+(1-zk)#mu)+((1-
zk)#y+par)/((par+(1-zk)#mu)#(par+(1-zk)#mu)));
hess=choose(hess=0,1E-23,hess);
par0=par;
par=par0-inv(hess)*gf;
dpar=par-par0;
if aux2>100 then
do;
dpar=0;
end;
aux2=aux2+1;
end;
k=1/par;
*print k par;
vmu=mu+k#(mu##2);
vmuderiv=1+2#k#mu;
z=eta+(y-mu)#gderiv;
w=((1-zk)#weight)/(gderiv#gderiv#vmu);
wo=w+weight#(y-mu)#(vmu#gderiv2 +
vmuderiv#gderiv)/(vmu#vmu#gderiv#gderiv#gderiv);
beta=inv((X#w)`*X)*(X#w)`*z;
eta=X*beta;
mu=exp(eta+offset);
olddev=dev;
gamma1=j(nrow(y),1,1);
do kkk=1 to nrow(y);
if y[kkk]>0 then
do;
do kkkk=1 to y[kkk];
if gamma1[kkk]<1E290 then
do;
if par>1E17 then
gamma1[kkk]=gamma1[kkk]*(y[kkk]+parkkkk)/( 1E300);
else gamma1[kkk]=gamma1[kkk]*(y[kkk]+parkkkk)/( par+mu[kkk])##y[kkk];
end;
end;
end;
end;
gamma1=choose(gamma1<1E-320,1E-320,gamma1);
log1=(par/(par+mu))##par;
dev=sum((1-zk)#(log(gamma1)-
log(gamma(y+1))+log(log1)+y#log(mu)));
ddev=dev-olddev;
*print beta dev olddev ddev;
if aux1>100 then
ddev=0;
aux1=aux1+1;
end;
H=-(X#wo)`*X;
dev=0;
ddev=1;
eta=G0*lambda;
pi=exp(eta)/(1+exp(eta));
do while (abs(ddev)>eps);
gderiv=mi/(pi#(mi-pi));
gderiv2=(mi#(2#pi-mi))/((pi#(mi-pi))##2);
vmu=(1/mi)#pi#(1-pi);
vmuderiv=(1/mi)#(1-2#pi);
z=eta+(zk-pi)#gderiv;
w=weight/(gderiv#gderiv#vmu);
lambda=inv((G0#w)`*G0)*(G0#w)`*z;
eta=G0*lambda;
pi=exp(eta)/(1+exp(eta));
olddev=dev;
dev=sum(zk#eta)-sum(log(1+exp(eta)));
ddev=dev-olddev;
*print lambda dev olddev ddev;
end;
zk=1/(1+exp(-
G0*lambda)#(par/(par+exp(X*beta+offset)))##par);
zk=choose(y>0,0,zk);
*llike=sum(log((exp(X[pos0,]*beta)+(par/(par+exp(X[pos0,
]*beta)))##par)/(1+exp(G0[pos0,]*lambda))))+
sum(log(gamma1[pos1])+par#log(par/(par+exp(X[pos1,]*beta
)))+
y[pos1]#(X[pos1,]*beta)-log(1+exp(G0[pos1,]*lambda)));
llike=sum(zk#(G0*lambda)-log(1+exp(G0*lambda)))+
sum((1-zk)#(log(gamma1)+par#log(par/(par+exp(X*beta)))+
y#(exp(X*beta)/(par+exp(X*beta)))));
dllike=llike-oldllike;
*print j beta k lambda llike dllike;
if j>150 then
dllike=0;
j=j+1;
end;
/**** fitting a zip model when k<1E10 ************/
if par>1E10 then
do;
k=1E-8;
par=1/k;
beta=inv(X[pos1,]`*X[pos1,])*X[pos1,]`*y_trans;
*print
beta;
lambda0=((ncol(pos0)-sum(exp(-exp(X*beta))))/nrow(y));
lambda0=log(lambda0/(1-lambda0));
%if &xzip ne %then
%do;
lambda=lambda0//j(ncol(G0)-1,1,0);
%end;
%else
%do;
lambda=lambda0;
%end;
*print lambda;
zk=1/(1+exp(G0*lambda+exp(X*beta+offset)));
zk1=zk;
zk=choose(y>0,0,zk);
*print y zk;
dllike=1;
llike=0;
j=1;
do while (abs(dllike)>eps);
oldllike=llike;
olddllike=dllike;
dev=0;
ddev=1;
eta=X*beta;
mu=exp(eta+offset);
aux1=1;
do while (abs(ddev)>eps);
gderiv=1/mu;
gderiv2=-1/(mu##2);
vmu=mu+k#(mu##2);
vmuderiv=1+2#k#mu;
z=eta+(y-mu)#gderiv;
w=((1-zk)#weight)/(gderiv#gderiv#vmu);
wo=w+weight#(y-mu)#(vmu#gderiv2 +
vmuderiv#gderiv)/(vmu#vmu#gderiv#gderiv#gderiv);
beta=inv((X#w)`*X)*(X#w)`*z;
eta=X*beta;
mu=exp(eta+offset);
olddev=dev;
gamma1=j(nrow(y),1,1);
do kkk=1 to nrow(y);
if y[kkk]>0 then
do;
do kkkk=1 to y[kkk];
if gamma1[kkk]<1E290 then
do;
if par>1E17 then
gamma1[kkk]=gamma1[kkk]*(y[kkk]+parkkkk)/( 1E300);
else gamma1[kkk]=gamma1[kkk]*(y[kkk]+parkkkk)/( par+mu[kkk])##y[kkk];
end;
end;
end;
end;
gamma1=choose(gamma1<1E-320,1E-320,gamma1);
log1=(par/(par+mu))##par;
dev=sum((1-zk)#(log(gamma1)-
log(gamma(y+1))+log(log1)+y#log(mu)));
ddev=dev-olddev;
*print beta dev olddev ddev;
if aux1>100 then
ddev=0;
aux1=aux1+1;
end;
H=-(X#wo)`*X;
dev=0;
ddev=1;
eta=G0*lambda;
pi=exp(eta)/(1+exp(eta));
do while (abs(ddev)>eps);
gderiv=mi/(pi#(mi-pi));
gderiv2=(mi#(2#pi-mi))/((pi#(mi-pi))##2);
vmu=(1/mi)#pi#(1-pi);
vmuderiv=(1/mi)#(1-2#pi);
z=eta+(zk-pi)#gderiv;
w=weight/(gderiv#gderiv#vmu);
lambda=inv((G0#w)`*G0)*(G0#w)`*z;
eta=G0*lambda;
pi=exp(eta)/(1+exp(eta));
olddev=dev;
dev=sum(zk#eta)-sum(log(1+exp(eta)));
ddev=dev-olddev;
*print lambda dev olddev ddev;
end;
zk=1/(1+exp(-G0*lambda)#(par/(par+exp(X*beta+offset)))##par);
zk=choose(y>0,0,zk);
*llike=sum(log((exp(G0[pos0,]*lambda)+(par/(par+exp(X[po
s0,]*beta+offset[pos0])))##par)/(1+exp(G0[pos0,]*lambda)
)))+
sum(log(gamma1[pos1])+par#log(par/(par+exp(X[pos1,]*beta
+offset[pos1])))+
y[pos1]#(X[pos1,]*beta+offset[pos1])-
log(1+exp(G0[pos1,]*lambda)));
llike=sum(zk#(G0*lambda)-log(1+exp(G0*lambda)))+
sum((1-zk)#(log(gamma1)+par#log(par/(par+exp(X*beta)))+
y#(exp(X*beta)/(par+exp(X*beta)))));
dllike=llike-oldllike;
*print j beta k lambda llike dllike;
if j>150 then
dllike=0;
j=j+1;
end;
end;
/********************* end **********************/
hess=sum(trigamma(par+(1-zk)#y)-trigamma(par)+1/par-
2/(par+(1-zk)#mu)+((1-zk)#y+par)/((par+(1-
zk)#mu)#(par+(1-zk)#mu)));
ser=sqrt(abs(1/hess));
sek=ser/(par**2);
b=beta;
phi=repeat(1, n, 1);
gderiv=1/mu;
gderiv2=-1/(mu##2);
vmu=mu+mu##2#(zk/(1-zk)+k/(1-zk));
vmuderiv=1+2#mu#(zk/(1-zk)+k/(1-zk));
we=weight/(phi#vmu#(gderiv##2));
wo=we+weight#(ymu)#((
vmu#gderiv2+vmuderiv#gderiv)/((vmu##2)#(gderiv##3)
#phi));
%end;
%end;
*b1=shapecol(b, ncol(x),ncol(y));
b1=T(shape(b`, ncol(y),ncol(x)));
%if %upcase(&dist)=BINOMIAL %then
%do;
*b = shapecol(shapecol(b,ncol(x),ncol(y))`,ncol(x)*ncol(y),1);
b = T(shape(T(shape(b`,ncol(y),
ncol(x)))`,1,ncol(x)*ncol(y)));
%end;
/*Predict and residuals*/
do ii=1 to ncol(y);
eta[,ii]=X*b1[,ii];
end;
%if %upcase(&domain) NE %then
%do;
y=y2;
x=x2;
G0=G2;
offset=offset2;
w=vvd[,ncol(domainclass)];
weight=w;
n=nrow(y);
phi=repeat(1, n, 1);
tab=y||x||w;
eta=j(nrow(y), ncol(y),0);
mu=j(nrow(y),ncol(y),0);
do ii=1 to ncol(y);
eta[,ii]=X*b[,ii];
end;
/*Link functions*/
%if %upcase(&link)=INV %then
%do;
yc = y + .000000000001;
y_trans = 1/yc;
eta=choose(eta=0,0.0000000000001,eta);
mu=1/eta;
/*Firt derivative - link function*/
gderiv=-1/mu##2;
/*Second derivative - link function*/
gderiv2=2/(mu##3);
%end;
%if %upcase(&link)=LOG %then
%do;
yc = (y + y[:])/2;
y_trans=log(yc);
mu=exp(eta);
/*Firt derivative - link function*/
gderiv=1/mu;
/*Second derivative - link function*/
gderiv2=-1/(mu##2);
%end;
%if %upcase(&link)=INV2 %then
%do;
yc = y + .000000000001;
y_trans = sqrt(1/yc);
eta=choose(eta<0,0.0000000000001,eta);
mu=sqrt(1/eta);
/*Firt derivative - link function*/
gderiv=-2/(mu##3);
/*Second derivative - link function*/
gderiv2=6/(mu##4);
%end;
%if %upcase(&link)=IDENTITY %then
%do;
yc = y;
y_trans=yc;
mu=eta;
/*Firt derivative - link function*/
gderiv=1;
/*Second derivative - link function*/
gderiv2=0;
%end;
%if %upcase(&link)=LOGIT %then
%do;
yc=j(nrow(y), ncol(y),0);
do f=1 to ncol(y);
yc[,f]=y[+,f]/n;
end;
y_trans = log(yc/(1-yc[,+]));
do ii=1 to ncol(y);
eta[,ii]=exp(X*oldb[,ii]);
end;
mu=mi#(eta/(1+eta[,+]));
/*Firt derivative - link function*/
gderiv=mi/(mu#(mi-mu));
/*Second derivative - link function*/
gderiv2=(mi#(2#mu-mi))/((mu#(mi-mu))##2);
%end;
%if %upcase(&dist)=GAMMA %then
%do;
/*Variance function*/
vmu=mu##2;
/*Firt derivative - variance function*/
vmuderiv=2#mu;
%end;
%if %upcase(&dist)=NORMAL %then
%do;
/*Variance function*/
vmu=1;
/*Firt derivative - variance function*/
vmuderiv=0;
%end;
%if %upcase(&dist)=POISSON %then
%do;
/*Variance function*/
vmu=mu;
/*Firt derivative - variance function*/
vmuderiv=1;
%end;
%if %upcase(&dist)=ZIP %then
%do;
pi=exp(G0*lambda)/(1+exp(G0*lambda));
zk=1/(1+exp(-G0*lambda-exp(X*beta+offset)));
/*Variance function*/
vmu=mu+mu##2#zk/(1-zk);
/*Firt derivative - variance function*/
vmuderiv=1+2#mu#zk/(1-zk);
%end;
%if %upcase(&dist)=ZINB %then
%do;
pi=exp(G0*lambda)/(1+exp(G0*lambda));
zk=1/(1+exp(-
G0*lambda)#(par/(par+exp(X*beta+offset)))##par);
/*Variance function*/
vmu=mu+(mu##2)#(zk/(1-zk)+k/(1-zk));
/*Firt derivative - variance function*/
vmuderiv=1+2#mu#(zk/(1-zk)+k/(1-zk));
%end;
%if %upcase(&dist)=IG %then
%do;
/*Variance function*/
vmu=mu##3;
/*Firt derivative - variance function*/
vmuderiv=3#(mu##2);
%end;
%if %upcase(&dist)=BINOMIAL %then
%do;
/*Variance function*/
vmu=(1/mi)#mu#(1-mu);
/*Firt derivative - variance function*/
vmuderiv=(1/mi)#(1-2#mu);
%end;
%if %upcase(&dist)=NEGBIN %then
%do;
/*Variance function*/
vmu=mu+k#(mu##2);
/*First derivative - variance function*/
vmuderiv=1+2#k#mu;
%end;
we=weight/(phi#vmu#(gderiv##2));
wo=we+weight#(ymu)#((
vmu#gderiv2+vmuderiv#gderiv)/((vmu##2)#(gderiv##3)
#phi));
%end;
e=y-mu;
e2=j(nrow(x), ncol(y)*ncol(x),0);
do ii=1 to ncol(y);
e2[,(ii-1)*ncol(x)+1:ii*ncol(x)]=(phi#we#e[,ii]#gderiv)#X;
end;
G=j(ncol(x)*ncol(y), ncol(x)*ncol(y), 0);
/*G matrix*/
*Strata-Cluster;
if nh>0 & ni>0 then
do;
_h1_=j(nrow(h1),1,0);
_h1_[1]=1;
_i_=j(nrow(i),1,0);
_i_[1]=1;
do kk=2 to nrow(h1);
if h1[kk]=h1[kk-1] then
_h1_[kk]=_h1_[kk-1];
else _h1_[kk]=_h1_[kk- 1]+1;
if i[kk]=i[kk-1] then
_i_[kk]=_i_[kk-1];
else _i_[kk]=_i_[kk-1]+1;
end;
tab=tab||_h1_||_i_;
do ii=1 to ncol(y);
m1=(ii-1)*ncol(x)+1;
m2=ii*ncol(x);
ehi=j(1,ncol(x)+4,1);
ehi[1,3:2+ncol(x)]=e2[1,(ii-
1)*ncol(x)+1:ii*ncol(x)];
ehi[1]=tab[1,ncol(tab)-1];
ehi[2]=tab[1,ncol(tab)];
ehi[ncol(ehi)-1]=fh[1];
do j=2 to nrow(tab);
if
tab[j,ncol(tab)-1]=tab[j-1,ncol(tab)-1] &
tab[j,ncol(tab)]=tab[j-1,ncol(tab)] then
do;
ehi[nrow(ehi),3:2+ncol(x)]=ehi[nrow(ehi),3:2+ncol(x)]+e2[j,(ii-1)*ncol(x)+1:ii*ncol(x)];
ehi[nrow(ehi),ncol(ehi)]=ehi[nrow(ehi),ncol(ehi)]+1;
ehi[nrow(ehi),ncol(ehi)-
1]=(ehi[nrow(ehi),ncol(ehi)-1]+fh[j])/2;
end;
else
ehi=ehi//(tab[j,ncol(tab)-
1]||tab[j,ncol(tab)]||e2[j,(ii-
1)*ncol(x)+1:ii*ncol(x)]||fh[j]||j(1,1,1));
end;
ee=ehi[,ncol(ehi)-
1:ncol(ehi)];
ehi=ehi[,1:ncol(ehi)-
2]||j(nrow(ehi),3+ncol(x),0);
ehi[1,3+ncol(x):2+2*ncol(x)]=ehi[1,3:2+ncol(x)
];
ehi[,ncol(ehi)-
2:ncol(ehi)-1]=ee;
ehi[1,ncol(ehi)]=ehi[1,ncol(ehi)-1];
count=0;
do j=2 to nrow(ehi);
if
ehi[j,1]=ehi[j-1,1] then
do;
ehi[j,3+ncol(x):2+2*ncol(x)]=ehi[j-
1,3+ncol(x):2+2*ncol(x)]+ehi[j,3:2+ncol(x)];
count=count+1;
ehi[j,ncol(ehi)]=ehi[j-
1,ncol(ehi)]+ehi[j,ncol(ehi)-1];
ehi[j,ncol(ehi)-2]=ehi[j-1,ncol(ehi)-2];
end;
else
do;
if
ehi[j-1,1]=1 then
do;
ehi[1:count+1,ncol(ehi)-1]=count+1;
ehi[1:count+1,3+ncol(x):2+2*ncol(x)]=repeat(ehi[j-1,3+ncol(x):2+2*ncol(x)]/ehi[j-1,ncol(ehi)-1],count+1);
in=ehi[j,2];
ehi[1:count+1,ncol(ehi)]=repeat(ehi[j-
1,ncol(ehi)],count+1);
end;
else
do;
ehi[in:in+count,ncol(ehi)-1]=count+1;
ehi[in:in+count,3+ncol(x):2+2*ncol(x)]=repeat(
ehi[j-1,3+ncol(x):2+2*ncol(x)]/ehi[j-1,ncol(ehi)-
1],count+1);
ehi[in:in+count,ncol(ehi)]=repeat(ehi[j-
1,ncol(ehi)],count+1);
in=ehi[j,2];
end;
ehi[j,3+ncol(x):2+2*ncol(x)]=ehi[j,3:2+ncol(x)];
ehi[j,ncol(ehi)]=ehi[j,ncol(ehi)-1];
count=0;
end;
if j=nrow(ehi) then
do;
ehi[in:in+count,ncol(ehi)-1]=count+1;
ehi[in:in+count,3+ncol(x):2+2*ncol(x)]=repeat(
ehi[j,3+ncol(x):2+2*ncol(x)]/ehi[j,ncol(ehi)-
1],count+1);
ehi[in:in+count,ncol(ehi)]=repeat(ehi[j,ncol(ehi)],count+1);
in=ehi[j,2];
end;
do jj=1 to nrow(ehi);
if
ehi[jj,ncol(ehi)-1]=1 then
do;
ehi[jj,ncol(ehi)-
1]=2;
ehi[jj,3:2+ncol(x)]=0;
ehi[jj,3+ncol(x):2+2*ncol(x)]
=0;
end;
end;
end;
G[m1:m2,
m1:m2]=((n-1)/(n-ncol(x)))*(((ehi[,ncol(ehi)-
1]/(ehi[,ncol(ehi)-1]-1))#(ehi[,3:2+ncol(x)]-
ehi[,3+ncol(x):2+2*ncol(x)]))`*(ehi[,3:2+ncol(x)]-
ehi[,3+ncol(x):2+2*ncol(x)]));
%if
%upcase(&fpc)^= %then
%do;
G[m1:m2,
m1:m2]=((n-1)/(n-ncol(x)))*((((ehi[,ncol(ehi)-1]#(1-
ehi[,ncol(ehi)-2]))/(ehi[,ncol(ehi)-1]-
1))#(ehi[,3:(2+ncol(x))]-
ehi[,3+ncol(x):2+2*ncol(x)]))`*(ehi[,3:(2+ncol(x))]-
ehi[,3+ncol(x):2+2*ncol(x)]));
%end;
*print G;
end;
%if %upcase(&dist)=BINOMIAL %then
%do;
Report="%upcase(&data)"//"%upcase(&dist)"//"%upcase(&link)"//"&y"//"&strata"//"&cluster"
//"&weight"//"Ridge Stabilized Newton-Raphson"//"Taylor Series";
ReportName={"Data Set" "Distribution" "Link Function" "Dependent Variable"
"Stratum Variable" "Cluster Variable" "Weight Variable"
"Optimization Technique" "Variance Estimation Method"};
print Report[label="Model
Information" rowname=reportname];
print (nrow(freq))[label='Number of Response Levels'];
report=(ehi[nrow(ehi),1])//(ehi[nrow(ehi),2]);
reportnames={"Number of Strata",
"Number of Clusters"};
print report[label="Design Summary"
rowname=reportnames];
%end;
%else
%do;
Report="%upcase(&data)"//"%upcase(&dist)"//"%upcase(&link)"//"&y"//"&strata"//
"&cluster"//"&weight"//"Ridge
Stabilized Newton-Raphson"//"Taylor Series";
ReportName={"Data Set"
"Distribution" "Link Function" "Dependent Variable"
"Stratum Variable" "Cluster Variable" "Weight Variable"
"Optimization Technique" "Variance Estimation Method"};
print Report[label="Model Information" rowname=reportname];
report=(ehi[nrow(ehi),1])//(ehi[nrow(ehi),2]);
reportnames={"Number of Strata",
"Number of Clusters"};
print report[label="Design Summary"
rowname=reportnames];
%end;
df=ehi[nrow(ehi),2]-ehi[nrow(ehi),1];
end;
*Strata;
else if nh>0 & ni=0 then
do;
_h1_=j(nrow(h1),1,0);
_h1_[1]=1;
do kk=2 to nrow(h1);
if h1[kk]=h1[kk-1] then
_h1_[kk]=_h1_[kk-1];
else _h1_[kk]=_h1_[kk- 1]+1;
end;
tab=tab||_h1_;
do ii=1 to ncol(y);
ehi=j(1,2*ncol(x)+2,0);
ehi[1,2:1+ncol(x)]=e2[1,(ii-
1)*ncol(x)+1:ii*ncol(x)];
ehi[1]=tab[1,ncol(tab)];
ehi[ncol(ehi)]=1;
do j=2 to nrow(tab);
if tab[j,ncol(tab)]=tab[j-
1,ncol(tab)] then
do;
ehi[nrow(ehi),2:1+ncol(x)]=ehi[nrow(ehi),2:1+ncol(x)]+e2[j,(ii-1)*ncol(x)+1:ii*ncol(x)];
ehi[nrow(ehi),ncol(ehi)]=ehi[nrow(ehi),ncol(ehi)]+1;
end;
else
do;
ehi=ehi//(tab[j,ncol(tab)]||e2[j,]||j(1,ncol(x)+1,1))
;
end;
end;
ehi[,2+ncol(x):1+2*ncol(x)]=ehi[,2:1+ncol(x)]/
ehi[,ncol(ehi)];
do jj=1 to nrow(ehi);
ehi2=ehi2//repeat(ehi[jj,2+ncol(x):2+2*ncol(x)
],ehi[jj,ncol(ehi)]);
end;
eh=e2[,(ii-
1)*ncol(x)+1:ii*ncol(x)]||ehi2[,1:ncol(ehi2)-
1]||fh||ehi2[,ncol(ehi2)];
G[m1:m2, m1:m2]=((n-1)/(nncol(
x)))*(((eh[,ncol(eh)]/(eh[,ncol(eh)]-
1))#(eh[,1:ncol(x)]-
eh[,1+ncol(x):2*ncol(x)]))`*(eh[,1:ncol(x)]-
eh[,1+ncol(x):2*ncol(x)]));
%if %upcase(&fpc) NE %then
%do;
G[m1:m2, m1:m2]=((n-1)/(nncol(
x)))*((((eh[,ncol(eh)]#(1-eh[,ncol(eh)-
1]))/(eh[,ncol(eh)]-1))#(eh[,1:ncol(x)]-
eh[,1+ncol(x):2*ncol(x)]))`*(eh[,1:ncol(x)]-
eh[,1+ncol(x):2*ncol(x)]));
%end;
*print G;
free ehi2;
end;
%if %upcase(&dist)=BINOMIAL %then
%do;
Report="%upcase(&data)"//"%upcase(&dist)"//"%upcase(&link)"//"&y"//"&strata"
//"&weight"//"Ridge Stabilized Newton-Raphson"//"Taylor Series";
ReportName={"Data Set"
"Distribution" "Link Function" "Dependent Variable"
"Stratum Variable" "Weight Variable"
"Optimization Technique" "Variance Estimation Method"};
print Report[label="Model Information" rowname=reportname];
print (nrow(freq))[label="Number of Response Levels"];
report=(ehi[<>,1]);
reportnames={"Number of Strata"};
print report[label="Design Summary"
rowname=reportnames];
%end;
%else
%do;
Report="%upcase(&data)"//"%upcase(&dist)"//"%upcase(&link)"//"&y"//"&strata"//
"&weight"//"Ridge Stabilized Newton-Raphson"//"Taylor Series";
ReportName={"Data Set"
"Distribution" "Link Function" "Dependent Variable"
"Stratum Variable" "Weight Variable"
"Optimization Technique" "Variance Estimation Method"};
print Report[label="Model
Information" rowname=reportname];
report=(ehi[<>,1]);
reportnames={"Number of Strata"};
print report[label="Design Summary"
rowname=reportnames];
%end;
df=n-ehi[<>,1];
end;
*Cluster;
else if ni>0 & nh=0 then
do;
_i_=j(nrow(i),1,0);
_i_[1]=1;
do kk=2 to nrow(i);
if i[kk]=i[kk-1] then
_i_[kk]=_i_[kk-1];
else _i_[kk]=_i_[kk-1]+1;
end;
tab=tab||_i_;
do ii=1 to ncol(y);
m1=(ii-1)*ncol(x)+1;
m2=ii*ncol(x);
ehi=j(1,ncol(x)+3,0);
ehi[1,2:1+ncol(x)]=e2[1,(ii-
1)*ncol(x)+1:ii*ncol(x)];
ehi[1]=tab[1,ncol(tab)];
ehi[ncol(ehi)-
1]=fh[1];
ehi[1,ncol(ehi)]=1;
do j=2 to nrow(tab);
if
tab[j,ncol(tab)]=tab[j-1,ncol(tab)] then
do;
ehi[nrow(ehi),2:1+ncol(x)]=ehi[nrow(ehi),2:1+ncol(x)]+e2
[j,(ii-1)*ncol(x)+1:ii*ncol(x)];
ehi[nrow(ehi),ncol(ehi)-
1]=(ehi[nrow(ehi),ncol(ehi)-
1]+fh[j])/2;
ehi[nrow(ehi),ncol(ehi)]=ehi[nrow(ehi),ncol(
ehi)]+1;
end;
else
do;
ehi=ehi//(tab[j,ncol(tab)]||e2[j,(ii-
1)*ncol(x)+1:ii*ncol(x)]||fh[j]||j(1,1,1));
end;
end;
G[m1:m2,
m1:m2]=((n-1)/(n-ncol(x)))*(ehi[<>,1]/(ehi[<>,1]-
1))*((ehi[,2:1+ncol(x)]-
ehi[+,2:1+ncol(x)]/ehi[+,ncol(ehi)])`*(ehi[,2:1+ncol(x)]
-ehi[+,2:1+ncol(x)]/ehi[+,ncol(ehi)]));
%if
%upcase(&fpc) NE %then
%do;
G[m1:m2,
m1:m2]=((n-1)/(n-ncol(x)))*(ehi[<>,1]/(ehi[<>,1]-
1))*(((1-ehi[,ncol(ehi)-1])#(ehi[,2:1+ncol(x)]-
ehi[+,2:1+ncol(x)]/ehi[+,ncol(ehi)]))`*(ehi[,2:1+ncol(x)
]-ehi[+,2:1+ncol(x)]/ehi[+,ncol(ehi)]));
%end;
*print G;
end;
%if
%upcase(&dist)=BINOMIAL %then
%do;
Report="%upcase(&data)"//"%upcase(&dist)"//"%upcase(&link)"//"&y"//"&cluster"
//"&weight"//"Ridge Stabilized Newton-Raphson"//"Taylor Series";
ReportName={"Data Set" "Distribution" "Link
Function" "Dependent Variable"
"Cluster Variable" "Weight Variable"
"Optimization Technique" "Variance Estimation Method"};
print
Report[label="Model Information" rowname=reportname];
report=(ehi[<>,1]);
reportnames={"Number of Clusters"};
print
report[label="Design Summary" rowname=reportnames];
%end;
%else
%do;
Report="%upcase(&data)"//"%upcase(&dist)"//"%upcase(&link)"//"&y"//
"&cluster"//"&weight"//"Ridge Stabilized Newton-Raphson"//"Taylor Series";
ReportName={"Data Set" "Distribution" "Link Function" "Dependent Variable"
"Cluster Variable" "Weight Variable"
"Optimization Technique" "Variance Estimation Method"};
print
Report[label="Model Information" rowname=reportname];
report=(ehi[<>,1]);
reportnames={"Number of Clusters"};
print
report[label="Design Summary" rowname=reportnames];
%end;
df=ehi[<>,1]-1;
end;
else
do;
do ii=1 to
ncol(y);
m1=(ii-
1)*ncol(x)+1;
m2=ii*ncol(x);
G[m1:m2,
m1:m2]=-H[(ii-1)*ncol(x)+1:ii*ncol(x), (ii-
1)*ncol(x)+1:ii*ncol(x)];
end;
%if
%upcase(&dist)=BINOMIAL %then
%do;
Report="%upcase(&data)"//"%upcase(&dist)"//"%upcase(&link)"//"&y"
//"&weight"//"Ridge Stabilized Newton-Raphson"//"Taylor Series";
ReportName={"Data Set" "Distribution" "Link
Function" "Dependent Variable"
"Weight Variable" "Optimization Technique" "Variance Estimation Method"};
print
Report[label="Model Information" rowname=reportname];
%end;
%else
%do;
Report="%upcase(&data)"//"%upcase(&dist)"//"%upcase(&link)"//"&y"//
"&weight"//"Ridge Stabilized Newton-Raphson"//"Taylor Series";
ReportName={"Data Set" "Distribution" "Link Function" "Dependent Variable"
"Weight Variable" "Optimization Technique" "Variance Estimation Method"};
print
Report[label="Model Information" rowname=reportname];
%end;
df=nrow(y)-
ncol(x);
end;
%if
%upcase(&dist)=BINOMIAL %then
%do;
print freq[colname={"Order" "&y" "Frequency"} label='Response Profile'],
"Logit have used &y ="(trim(left(char(freq[1,2]))))"as the reference category.";
%end;
%if %upcase(&scale)=DEVIANCE %then
%do;
/*Deviance and Scale deviance*/
%if %upcase(&dist)=NORMAL %then
%do;
deviance=sum(weight#(ymu)##2);
%end;
%if %upcase(&dist)=POISSON %then
%do;
mu=choose(mu<=0,0.0000000000001,mu);
deviance=2#sum(weight#(yc#log(yc/mu)-(ycmu)));
%end;
%if %upcase(&dist)=NEGBIN %then
%do;
mu=choose(mu<=0,0.0000000000001,mu);
deviance=2#sum(yc#log(k*mu/weight)-
(yc+weight/k)#log((yc+weight/k)/(mu+weight/k)));
%end;
%if %upcase(&dist)=IG %then
%do;
mu=choose(mu=0,0.0000000000001,mu);
deviance=sum(weight#((ycmu)##2)/(yc#mu##2));
%end;
%if %upcase(&dist)=BINOMIAL %then
%do;
deviance=2#sum(weight#mi#(yc#log(yc/mu))+(1-yc)#log((1-yc)/(1-mu)));
%end;
%if %upcase(&dist)=ZIP %then
%do;
dev0=(weight#log(zk+(1-zk)#exp(-mu)));
dev1=(weight#(log(1-zk)+y#log(mu)-mu-log(fact(y))));
dev01=choose(y>0,dev1,dev0);
deviance=-2#sum(dev01);
%end;
%if %upcase(&dist)=ZINB %then
%do;
dev0=(log(zk+(1-
zk)#(1+mu#k/weight)##(-1/k)));
dev1=((log(1-
zk)+y#log(mu#k/weight)-
(y+weight/k)#log(1+mu#k/weight)));
dev01=choose(y>0,dev1,dev0);
deviance=-2#sum(dev01);
%end;
scale1=deviance/df;
%if &strata NE | &cluster NE %then
%do;
scale1=n#deviance/((np)#sum(weight));
%end;
scaled=sqrt(scale1);
%if %upcase(&dist)=GAMMA %then
%do;
mu=choose(mu<=0,0.0000000000001,mu);
deviance=2#sum(weight#(-log(yc/mu)+(yc-mu)/mu));
scale1=deviance/df;
%if &strata NE | &cluster NE %then
%do;
scale1=n#deviance/((np)#sum(weight));
%end;
scaled=1/scale1;
%end;
sdeviance=deviance/scale1;
/*Pearson and Scale pearson*/
pearson=j(ncol(y),1,0);
do ii=1 to ncol(y);
pearson[ii,]=sum((weight#(y[,ii]-
mu[,ii])##2)/vmu[,ii]);
end;
pearson=sum(pearson);
scale2=pearson/df;
%if &strata NE | &cluster NE %then
%do;
scale2=n#pearson/((np)#sum(weight));
%end;
spearson=pearson/scale1;
valeudf=sdeviance/df;
valeudf1=spearson/df;
dev=df||deviance||scale1;
sdev=df||sdeviance||valeudf;
pcs=df||pearson||scale2;
sp=df||spearson||valeudf1;
phi=scale1;
%end;
%if %upcase(&scale)=PEARSON %then
%do;
/*Pearson*/
pearson=j(ncol(y),1,0);
do ii=1 to ncol(y);
pearson[ii,]=sum((weight#(y[,ii]-
mu[,ii])##2)/vmu[,ii]);
end;
pearson=sum(pearson);
scale1=pearson/df;
%if &strata NE | &cluster NE %then
%do;
scale1=n#pearson/((np)#sum(weight));
%end;
/*Deviance and Scale pearson*/
%if %upcase(&dist)=NORMAL %then
%do;
deviance=sum(weight#(ymu)##2);
%end;
%if %upcase(&dist)=POISSON %then
%do;
mu=choose(mu<=0,0.0000000000001,mu);
deviance=2#sum(weight#(yc#log(yc/mu)-(ycmu)));
%end;
%if %upcase(&dist)=NEGBIN %then
%do;
mu=choose(mu<=0,0.0000000000001,mu);
deviance=2#sum(yc#log(k*mu/weight)-
(yc+weight/k)#log((yc+weight/k)/(mu+weight/k)));
%end;
%if %upcase(&dist)=IG %then
%do;
mu=choose(mu=0,0.0000000000001,mu);
deviance=sum(weight#((ycmu)##
2)/(yc#mu##2));
%end;
%if %upcase(&dist)=BINOMIAL %then
%do;
deviance=2*sum(weight#mi#(yc#log(yc/mu))+(1-
y)#log((1-yc)/(1-mu)));
%end;
%if %upcase(&dist)=ZIP %then
%do;
dev0=(weight#log(zk+(1-
zk)#exp(-mu)));
dev1=(weight#(log(1-
zk)+y#log(mu)-mu-log(fact(y))));
dev01=choose(y>0,dev1,dev0);
deviance=-2#sum(dev01);
%end;
%if %upcase(&dist)=ZINB %then
%do;
dev0=(log(zk+(1-
zk)#(1+mu#k/weight)##(-1/k)));
dev1=((log(1-
zk)+y#log(mu#k/weight)-
(y+weight/k)#log(1+mu#k/weight)));
dev01=choose(y>0,dev1,dev0);
deviance=-2#sum(dev01);
%end;
scale2=deviance/df;
%if &strata NE | &cluster NE %then
%do;
scale2=n*deviance/(np)#
sum(weight);
%end;
scaled=sqrt(scale1);
%if %upcase(&dist)=GAMMA %then
%do;
mu=choose(mu<=0,0.0000000000001,mu);
deviance=2#sum(weight#(-
log(yc/mu)+(yc-mu)/mu));
scale2=deviance/df;
%if &strata NE | &cluster NE %then
%do;
scale2=n*deviance/(n-p)#sum(weight);
%end;
scaled=1/scale1;
%end;
spearson=pearson/scale1;
sdeviance=deviance/scale1;
valeudf=sdeviance/df;
valeudf1=spearson/df;
dev=df||deviance||scale2;
sdev=df||sdeviance||valeudf;
pcs=df||pearson||scale1;
sp=df||spearson||valeudf1;
phi=scale1;
%end;
%else;
%if %upcase(&scale)= %then
%do;
/*Deviance and Scale deviance*/
%if %upcase(&dist)=POISSON %then
%do;
mu=choose(mu<=0,0.0000000000001,mu);
deviance=2#sum(weight#(yc#log(yc/mu)-(ycmu)));
%end;
%if %upcase(&dist)=NEGBIN %then
%do;
mu=choose(mu<=0,0.0000000000001,mu);
deviance=2#sum(yc#log(k*mu/weight)-
(yc+weight/k)#log((yc+weight/k)/(mu+weight/k)));
%end;
%if %upcase(&dist)=BINOMIAL %then
%do;
deviance=2#sum(weight#mi#(yc#log(yc/mu))+(1-
yc)#log((1-yc)/(1-mu)));
%end;
%if %upcase(&dist)=ZIP %then
%do;
dev0=(weight#log(zk+(1-
zk)#exp(-mu)));
dev1=(weight#(log(1-
zk)+y#log(mu)-mu-log(fact(y))));
dev01=choose(y>0,dev1,dev0);
deviance=-2#sum(dev01);
%end;
%if %upcase(&dist)=ZINB %then
%do;
dev0=(log(zk+(1-
zk)#(1+mu#k/weight)##(-1/k)));
dev1=((log(1-
zk)+y#log(mu#k/weight)-
(y+weight/k)#log(1+mu#k/weight)));
dev01=choose(y>0,dev1,dev0);
deviance=-2#sum(dev01);
%end;
scale1=deviance/df;
%if &strata NE | &cluster NE %then
%do;
scale1=n*deviance/((np)#
sum(weight));
%end;
sdeviance=deviance/scale1;
/*Pearson and Scale pearson*/
pearson=j(ncol(y),1,0);
do ii=1 to ncol(y);
pearson[ii,]=sum((weight#(y[,ii]-
mu[,ii])##2)/vmu[,ii]);
end;
pearson=sum(pearson);
scale2=pearson/df;
%if &strata NE | &cluster NE %then
%do;
scale2=n*pearson/((np)#
sum(weight));
%end;
spearson=pearson/scale1;
valeudf=sdeviance/df;
valeudf1=spearson/df;
dev=df||deviance||scale1;
sdev=df||sdeviance||valeudf;
pcs=df||pearson||scale2;
sp=df||spearson||valeudf1;
phi=1;
scaled=1;
%end;
/*Covariance*/
/*Weights*/
do ii=1 to ncol(y);
we[,ii]=weight/(phi#vmu[,ii]#gderiv[,ii]##2);
wo[,ii]=we[,ii]+weight#(y[,ii]-
mu[,ii])#((vmu[,ii]#gderiv2[,ii]+vmuderiv[,ii]#gderiv[,ii])/((vmu[,ii]##2)#(gderiv[,ii]##3)#phi));
end;
/*Hessian matrix*/
do ii=1 to ncol(y);
m1=(ii-1)*ncol(x)+1;
m2=ii*ncol(x);
H[m1:m2, m1:m2]=-phi#(X#we[,ii])`*X;
end;
if nh=0 & ni=0 then
do;
G=phi#G;
end;
/*Q matrix*/
Q=-H;
Sigma=j(ncol(x),ncol(y),0);
do ii=1 to ncol(y);
%if %upcase(&vadjust)=N %then
%do;
Sigma[,ii]=vecdiag(ginv(Q[(ii-
1)*ncol(x)+1:ii*ncol(x), (ii-
1)*ncol(x)+1:ii*ncol(x)])*(G[(ii-
1)*ncol(x)+1:ii*ncol(x), (ii-
1)*ncol(x)+1:ii*ncol(x)])*ginv(Q[(ii-
1)*ncol(x)+1:ii*ncol(x), (ii-1)*ncol(x)+1:ii*ncol(x)]));
%end;
%if
%upcase(&vadjust)=Y %then
%do;
/*Matrix
correction*/
traceqg=trace(ginv(Q[(ii-
1)*ncol(x)+1:ii*ncol(x), (ii-
1)*ncol(x)+1:ii*ncol(x)])*G[(ii-1)*ncol(x)+1:ii*ncol(x),
(ii-1)*ncol(x)+1:ii*ncol(x)]);
traceqgp=traceqg/p;
km=max(&delta,traceqgp);
lambda=min(&fi,
p/(ehi[nrow(ehi),2]-p));
correction=(km#lambda)*ginv(Q[(ii-
1)*ncol(x)+1:ii*ncol(x), (ii-1)*ncol(x)+1:ii*ncol(x)]);
Sigma[,ii]=vecdiag(ginv(Q[(ii-
1)*ncol(x)+1:ii*ncol(x), (ii-
1)*ncol(x)+1:ii*ncol(x)])*G[(ii-1)*ncol(x)+1:ii*ncol(x),
(ii-1)*ncol(x)+1:ii*ncol(x)]*ginv(Q[(ii-
1)*ncol(x)+1:ii*ncol(x), (ii-
1)*ncol(x)+1:ii*ncol(x)])+correction);
%end;
end;
/*Standard Error*/
%if %upcase(&dist)=BINOMIAL %then
%do;
std=sqrt(Sigma);
*std=shapecol(shapecol(std,ncol(x),ncol(y))`,n
col(x)*ncol(y),1);
std=T(shape(T(shape(std`,ncol(y),ncol(x)))`,1,ncol(x)*ncol(y)));
%end;
%else
%do;
std=sqrt(Sigma);
%end;
%if %upcase(&dist)=ZIP or %upcase(&dist)=ZINB %then
%do;
tab=zk||G0||w;
/*logit link function */
eta=exp(G0*lambda);
mul=mi#(eta/(1+eta));
*mul=(1-pi)#eta;
/*Variance function (binomial)*/
vmu=(1/mi)#mul#(1-mul);
/*Firt derivative - variance function*/
vmuderiv=(1/mi)#(1-2#mul);
%if %upcase(&dist)=ZIP %then
%do;
%end;
%if %upcase(&dist)=ZINB %then
%do;
%end;
/*Firt derivative - link function*/
gderiv=mi/(mul#(mi-mul));
/*Second derivative - link function*/
gderiv2=(mi#(2#mul-mi))/((mul#(mi-mul))##2);
we=weight/(phi#vmu#(gderiv##2));
wo=we+weight#(zkmul)#((
vmu#gderiv2+vmuderiv#gderiv)/((vmu##2)#(gderiv##3
)#phi));
H=-phi#(G0#wo)`*G0;
e=zk-mul;
e2=j(nrow(G0), ncol(zk)*ncol(G0),0);
do ii=1 to ncol(zk);
e2[,(ii-
1)*ncol(G0)+1:ii*ncol(G0)]=(phi#we#e[,ii]#gderiv)#G0;
end;
G=j(ncol(G0)*ncol(zk), ncol(G0)*ncol(zk), 0);
/*G matrix*/
*Strata-Cluster;
if nh>0 & ni>0 then
do;
_h1_=j(nrow(h1),1,0);
_h1_[1]=1;
_i_=j(nrow(i),1,0);
_i_[1]=1;
do kk=2 to nrow(h1);
if h1[kk]=h1[kk-1] then
_h1_[kk]=_h1_[kk-1];
else _h1_[kk]=_h1_[kk- 1]+1;
if i[kk]=i[kk-1] then
_i_[kk]=_i_[kk-1];
else _i_[kk]=_i_[kk-1]+1;
end;
tab=tab||_h1_||_i_;
do ii=1 to ncol(zk);
m1=(ii-1)*ncol(G0)+1;
m2=ii*ncol(G0);
ehi=j(1,ncol(G0)+4,1);
ehi[1,3:2+ncol(G0)]=e2[1,(ii-
1)*ncol(G0)+1:ii*ncol(G0)];
ehi[1]=tab[1,ncol(tab)-1];
ehi[2]=tab[1,ncol(tab)];
ehi[ncol(ehi)-1]=fh[1];
do j=2 to nrow(tab);
if
tab[j,ncol(tab)-1]=tab[j-1,ncol(tab)-1] &
tab[j,ncol(tab)]=tab[j-1,ncol(tab)] then
do;
ehi[nrow(ehi),3:2+ncol(G0)]=ehi[nrow(ehi),3:2+
ncol(G0)]+e2[j,(ii-1)*ncol(G0)+1:ii*ncol(G0)];
ehi[nrow(ehi),ncol(ehi)]=ehi[nrow(ehi),ncol(eh
i)]+1;
ehi[nrow(ehi),ncol(ehi)-
1]=(ehi[nrow(ehi),ncol(ehi)-1]+fh[j])/2;
end;
else
ehi=ehi//(tab[j,ncol(tab)-
1]||tab[j,ncol(tab)]||e2[j,(ii-
1)*ncol(G0)+1:ii*ncol(G0)]||fh[j]||j(1,1,1));
end;
ee=ehi[,ncol(ehi)-
1:ncol(ehi)];
ehi=ehi[,1:ncol(ehi)-
2]||j(nrow(ehi),3+ncol(G0),0);
ehi[1,3+ncol(G0):2+2*ncol(G0)]=ehi[1,3:2+ncol(
G0)];
ehi[,ncol(ehi)-
2:ncol(ehi)-1]=ee;
ehi[1,ncol(ehi)]=ehi[1,ncol(ehi)-1];
count=0;
do j=2 to nrow(ehi);
if
ehi[j,1]=ehi[j-1,1] then
do;
ehi[j,3+ncol(G0):2+2*ncol(G0)]=ehi[j-
1,3+ncol(G0):2+2*ncol(G0)]+ehi[j,3:2+ncol(G0)];
count=count+1;
ehi[j,ncol(ehi)]=ehi[j-
1,ncol(ehi)]+ehi[j,ncol(ehi)-1];
ehi[j,ncol(ehi)-2]=ehi[j-1,ncol(ehi)-2];
end;
else
do;
if
ehi[j-1,1]=1 then
do;
ehi[1:count+1,ncol(ehi)-1]=count+1;
ehi[1:count+1,3+ncol(G0):2+2*ncol(G0)]=repeat(
ehi[j-1,3+ncol(G0):2+2*ncol(G0)]/ehi[j-1,ncol(ehi)-
1],count+1);
in=ehi[j,2];
ehi[1:count+1,ncol(ehi)]=repeat(ehi[j-
1,ncol(ehi)],count+1);
end;
else
do;
ehi[in:in+count,ncol(ehi)-1]=count+1;
ehi[in:in+count,3+ncol(G0):2+2*ncol(G0)]=repeat(ehi[j-1,3+ncol(G0):2+2*ncol(G0)]/ehi[j-1,ncol(ehi)-
1],count+1);
ehi[in:in+count,ncol(ehi)]=repeat(ehi[j-
1,ncol(ehi)],count+1);
in=ehi[j,2];
end;
ehi[j,3+ncol(G0):2+2*ncol(G0)]=ehi[j,3:2+ncol(G0)];
ehi[j,ncol(ehi)]=ehi[j,ncol(ehi)-1];
count=0;
end;
if j=nrow(ehi) then
do;
ehi[in:in+count,ncol(ehi)-1]=count+1;
ehi[in:in+count,3+ncol(G0):2+2*ncol(G0)]=repeat(ehi[j,3+ncol(G0):2+2*ncol(G0)]/ehi[j,ncol(ehi)-
1],count+1);
ehi[in:in+count,ncol(ehi)]=repeat(ehi[j,ncol(ehi)],count+1);
in=ehi[j,2];
end;
do jj=1 to nrow(ehi);
if
ehi[jj,ncol(ehi)-1]=1 then
do;
ehi[jj,ncol(ehi)-1]=2;
ehi[jj,3:2+ncol(G0)]=0;
ehi[jj,3+ncol(G0):2+2*ncol(G0)]=0;
end;
end;
end;
G[m1:m2,
m1:m2]=((n-1)/(n-ncol(G0)))*(((ehi[,ncol(ehi)-
1]/(ehi[,ncol(ehi)-1]-1))#(ehi[,3:2+ncol(G0)]-
ehi[,3+ncol(G0):2+2*ncol(G0)]))`*(ehi[,3:2+ncol(G0)]-
ehi[,3+ncol(G0):2+2*ncol(G0)]));
%if
%upcase(&fpc)^= %then
%do;
G[m1:m2,
m1:m2]=((n-1)/(n-ncol(G0)))*((((ehi[,ncol(ehi)-1]#(1-
ehi[,ncol(ehi)-2]))/(ehi[,ncol(ehi)-1]-
1))#(ehi[,3:(2+ncol(G0))]-
ehi[,3+ncol(G0):2+2*ncol(G0)]))`*(ehi[,3:(2+ncol(G0))]-
ehi[,3+ncol(G0):2+2*ncol(G0)]));
%end;
*print G;
end;
df=ehi[nrow(ehi),2]-ehi[nrow(ehi),1];
end;
*Strata;
else if nh>0 & ni=0 then
do;
_h1_=j(nrow(h1),1,0);
_h1_[1]=1;
do kk=2 to nrow(h1);
if h1[kk]=h1[kk-1] then
_h1_[kk]=_h1_[kk-1];
else _h1_[kk]=_h1_[kk- 1]+1;
end;
tab=tab||_h1_;
do ii=1 to ncol(zk);
ehi=j(1,2*ncol(G0)+2,0);
ehi[1,2:1+ncol(G0)]=e2[1,(ii-
1)*ncol(G0)+1:ii*ncol(G0)];
ehi[1]=tab[1,ncol(tab)];
ehi[ncol(ehi)]=1;
do j=2 to nrow(tab);
if tab[j,ncol(tab)]=tab[j-1,ncol(tab)] then
do;
ehi[nrow(ehi),2:1+ncol(G0)]=ehi[nrow(ehi),2:1+
ncol(G0)]+e2[j,(ii-1)*ncol(G0)+1:ii*ncol(G0)];
ehi[nrow(ehi),ncol(ehi)]=ehi[nrow(ehi),ncol(ehi)]+1;
end;
else
do;
ehi=ehi//(tab[j,ncol(tab)]||e2[j,]||j(1,ncol(G0)+1,1));
end;
end;
ehi[,2+ncol(G0):1+2*ncol(G0)]=ehi[,2:1+ncol(G0)]/ehi[,ncol(ehi)];
do jj=1 to nrow(ehi);
ehi2=ehi2//repeat(ehi[jj,2+ncol(G0):2+2*ncol(G0)],ehi[jj,ncol(ehi)]);
end;
eh=e2[,(ii-
1)*ncol(G0)+1:ii*ncol(G0)]||ehi2[,1:ncol(ehi2)-
1]||fh||ehi2[,ncol(ehi2)];
G[m1:m2, m1:m2]=((n-1)/(nncol(
G0)))*(((eh[,ncol(eh)]/(eh[,ncol(eh)]-
1))#(eh[,1:ncol(G0)]-
eh[,1+ncol(G0):2*ncol(G0)]))`*(eh[,1:ncol(G0)]-
eh[,1+ncol(G0):2*ncol(G0)]));
%if %upcase(&fpc) NE %then
%do;
G[m1:m2, m1:m2]=((n-1)/(nncol(
G0)))*((((eh[,ncol(eh)]#(1-eh[,ncol(eh)-
1]))/(eh[,ncol(eh)]-1))#(eh[,1:ncol(G0)]-
eh[,1+ncol(G0):2*ncol(G0)]))`*(eh[,1:ncol(G0)]-
eh[,1+ncol(G0):2*ncol(G0)]));
%end;
*print G;
free ehi2;
end;
df=n-ehi[<>,1];
end;
*Cluster;
else if ni>0 & nh=0 then
do;
_i_=j(nrow(i),1,0);
_i_[1]=1;
do kk=2 to nrow(i);
if i[kk]=i[kk-1] then
_i_[kk]=_i_[kk-1];
else _i_[kk]=_i_[kk-1]+1;
end;
tab=tab||_i_;
do ii=1 to ncol(zk);
m1=(ii-1)*ncol(G0)+1;
m2=ii*ncol(G0);
ehi=j(1,ncol(G0)+3,0);
ehi[1,2:1+ncol(G0)]=e2[1,(ii-
1)*ncol(G0)+1:ii*ncol(G0)];
ehi[1]=tab[1,ncol(tab)];
ehi[ncol(ehi)-
1]=fh[1];
ehi[1,ncol(ehi)]=1;
do j=2 to nrow(tab);
if
tab[j,ncol(tab)]=tab[j-1,ncol(tab)] then
do;
ehi[nrow(ehi),2:1+ncol(G0)]=ehi[nrow(ehi),2:1+ncol(G0)]+
e2[j,(ii-
1)*ncol(G0)+1:ii*ncol(G0)];
ehi[nrow(ehi),ncol(ehi)-
1]=(ehi[nrow(ehi),ncol(ehi)-
1]+fh[j])/2;
ehi[nrow(ehi),ncol(ehi)]=ehi[nrow(ehi),ncol(
ehi)]+1;
end;
else
do;
ehi=ehi//(tab[j,ncol(tab)]||e2[j,(ii-
1)*ncol(G0)+1:ii*ncol(G0)]||fh[j]||j(1,1,1));
end;
end;
G[m1:m2,m1:m2]=((n-1)/(n-ncol(G0)))*(ehi[<>,1]/(ehi[<>,1]-
1))*((ehi[,2:1+ncol(G0)]-
ehi[+,2:1+ncol(G0)]/ehi[+,ncol(ehi)])`*(ehi[,2:1+ncol(G0
)]-ehi[+,2:1+ncol(G0)]/ehi[+,ncol(ehi)]));
%if
%upcase(&fpc) NE %then
%do;
G[m1:m2,
m1:m2]=((n-1)/(n-ncol(G0)))*(ehi[<>,1]/(ehi[<>,1]-
1))*(((1-ehi[,ncol(ehi)-1])#(ehi[,2:1+ncol(G0)]-
ehi[+,2:1+ncol(G0)]/ehi[+,ncol(ehi)]))`*(ehi[,2:1+ncol(G0)]-ehi[+,2:1+ncol(G0)]/ehi[+,ncol(ehi)]));
%end;
*print G;
end;
df=ehi[<>,1]-1;
end;
else
do;
do ii=1 to
ncol(zk);
m1=(ii-
1)*ncol(G0)+1;
m2=ii*ncol(G0);
G[m1:m2,
m1:m2]=-H[(ii-1)*ncol(G0)+1:ii*ncol(G0), (ii-
1)*ncol(G0)+1:ii*ncol(G0)];
end;
df=nrow(zk)-ncol(G0);
end;
/*Covariance*/
/*Weights*/
do ii=1 to ncol(zk);
we[,ii]=weight/(phi#vmu[,ii]#gderiv[,ii]##2);
end;
/*Hessian matrix*/
do ii=1 to ncol(zk);
m1=(ii-1)*ncol(G0)+1;
m2=ii*ncol(G0);
H[m1:m2, m1:m2]=-
phi#(G0#we[,ii])`*G0;
end;
if nh=0 & ni=0 then
do;
G=phi#G;
end;
/*Q matrix*/
Q=-H;
Sigmal=j(ncol(G0),ncol(zk),0);
do ii=1 to ncol(zk);
%if %upcase(&vadjust)=N %then
%do;
Sigmal[,ii]=vecdiag(ginv(Q[(ii-
1)*ncol(G0)+1:ii*ncol(G0), (ii-
1)*ncol(G0)+1:ii*ncol(G0)])*(G[(ii-
1)*ncol(G0)+1:ii*ncol(G0), (ii-
1)*ncol(G0)+1:ii*ncol(G0)])*ginv(Q[(ii-
1)*ncol(G0)+1:ii*ncol(G0), (ii-
1)*ncol(G0)+1:ii*ncol(G0)]));
%end;
%if
%upcase(&vadjust)=Y %then
%do;
/*Matrix
correction*/
traceqg=trace(ginv(Q[(ii-
1)*ncol(G0)+1:ii*ncol(G0), (ii-
1)*ncol(G0)+1:ii*ncol(G0)])*G[(ii-
1)*ncol(G0)+1:ii*ncol(G0), (ii-
1)*ncol(G0)+1:ii*ncol(G0)]);
traceqgp=traceqg/p;
km=max(&delta,traceqgp);
lambda=min(&fi,p/(ehi[nrow(ehi),2]-p));
correction=(km#lambda)*ginv(Q[(ii-
1)*ncol(G0)+1:ii*ncol(G0), (ii-
1)*ncol(G0)+1:ii*ncol(G0)]);
Sigmal[,ii]=vecdiag(ginv(Q[(ii-
1)*ncol(G0)+1:ii*ncol(G0), (ii-
1)*ncol(G0)+1:ii*ncol(G0)])*G[(ii-
1)*ncol(G0)+1:ii*ncol(G0), (ii-
1)*ncol(G0)+1:ii*ncol(G0)]*ginv(Q[(ii-
1)*ncol(G0)+1:ii*ncol(G0), (ii-
1)*ncol(G0)+1:ii*ncol(G0)])+correction);
%end;
end;
stdlambda=sqrt(abs(Sigmal));
waldlambda=(lambda/stdlambda)##2;
pvaluelambda=1-probchi(waldlambda,1);
%end;
%if %upcase(&dist)=GAMMA %then
%do;
lli=(weight/phi)#log(weight#yc/(phi#mu))-
(weight#yc/(phi#mu))-log(yc)-log(gamma(weight/phi));
*lli=li;
%end;
%if %upcase(&dist)=NORMAL %then
%do;
pi=constant("pi");
logphiweight=j(nrow(weight),1,0);
do jj=1 to nrow(weight);
if weight[jj]=0 then
logphiweight[jj]=0;
else logphiweight[jj]=log(phi/weight[jj]);
end;
lli=-0.5#((weight#((yc-mu)##2))/phi + logphiweight +log(2#pi));
*lli=li;
%end;
%if %upcase(&dist)=POISSON %then
%do;
*li=(weight/phi)#(y#log(mu)-mu);
fy=fact(y);
lli=weight#(y#log(mu)-mu-log(fy));
%end;
%if %upcase(&dist)=NEGBIN %then
%do;
gamma1=y+weight/k;
gamma1=choose(gamma1>170,100,gamma1);
gamma1=choose(gamma1=0,0.0001,gamma1);
gamma1=gamma(gamma1);
gamma2=(y+1);
gamma2=choose(gamma2>170,100,gamma2);
gamma2=choose(gamma2=0,0.0001,gamma2);
gamma2=gamma(gamma2);
gamma3=weight/k;
gamma3=choose(gamma3>170,100,gamma3);
gamma3=choose(gamma3=0,0.0001,gamma3);
gamma3=gamma(gamma3);
lli=y#log(k*mu/weight)-(y+weight/k)#log(1+k*mu/weight)+log(gamma1/(gamma2#gamma3));
%end;
%if %upcase(&dist)=IG %then
%do;
pi=constant('pi');
lli=-0.5#((weight#((ycmu)##2)/(phi#yc#(mu##2)))+log((phi#(yc##3))/weight)+log(2#pi));
*lli=li;
%end;
%if %upcase(&dist)=BINOMIAL %then
%do;
xbetas=0;
yxbetas=0;
do ii=1 to ncol(y);
*li=(weight/phi)#(y#log(mu)+(miy)#log(1-mu));
*lli=weight#(log(comb(mi,y))+y#log(mu)+(miy)#log(1-mu));
xbetas=xbetas+exp(X*b1[,ii]);
yxbetas=yxbetas+y[,ii]#(X*b1[,ii]);
end;
lli=sum(weight#(yxbetaslog(1+xbetas)));
%end;
%if %upcase(&dist)=ZIP %then
%do;
ll0=(weight#log(zk+(1-zk)#exp(-mu)));
ll1=(weight#(log(1-zk)+y#log(mu)-mulog(fact(y))));
lli=choose(y>0,ll1,ll0);
%end;
%if %upcase(&dist)=ZINB %then
%do;
ll0=(log(zk+(1-zk)#(1+mu#k/weight)##(-1/k)));
ll1=((log(1-zk)+y#log(mu#k/weight)-
(y+weight/k)#log(1+mu#k/weight)));
lli=choose(y>0,ll1,ll0);
%end;
*l=sum(li);
p=nrow(b);
ll=-2#sum(lli);
AIC=.||ll+2#p||.;
AICC=.||ll+2#p#(n/(n-p-1))||.;
BIC=.||ll+p#log(n)||.;
*logl=.||l||.;
_2flogl=.||ll||.;
Report=dev//sdev//pcs//sp//_2flogl//AIC//AICC//BIC;
ReportName={ "Deviance" "Scaled Deviance" "Pearson Chi-Square" "Scaled Pearson X2" "-2 Log Likelihood" "AIC" "AICC" "BIC"};
ReportcolName={ "DF" "Value" "Value/DF"};
%if %upcase(&dist)=BINOMIAL %then
%do;
AIC=ll+2#p;
AICC=ll+2#p#(n/(n-p-1));
BIC=ll+p#log(n);
_2flogl=ll;
Report=_2flogl//AIC//AICC//BIC;
ReportName={"-2 Log Likelihood" "AIC" "AICC" "BIC"};
ReportcolName={ "Value"};
%end;
print Report[label="Criteria For Assessing Goodness Of Fit" rowname=reportname
colname=reportcolname format=12.4];
/*t test and Wald test*/
t=b/std;
pvalue=2#(1-probt(abs(t),df));
%if %upcase(&dist)=BINOMIAL %then
%do;
wald=(b/std)##2;
pvalue=1-probchi(wald, 1);
t=wald;
%end;
*b1=b`||scaled;
*b1=shapecol(b,0,1)`||scaled;
%if %upcase(&dist)=NEGBIN or
%upcase(&dist)=ZINB %then
%do;
scaled=k;
tk=k/sek;
pvaluek=2#(1-probt(abs(tk),df));
std=std//sek;
t=t//tk;
pvalue=pvalue//pvaluek;
%end;
b1=T(shape(b`,1,0))`||scaled;
b2=b1`;
if ncol(y)>1 then
do;
varname1=varname`||repeat(varname`,ncol(y)-1);
varname1=rowvec(varname1)`//"Scale";
%if %upcase(&dist)=NEGBIN or
%upcase(&dist)=ZINB %then
%do;varname1=rowvec(varname1)`//"Dispersion";%end;
end; else do;
varname1=varname`//"Scale";
%if %upcase(&dist)=NEGBIN or %upcase(&dist)=ZINB %then
%do;varname1=rowvec(varname1)`//"Dispersion";%end;
end;
%if %upcase(&intercept)=Y %then %do;
if ncol(y)>1 then do;
varname1="Intercept"//varname`||repeat("Intercept"//varname`,ncol(y)-1);
varname1=rowvec(varname1)`//"Scale";
%if %upcase(&dist)=NEGBIN or
%upcase(&dist)=ZINB %then
%do;varname1=rowvec(varname1)`//"Dispersion";%end;
end; else do;
varname1="Intercept"//varname`//"Scale";
%if %upcase(&dist)=NEGBIN or %upcase(&dist)=ZINB %then
%do;varname1="Intercept"//varname`//"Dispersion";%end;
end;
%end;
print "Analysis of Maximum Likelihood Parameter Estimates";
%if %upcase(&dist)=BINOMIAL %then %do;
varnamey=repeat(freq[2:nrow(freq),2],ncol(x));
print varname1[label="Parameter"]
varnamey[label="&y"] b2[label="Estimate" format=12.4]
std[label="Standard Error" format=12.4] t[label="Wald Chi-Square" format=12.2]
pvalue[label="Pr > ChiSq"
format=pvalue6.4];
%end; %else %do;
print varname1[label="Parameter"]
b2[label="Estimate" format=12.4] std[label="Standard Error" format=12.4] t[label="t Value" format=12.2]
pvalue[label="Pr > |t|"
format=pvalue6.4];
%end;
%if %upcase(&dist) NE BINOMIAL %then
%do;
print "Note: The denominator degrees of freedom for the t tests is" df[label=""];
%end;
%if %upcase(&scale)=DEVIANCE %then %do;
%if %upcase(&dist)=GAMMA %then %do;
print "Note: The Gamma scale parameter was estimated by DOF/DEVIANCE.";
%end; %if %upcase(&dist)=NORMAL |
%upcase(&dist)=POISSON | %upcase(&dist)=IG |
%upcase(&dist)=BINOMIAL %then %do;
print "Note: The scale parameter was estimated by the square root of DEVIANCE/DOF.";
%end;
%end; %if %upcase(&scale)=PEARSON %then %do;
%if %upcase(&dist)=GAMMA %then %do;
print "Note: The Gamma scale parameter was estimated by DOF/Pearson's Chi-Square.";
%end; %if %upcase(&dist)=NORMAL |
%upcase(&dist)=POISSON | %upcase(&dist)=IG |
%upcase(&dist)=BINOMIAL %then %do;
print "Note: The scale parameter was estimated by the square root of Pearson's Chi-Square/DOF. ";
%end;
%end;
%if %upcase(&dist)=ZIP or %upcase(&dist)=ZINB
%then %do;
print "Analysis Of Maximum Likelihood Zero Inflation Parameter Estimate";
%if &xzip ne %then
%do;
varnamezip1="Intercept"//varnamezip`;
%end;
%else
%do;
varnamezip1={"Intercept"};
%end;
print varnamezip1[label="Parameter"]
lambda[label="Estimate" format=12.4]
stdlambda[label="Standard Error" format=12.4]
waldlambda[label="Wald Chi-Square"
format=12.2] pvaluelambda[label="Pr > ChiSq"
format=pvalue6.4];;
%end;
/*Odds Ratio*/
%if %upcase(&dist)=BINOMIAL %then
%do;
if ncol(y)>1 then
do;
varnameodds=rowvec(varname`||repeat(varname`,ncol(y)-1))`;
varnameoddsy=repeat(freq[2:nrow(freq),2],ncol(x));
end;
else
do;
varnameodds=varname`;
varnameoddsy=repeat(freq[2:nrow(freq),2],ncol(x));
end;
%if %upcase(&intercept)=Y %then
%do;
if ncol(y)>1 then
do;
varnameodds=rowvec(varname`||repeat(varname`,ncol(y)-1))`;
varnameoddsy=repeat(freq[2:nrow(freq),2],ncol(x)-1);
end;
else
do;
varnameodds=varname`;
varnameoddsy=repeat(freq[2:nrow(freq),2],ncol(x)-1);
end;
%end;
%end;
%if %upcase(&dist)=BINOMIAL %then
%do;
*beta_=shapecol(b,0,1);
beta_=T(shape(b`,1,0));
odds_ratio=exp(beta_);
odds_ratioclm=j(nrow(odds_ratio),2,0);
odds_ratioclm[,1]=exp(beta_[1:nrow(beta_)]-
probit(1-&alpha/2)*std[1:nrow(std)]);
odds_ratioclm[,2]=exp(beta_[1:nrow(beta_)]+probit(1-&alpha/2)*std[1:nrow(std)]);
%if %upcase(&intercept)=Y %then
%do;
odds_ratio=exp(beta_[ncol(y)+1:nrow(beta_)]);
odds_ratioclm=j(nrow(odds_ratio),2,0);
odds_ratioclm[,1]=exp(beta_[ncol(y)+1:nrow(beta_)]-probit(1-&alpha/2)*std[ncol(y)+1:nrow(std)]);
odds_ratioclm[,2]=exp(beta_[ncol(y)+1:nrow(beta_)]+probit(1-&alpha/2)*std[ncol(y)+1:nrow(std)]);
%end;
print "Odds Ratio Estimates";
print varnameodds[colname='Effect' label=""]
varnameoddsy[colname="&y" label=""] odds_ratio[label=""
colname='Point Estimate' format=12.3]
odds_ratioclm[label="" colname={"Wald %sysevalf(100*(1-&alpha))% CL Lower" "Wald %sysevalf(100*(1-&alpha))% CL Upper"} format=12.3];
%end;
quit;
%end;
%mend surveygenmod;
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment