%Macro GEEGropOutMissing(data,id,response,complete,Model_covariate,corr_type,dist,offset,Mdy,Mdx); proc iml; use &data; %let VarListF = %cmpres(&Model_covariate) | ; %let AString = ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789_*|#; %Let n_model = 0; %let VarRead = ; %let VarTest = ; %let M_ModelXP =; %let Max_m_p = 0; **************************************************************************************************; %do %until(%index(&VarListF, | ) < 1); %let n_model = %eval(&n_model + 1); %let VarList1=%substr(&VarListF,1,%index(&VarListF, | )-1); %let VarList2 = &VarList1; %let VarList3 = &VarList1 _; %do %until(%index(&VarList2,* ) < 1); %if %index(&VarList2,* ) >= 1 %then %do; %let VarList2=%substr(&VarList2,1,%index(&VarList2, * )-1) %substr(&VarList2,%index(&VarList2, * )+1); %end; %end; %let VarList2 = %cmpres(&VarList2) _ ; %do %until( %verify(&VarList2, &AString)<1); %let pos = %verify(&VarList2, &AString); %if %index(&VarRead,%substr(&VarList2,1,&pos-1)) <= 0 %then %do; %Let VarRead = &VarRead %substr(&VarList2,1,&pos-1); %end; %let VarList2 = %substr(&VarList2,&pos+1); %end; %do %until(%index(&VarList1,* ) < 1); %if %index(&VarList1,* ) >= 1 %then %do; %let VarList1=%substr(&VarList1,1,%index(&VarList1, * )-1)#%substr(&VarList1,%index(&VarList1, * )+1); %end; %end; %if %verify(&VarList1, &AString) > 0 %then %do; %do %until( %verify(&VarList1, &AString)<1); %let pos = %verify(&VarList1, &AString); %if %index(%substr(&VarList1,&pos-1,1),#) <= 0 and %index(%substr(&VarList1,&pos+1,1),#) <= 0 %then %do; %let VarList1=%substr(&VarList1,1,&pos-1)||%substr(&VarList1,&pos+1); %end; %else %do; %let VarList1=%substr(&VarList1,1,&pos-1)%substr(&VarList1,&pos+1); %end; %end; %end; %let Xmatric&n_model = &VarList1; %let VarList1 = %cmpres(&VarList1)|| ; %let cur_m_p = 1; %do %until(%index(&VarList1,|| ) < 1); %if %index(&VarList1,|| ) >= 1 %then %do; %let cur_m_p = %eval(&cur_m_p + 1); %let VarList1=%substr(&VarList1,%index(&VarList1, || )+2); %end; %end; %if &cur_m_p > &Max_m_p %then %do; %let Max_m_p = &cur_m_p; %let Max_P_model_num = &n_model; %end; %let XmLabel&n_model= %substr(1 1,2,1) ; %do %until( %verify(&VarList3, &AString)<1); %let pos = %verify(&VarList3, &AString); %if %index(%substr(&VarList3,&pos-1,1),*) <= 0 and %index(%substr(&VarList3,1,1),*) <= 0 %then %do; %let XmLabel&n_model= &&XmLabel&n_model %substr(&VarList3,1,&pos-1); %end; %if %index(%substr(&VarList3,1,1),*) > 0 and %index(%substr(&VarList3,&pos-1,1),*) <= 0 %then %do; %let XmLabel&n_model= &&XmLabel&n_model%substr(&VarList3,1,&pos-1); %end; %if %index(%substr(&VarList3,1,1),*) > 0 and %index(%substr(&VarList3,&pos-1,1),*) > 0 %then %do; %let XmLabel&n_model= &&XmLabel&n_model%substr(&VarList3,1,&pos-2); %end; %if %index(%substr(&VarList3,&pos-1,1),*) > 0 %then %do; %let VarList3= *%substr(&VarList3,&pos+1); %end; %else %do; %let VarList3= %substr(&VarList3,&pos+1); %end; %end; %let VarListF=%substr(&VarListF,%index(&VarListF, | ) + 1); %end; **************************************************************************************************; %let tempmodel = &Xmatric1; %let tempmodellabel = &XmLabel1; %let Xmatric1 = &&Xmatric&Max_P_model_num; %let XmLabel1 = &&XmLabel&Max_P_model_num; %let Xmatric&Max_P_model_num = &tempmodel; %let XmLabel&Max_P_model_num = &tempmodellabel; %Let n_corr = 0; %let CorrListF = %cmpres(&corr_type) _ ; %do %until( %verify(&CorrListF, &AString)<1); %Let n_corr = %eval(&n_corr + 1); %let pos = %verify(&CorrListF, &AString); %let corr&n_corr = %lowcase("%substr(&CorrListF,1,&pos-1)"); %let CorrListF = %substr(&CorrListF,&pos+1); %end; %let MdxP = &Mdx; %let VarTest = %cmpres(&Mdx) ; %if %length(&Mdx) > 0 %then %do; %do %until(%index(&Mdx,* ) < 1); %if %index(&Mdx,* ) >= 1 %then %do; %let Mdx=%substr(&Mdx,1,%index(&Mdx, * )-1) %substr(&Mdx,%index(&Mdx, * )+1); %end; %end; %put *********************Mdx &Mdx; %do %until(%index(&VarTest,* ) < 1); %if %index(&VarTest,* ) >= 1 %then %do; %let VarTest=%substr(&VarTest,1,%index(&VarTest, * )-1)#%substr(&VarTest,%index(&VarTest, * )+1); %end; %end; %put *********************VarTest &VarTest; %if %verify(&VarTest, &AString) > 0 %then %do; %do %until( %verify(&VarTest, &AString)<1); %let pos = %verify(&VarTest, &AString); %if %index(%substr(&VarTest,&pos-1,1),#) <= 0 and %index(%substr(&VarTest,&pos+1,1),#) <= 0 %then %do; %let VarTest=%substr(&VarTest,1,&pos-1)||%substr(&VarTest,&pos+1); %end; %else %do; %let VarTest=%substr(&VarTest,1,&pos-1)%substr(&VarTest,&pos+1); %end; %end; %end; %put *********************VarTest &VarTest; %end; MParameter = "Intercept"; %if (%length(&Mdy) = 0 or &Mdy = 0) %then %do; %let M_ModelYP =;%end; %else %do; %let M_ModelYP =&response(t-&Mdy) ... &response(t-1); %do i=&Mdy %to 1 %by-1; MPa="&response(t-&i)"; %if i=1 %then %do; MParameter=MParameter//Mpa; %end; %else %do; MParameter = MParameter//MPa; %end; %end; %end; /* %if %lowcase(&Mdy) = y1 %then %do; MParameter = MParameter//"&response(t-1)";%let M_ModelYP = &response(t-1);%end; %if %lowcase(&Mdy) = y2 %then %do; MParameter = MParameter//"&response(t-1)"//"&response(t-2)";%let M_ModelYP = &response(t-1) &response(t-2);%end; */ %let OutputMPLab = %cmpres(&MdxP) _ ; %do %until(%verify(&&OutputMPLab, &AString) < 1); %if %verify(&&OutputMPLab, &AString) >= 1 %then %do; %let pos = %verify(&&OutputMPLab, &AString); nowParameter = "%substr(&&OutputMPLab,1,&pos-1)"; MParameter = MParameter//nowParameter; %let M_ModelXP = &M_ModelXP %substr(&&OutputMPLab,1,&pos-1); %let OutputMPLab = %substr(&&OutputMPLab,&pos+1); %end; %end; read all var { &id &response &complete &VarRead &offset &Mdx}; idd=unique(&id); n=ncol(idd); nt=nrow(&id); do i=1 to n; ipp=i*(&id=idd[i]); if i=1 then ip=ipp; else ip=ip+ipp; end; %if %length(&complete) > 0 %then %do; r = &complete; %end; start cri; do i=1 to n; ind=loc(ip=i); ni=ncol(ind); yi=y[ind]; xi=x[ind,]; ri=r[ind]; pii=pi[ind]; wi=(ri/pii); mui0=mu0[ind]; offi=off[ind]; tsci=tsc[,i]; etai=xi*b+offi; if d='nor' then do; mui=etai; di=i(ni); ai=di*ophi; end; if d='po' then do; mui=exp(etai); di=diag(mui);ai=di*ophi;end; if d='bin' then do; mui=1/(1+exp(-etai)); di=diag(mui#(1-mui));ai=di*ophi; end; if wc='c' then cori=(1-orho)*i(ni)+j(ni,ni,orho); if wc='i' then cori=i(ni); if wc='ar' then cori=toeplitz(orho##(0:ni-1)); vi=ai##(0.5)*cori*ai##(0.5); /* sci=xi`*di*inv(vi)*(wi#(yi-mui));*/ resi=wi#(yi-mui0); rri=resi*resi`; omegai=j(ni,ni,pii[1]); do jj=2 to ni; omegai[jj:ni, jj:ni]=j(ni-jj+1,ni-jj+1,pii[jj]); end; cr1i=(yi-mui)`*(wi#(yi-mui)); ccr2i=xi`*di*inv(vi)*(omegai#rri)*di*xi; cr2i=xi`*di*inv(vi)*(rri)*di*xi-tsci*resi`*di*xi; hti=xi`*di*ginv(vi)*(wi#vecdiag(di)#xi); if i=1 then do; cr1=cr1i; cr2=cr2i; ccr2=ccr2i; ht=hti;end; else do; cr1=cr1+cr1i; cr2=cr2+cr2i; ccr2=ccr2+ccr2i; ht=ht+hti;end; end; ccri0=cr1+2*trace(inv(ht)*ccr2); cri0=cr1+2*trace(inv(ht)*cr2); finish; start gee; run estsc; oldb=b+1; do iter=1 to 20 while (max(abs(/*oldb-b*/es))>0.00001 & max(abs(es))^=1E10); oldb=b; dh=abs(det(h));/*print dh;*/ if dh<=0.00000001 then es=1E10; else do; delta=-solve(H,es); b=b+delta; run estsc; /*print b es iter;*/ end; end; finish; start estsc; do i=1 to n; ind=loc(ip=i); ni=ncol(ind); yi=y[ind]; xi=x[ind,]; ri=r[ind]; pii=pi[ind]; wi=(ri/pii); offi=off[ind]; etai=xi*b+offi; if ite=1 then ophi=1; if d='nor' then do; mui=etai; di=i(ni); ai=di*(ophi);end; if d='po' then do; mui=exp(etai); di=diag(mui);ai=di*ophi;end; if d='bin' then do; mui=1/(1+exp(-etai)); ai=diag(mui#(1-mui));di=ai; end; if ite=1 then cori=i(ni); else do; if wc='c' then cori=(1-orho)*i(ni)+j(ni,ni,orho); if wc='i' then cori=i(ni); if wc='ar' then cori=toeplitz(orho##(0:(ni-1))); end; vi=ai##(0.5)*cori*ai##(0.5); sci=xi`*di*inv(vi)*(wi#(yi-mui)); hi=-xi`*di*inv(vi)*(wi#vecdiag(di)#xi); resi=wi#(yi-mui)/(vecdiag(sqrt(ai))); rri=resi*resi`; if ni=1 then do; srri=0;assri=0;end; if ni>1 then do; omegai=j(ni,ni,pii[1]); do jj=2 to ni; omegai[jj:ni, jj:ni]=j(ni-jj+1,ni-jj+1,pii[jj]); end; orri=omegai#rri; srri=sum(orri)-/*resi`*(pii#resi)*/trace(orri); asrri=(vecdiag(orri[1:ni-1,2:ni]))[+]; end; ki=ni*(ni-1); phii=((yi-mui)#wi#(yi-mui)/(vecdiag(di)))[+]; if i=1 then do; sc=sci; h=hi; srr=srri; k=ki; asrr=asrri; mu=mui; phi=phii;res=resi; end; else do; sc=sc||sci; h=h+hi; srr=srr+srri; k=k+ki; asrr=asrr+asrri; mu=mu//mui; phi=phi+phii;res=res//resi; end; es=sc[,+]; end; phi=phi/(nt-p); if wc='i' then rho=0; if wc='c' then rho=srr/(k-p); if wc='ar' then rho=asrr/(nt-n-p); finish; start lik; run estscp; oldbp=bp+1; do iterp=1 to 20 while (max(abs(/*oldbp-bp*/esp))>0.00001 & max(abs(esp))^=1E10); oldbp=bp; dh=abs(det(hp));/*print dh;*/ if dh<=0.00000001 then esp=1E10; else do; delta=-solve(Hp,esp); bp=bp+delta; run estscp;/*print bp esp iterp;*/ end; end; finish; start estscp; do i=1 to n; ind=loc(ip=i); ni=ncol(ind); yi=y[ind]; ri=r[ind]; do j=1 to ni; if j=1 then do; piij=1; scpij=0;hpij=0; end; if (j>=2) then do; if mpi=0 then zij=1; else do; if j>mpi then zij=1||(yi[j-mpi:j-1])`;else zij=1||j(1,mpi-j+1,0)||(yi[1:j-1])`; end; /* if mpi='y1' then zij=1||yi[j-1]; if mpi='y2' then do; if j>2 then zij=1||(yi[j-2:j-1])`;else zij=1||0||yi[j-1];end; if mpi='y3' then do; if j>2 then zij=1||(yi[j-2:j-1])`;else zij=1||0||yi[j-1];end; */ if mpix=1 then do; pxi=px[ind,];zij=zij||pxi[j,];end; lij=zij*bp; aiij=1/(1+exp(-lij)); scpij=ri[j-1]*(ri[j]-aiij)*zij; hpij=-ri[j-1]*aiij*(1-aiij)*zij`*zij; piij=piij*aiij; end; if j=1 then do;scpi=scpij;hpi=hpij;pii=piij;end; else do; scpi=scpi+scpij;hpi=hpi+hpij;pii=pii//piij;end; end; if (i=1) then do; scp=scpi; hp=hpi; pi=pii;end; else do; scp=scp//scpi; hp=hp+hpi; pi=pi//pii; end; esp=(scp[+,])`; end; finish; start ana; oldbe=be+1; do ite=1 to 20 while (max(abs(oldbe-be))>0.0001 ); oldbe=be; b=j(p,1,0); run gee; be=b; /*if wc='c' then orho=rho; if wc='ar' then orho= rho;*/ orho=rho; ophi=phi; /*if (wc='c'|wc='ar') then print es oldbe be phi rho ite; else print es oldbe be phi ite;*/ if ite=1 then do; if (max(abs(es))<0.00001) then do; hind=-h /* /phi*/; /*if d='po' then do; yy=y[loc(y>0)]; quasi=((r[loc(y>0)]#(yy#log(mu[loc(y>0)])-mu[loc(y>0)]-yy#(log(yy)-1)))[+])/phi; end; if d='bin' then quasi=((r#((y#log(mu/(1-mu)))+log(1-mu)))[+])/phi; if d='nor' then quasi=-0.5*(y-mu)`*(r#(y-mu))/phi;*/ end; else do; hind=.;end; end; end; if (max(abs(oldbe-be))<0.0001 & max(abs(es))<0.00001) then do; b=be; ht=-h; %if %length(&complete) > 0 %then %do; tsc=(sc*scp)*inv(/*scp`*scp*/-hp)*scp`;%end; %else tsc=j(p,n,0); sc=sc-tsc; se=inv(ht); se2=se*sc*sc`*se; se=sqrt(vecdiag(se2)); rrho=0; if (wc^='i') then rrho=rho; print Parameter oldbe be se ite rrho; if (p=pf/* & wc=&corr1*/) then mu0=mu; run cri; print ccri0 cri0; end; else do; ccri0=1E10;cri0=1E10; print "not converged"; end; finish; **********************************************************; y=&response; %let dist = %lowcase("&dist"); d=&dist; %if %length(&offset) = 0 %then %let offset = 0; off=&offset; if off=0 then off=j(nt,1,0); %if %length(&Mdy) = 0 %then %let Mdy = y0; %let Mdy = %lowcase(&Mdy); mpi=&Mdy; %if %length(&Mdx) = 0 %then %do;mpix= 0;%end; %else %do;mpix= 1;%end; %if %length(&Mdx) > 0 %then %do; px=&VarTest; %end; p=ncol(px); dimp=mpi+1+p; /* if (mpi='y0' & mpix=0) then dimp=1; if (mpi='y0' & mpix=1) then dimp=1+p; if (mpi='y1' & mpix=0) then dimp=2; if (mpi='y1' & mpix=1) then dimp=2+p; if (mpi='y2' & mpix=0) then dimp=3; if (mpi='y2' & mpix=1) then dimp=3+p; */ bp=j(dimp,1,0);bp[1]=0; oldbp=bp+1; %if %length(&complete) = 0 %then %do; pi=j(nt,1,1);r=j(nt,1,1);esp=0; %end; %else %do; run lik; %end; if (/*max(abs(oldbp-bp))<=0.0001 &*/ max(abs(esp))<0.00001) then do; %if %length(&complete) > 0 %then %do; sebp=sqrt(vecdiag(inv(-hp))); print "&complete = &M_ModelYP &M_ModelXP"; print MParameter bp sebp esp; %end; *************************************************************************************; **********Run model*********; %let nowmodel = 0; %do model_Lab = 1 %to &n_model; %do model_corr = 1 %to &n_corr; wc=&&corr&model_corr; x=&&Xmatric&model_Lab; nx= "&&XmLabel&model_Lab"; Parameter = "Intercept"; %let OutputPLab = %cmpres(&&XmLabel&model_Lab) _ ; %do %until(%verify(&&OutputPLab, &AString) < 1); %if %verify(&&OutputPLab, &AString) >= 1 %then %do; %let pos = %verify(&&OutputPLab, &AString); nowParameter = "%substr(&&OutputPLab,1,&pos-1)"; Parameter = Parameter//nowParameter; %let OutputPLab = %substr(&&OutputPLab,&pos+1); %end; %end; x=j(nrow(x),1,1)||x; p=ncol(x); pf=&Max_m_p; be=j(p,1,0); %let nowmodel = %eval(&nowmodel + 1); %if &model_Lab =1 and &model_corr =1 %then %do; print "--------------------------------------------------------------dropout model"; %end; %else %do; /* print "-----------------------------------------------------------------------------";*/ %end; print "Dist:" %upcase(&dist) ; print"wc:" (wc[1]); print %upcase("&response") " =" (nx[1]) ; run ana; print "---------------------------------------------------------------Candidate Model &nowmodel"; %if &model_Lab =1 and &model_corr =1 %then %do; /*if p=pf then do;*/MLICC=ccri0;MLIC=cri0;mo=nx; tw=wc;/*end; else do;tc=tc//cri0;mo=mo//nx;tw=tw//wc;end;*/ %end; %if &model_Lab >1 or &model_corr >1 %then %do; MLICC=MLICC//ccri0;MLIC=MLIC//cri0;mo=mo//nx;tw=tw//wc; %end; %end; %end; **********Run model*********; print "total # of models: " (nrow(tc)); print mo MLICC MLIC; sel=loc(MLICC=min(MLICC)); msel=mo[sel];tsel=MLICC[sel];wsel=tw[sel]; print "selected= " (msel[1]) "wc= " (wsel[1]); print "ccri0= " (tsel[1]); sel=loc(MLIC=min(MLIC)); msel=mo[sel];tsel=MLIC[sel];wsel=tw[sel]; print "selected= " (msel[1]) "wc= " (wsel[1]); print "cri0= " (tsel[1]); end; quit iml; %mend;