option nocenter ls=72; %macro qif(data=, yvar=, xvar=, id=, weight=1, dist=, corr=, descend=Y, print=Y, outpar=, outqif=, outcov=,outres= ); %macro outprint(print); %if &print^=N %then %do; title "Information Criteria of QIF model"; proc print data = par_qif; id test; run; title "Estimates Of The Model Parameters"; proc print data = par_est; id parameter; run; title "Covariance Matrix Of The Model Parameters"; proc print data = par_cov; id parameter; run; %end; %mend outprint; %let stop=0; %if (%length(&data)=0) %then %do; %put ERROR: DATA is not provided.; %let stop = 1; %end; %if (%length(&yvar)=0) %then %do; %put ERROR: Response variable YVAR is not provided.; %let stop = 1; %end; %if (%length(&xvar)=0) %then %do; %put ERROR: Covariate variable XVAR is not provided.; %let stop = 1; %end; %if (%length(&id)=0) %then %do; %put ERROR: ID is not provided.; %let stop = 1; %end; %let dist=%upcase(&dist); %if (&dist^=NORMAL and &dist^=GAMMA and &dist^=POISSON and &dist^=BIN) %then %do; %put ERROR: Distribution variable DIST must be one of NORMAL, POISSON, GAMMA or BIN.; %let stop = 1; %end; %let corr=%upcase(&corr); %if (&corr^=AR and &corr^=EXCH and &corr^=IND and &corr^=UNSTR) %then %do; %put ERROR: Correlation structure variable CORR must be one of AR, EXCH, IND or UNSTR.; %let stop = 1; %end; %if &stop >0 %then %do; %put NOTE: Macro SLMM stopped because of errors!; %goto done; %end; %let print=%upcase(&print); ods listing close; %if (&dist=BIN and %upcase(&descend=Y)) %then %do; proc genmod data=&data descend; class &id; model &yvar=&xvar / d=&dist; repeated subject = &id / type=&corr; ods output GEEEmpPEst = GEEEmpPEst; run; %end; %else %do; proc genmod data=&data; class &id; model &yvar=&xvar / d=&dist; repeated subject = &id / type=&corr; ods output GEEEmpPEst = GEEEmpPEst; run; %end; ods listing; proc iml; start getname; parameter = j(np,1,char(20)); parameter[1,1] = "intercept"; do i=1 to (np-1); parameter[i+1,1] = scan("&xvar", i); end; finish; start outinf; %if &print^=N %then %do; col_name = {"Iteration","Statistic","DF", "Pr>Q","AIC","BIC"}; par_qif = iteration || Q || DF || qifpvalue || AIC || BIC; create par_qif from par_qif[rowname=Test colname=col_name]; append from par_qif[rowname=Test]; col_name = {"Estimate", "Standard Error", "Z", "Pr>|Z|"}; par_est = beta || parse || Z || parpvalue; create par_est from par_est[rowname=parameter colname=col_name]; append from par_est[rowname=parameter]; create par_cov from invarqif2dev[rowname=parameter colname=parameter]; append from invarqif2dev[rowname=parameter]; %end; finish; start outdata; qifpvalue= 1 - probchi(Q,np); parvar = vecdiag(invarqif2dev); parse = sqrt(parvar); Z = beta/parse; parpvalue = 2*(1-probnorm(abs(Z))); Test = "QIF"; DF = np; Statistic = Q; %if &corr=IND %then %do; AIC = .; BIC = .; %end; %else %do; AIC = Q+2*np; BIC = Q+np*log(nsub); *BIC = Q+np*log(nrow(y)); %end; %if %length(&outqif)>0 %then %do; create &outqif var {Test Iteration Statistic DF qifpvalue AIC BIC}; append; %end; %if %length(&outpar)>0 %then %do; estimate = beta; stderr = parse; create &outpar var {parameter estimate stderr Z parpvalue}; append; %end; %if %length(&outcov)>0 %then %do; create &outcov from invarqif2dev[colname=parameter]; append from invarqif2dev; %end; finish; start outresidual; %if %length(&outres)>0 %then %do; %if &dist = NORMAL %then %do; u = x*beta; pear_res = y - u; devi_res = sign(y-u)#abs(y-u); %end; %else %if &dist = POISSON %then %do; u = exp(x*beta); pear_res = (y - u)/sqrt(u); devi_res = j(nrow(y),1,.); do i = 1 to nrow(y); if y[i,1]=0 then devi_res[i,1] = 2*u[i,1]; else devi_res[i,1] = 2*(y[i,1]*log(y[i,1]/u[i,1])-(y[i,1]-u[i,1])); end; devi_res = sign(y-u)#sqrt(devi_res); %end; %else %if &dist = GAMMA %then %do; u = 1/(x*beta); pear_res = (y - u)/u; devi_res = j(nrow(y),1,.); do i = 1 to nrow(y); if y[i,1]^=0 then devi_res[i,1] = 2*(-log(y[i,1]/u[i,1])+y[i,1]/u[i,1]-1); end; devi_res = sign(y-u)#sqrt(devi_res); %end; %else %if &dist = BIN %then %do; u = 1/(1+exp(-x*beta)); pear_res = (y - u)/sqrt(u#(1-u)); devi_res = j(nrow(y),1,.); do i = 1 to nrow(y); if y[i,1]=0 then devi_res[i,1] = 2*(log(1/(1-u[i,1]))); else devi_res[i,1] = 2*(log(1/u[i,1])); end; devi_res = sign(y-u)#sqrt(devi_res); %end; col_name={"id","y","uhat","pearson residuals","deviance residuals"}; residual = id || y || u || pear_res || devi_res; create &outres from residual[colname=col_name]; append from residual; %end; finish; /******* main QIF program ******/ use GEEEmpPEst; read all var{estimate} into beta; use &data; read all var{&xvar} into x; read all var{&yvar} into y; read all var{&id} into id; %if &weight=1 %then %do; wm = j(nrow(y),1,1); %end; %else %do; read all var{&weight} into wm; %end; /*** deal with missing data by DELETE under MCAR* ***/ xyncol = ncol(x)+3; xy = j(nrow(x),xyncol,0); xy[,1]=y; xy[,2:(xyncol-2)] = x; xy[,xyncol-1] = wm; xy[,xyncol] = id; nxy = nrow(xy); print (nxy); do i = 1 to xyncol; xy = xy[loc(xy[,i] ^= .),]; end; nxym = nrow(xy); print (nxym); y=j(nrow(xy),1,0); x=j(nrow(xy),ncol(xy)-2,0); wm = j(nrow(xy),1,0); id = j(nrow(xy),1,0); y=xy[,1]; x=xy[,1:(xyncol-2)]; x[,1] = 1; wm=xy[,xyncol-1]; id = xy[,xyncol]; /* end missing data process */ /* xx=j(nrow(x),ncol(x)+1,0); xx[,1]=1; xx[,2:ncol(x)+1]=x; x=xx; */ np=ncol(x); uid = unique(id); nsub = ncol(uid); betadiff = 1; iteration=0; diff = 1e-8; betanew=beta; do until(betadiff <= diff | iteration > 1000); beta = betanew; %if &corr = IND %then %do; sumg=j(np,1,0); sumc=j(np,np,0); arsumg=j(np,1,0); arsumc=j(np,np,0); gi=j(np,1,0); arsumgfirstdev=j(np,np,0); firstdev=j(np,np,0); %end; %else %do; sumg=j(2*np,1,0); sumc=j(2*np,2*np,0); arsumg=j(2*np,1,0); arsumc=j(2*np,2*np,0); gi=j(2*np,1,0); arsumgfirstdev=j(2*np,np,0); firstdev=j(2*np,np,0); %end; %if &corr = UNSTR %then %do; *print (yi); xi = x[loc(id=uid[1,1]),]; ni=nrow(xi); m0=I(ni); m1=j(ni,ni,0); do i=1 to nsub; xi = x[loc(id=uid[1,i]),]; yi = y[loc(id=uid[1,i])]; %if &dist = NORMAL %then %do; ui = xi*beta; %end; %else %if &dist = POISSON %then %do; ui = exp( (xi*beta) ); %end; %else %if &dist = GAMMA %then %do; ui = 1/(xi*beta); %end; %else %if &dist = BIN %then %do; ui = 1/(1+exp(-xi*beta)); %end; m1= m1 + (yi-ui)*t(yi-ui); end; m1=1/nsub * m1; %end; do i = 1 to nsub; xi = x[loc(id=uid[1,i]),]; yi = y[loc(id=uid[1,i])]; ni = nrow(xi); wmi = diag(wm[loc(id=uid[1,i])]); m0=I(ni); %if &corr = AR %then %do; m1=j(ni,ni,0); do k = 1 to ni; do l= 1 to ni; if abs(k-l)=1 then m1[k,l]=1; end; end; %end; %else %if &corr = EXCH %then %do; m1=j(ni,ni,1)-m0; %end; %else %if &corr = IND %then %do; m1 = j(ni,ni,0); %end; %if &dist = NORMAL %then %do; ui = xi*beta; fui=ui; fui_dev = diag(repeat(1,ni)); vui = diag(repeat(1,ni)); %end; %else %if &dist = POISSON %then %do; ui = exp( (xi*beta) ); fui = log(ui); fui_dev = diag(ui); vui = diag(sqrt(1/ui)); %end; %else %if &dist = GAMMA %then %do; ui = 1/(xi*beta); fui = 1/ui; fui_dev = -diag(ui)*diag(ui); vui = diag(1/ui); %end; %else %if &dist = BIN %then %do; ui = 1/(1+exp(-xi*beta)); fui = log(ui)-log(1-ui); fui_dev = diag(ui)*diag(1-ui); vui = diag(sqrt(1/ui))*diag(sqrt(1/(1-ui))); %end; %if &corr=IND %then %do; wi = t(xi)*fui_dev*vui*m0*vui*wmi; gi0 = (1/nsub)*wi*(yi-ui); gi[1:np,]=gi0; arsumc=arsumc+gi*t(gi); arsumg=arsumg+gi; di0 = -(1/nsub)*wi*fui_dev*xi; firstdev[1:np,]=di0; arsumgfirstdev=arsumgfirstdev+firstdev; %end; %else %if &corr=UNSTR %then %do; wi = t(xi)*fui_dev*m0*wmi; zi = t(xi)*fui_dev*m1*wmi; gi0 = (1/nsub)*wi*(yi-ui); gi1 = (1/nsub)*zi*(yi-ui); gi[1:np,]=gi0; gi[np+1:2*np,]=gi1; arsumc=arsumc+gi*t(gi); arsumg=arsumg+gi; di0 = -(1/nsub)*wi*fui_dev*xi; di1 = -(1/nsub)*zi*fui_dev*xi; firstdev[1:np,]=di0; firstdev[np+1:2*np,]=di1; arsumgfirstdev=arsumgfirstdev+firstdev; %end; %else %do; wi = t(xi)*fui_dev*vui*m0*vui*wmi; zi = t(xi)*fui_dev*vui*m1*vui*wmi; gi0 = (1/nsub)*wi*(yi-ui); gi1 = (1/nsub)*zi*(yi-ui); gi[1:np,]=gi0; gi[np+1:2*np,]=gi1; arsumc=arsumc+gi*t(gi); arsumg=arsumg+gi; di0 = -(1/nsub)*wi*fui_dev*xi; di1 = -(1/nsub)*zi*fui_dev*xi; firstdev[1:np,]=di0; firstdev[np+1:2*np,]=di1; arsumgfirstdev=arsumgfirstdev+firstdev; %end; end; arcinv=ginv(arsumc); Q=t(arsumg)*arcinv*arsumg; arqif1dev=t(arsumgfirstdev)*arcinv*arsumg; arqif2dev=t(arsumgfirstdev)*arcinv*arsumgfirstdev; invarqif2dev=ginv(arqif2dev); betanew=beta-invarqif2dev*arqif1dev; betadiff=sum(abs(betanew-beta)); iteration=iteration+1; end; run getname; run outresidual; run outdata; run outinf; quit; %outprint(&print); %done: %mend qif;