/* Functions that must be used with programs qpl3.txt and add3.txt. They must be included with the appropriate code in the files add3.txt and qpl3.txt. See files for specific line reading: %INCLUDE "..\IML\LIB2.IML"/ SOURCE2; Change the INCLUDE statement to point to where you have saved the lib2.iml code. */ /*******************************************************************************/ START qpl_will(x,kid,family,nkids,missing,nparms,lambda,covwill); sample=nrow(family); lambda=j(nparms,1,0); count=0; p=j(3,1,0); dl=j(nparms,1,0); addvars=j(nparms,nparms,0); addhs = j(nparms,nparms,0); v = j(nparms,nparms,0); h = j(nparms,nparms,0); hfill = j(nparms,nparms,0); dlbydlt = j(nparms,nparms,0); halfcov= j(nparms,nparms,0); covlambda= j(nparms,nparms,0); lambda=j(nparms,1,0); covwill=j(nparms,nparms,0); do kk=1 to sample; do ii=1 to nkids[kk]; if (missing[kk]=0) then count=count+1; end; end; likepre=0.0; loglike=0.0; do ll=1 to 25; addparms=j(nparms,1,0); addhs=j(nparms,nparms,0); addvars=j(nparms,nparms,0); v=j(nparms,nparms,0); covlambda=j(nparms,nparms,0); pluslambda=j(nparms,1,0); if (ll=1) then do; b0=lambda[1]; b2=lambda[2]; a000=lambda[3]; a002=lambda[4]; a010=lambda[5]; a012=lambda[6]; a020=lambda[7]; a022=lambda[8]; a110=lambda[9]; a112=lambda[10]; a120=lambda[11]; a122=lambda[12]; a220=lambda[13]; a222=lambda[14]; loglike=0; do kk=1 to sample; do ii=1 to nkids[kk]; if (missing[kk]=0) then do; num0= exp(b0*x[kk,ii]+a000*(family[kk]=1111)+a020*(family[kk]=1122)+a220*(family[kk]=2222)+ a010*(family[kk]=1112)+a110*(family[kk]=1212)+a120*(family[kk]=1222)); num2=exp(b2*x[kk,ii]+a002*(family[kk]=1111)+a022*(family[kk]=1122)+a222*(family[kk]=2222)+ a012*(family[kk]=1112)+a112*(family[kk]=1212)+a122*(family[kk]=1222)); denom=(1+num0+num2); p[1]=num0/denom; p[2]=1/denom; p[3]=num2/denom; loglike = loglike + log(((kid[kk,ii]=0)*p[1])+((kid[kk,ii]=1)*p[2])+((kid[kk,ii]=2)*p[3])); end; end; end; end; lastone=0; do kk=1 to sample; do ii=1 to nkids[kk]; if (missing[kk]=0) then do; lastone=lastone+1; num0=exp(b0*x[kk,ii]+a000*(family[kk]=1111)+a020*(family[kk]=1122)+ a220*(family[kk]=2222)+a010*(family[kk]=1112)+ a110*(family[kk]=1212)+a120*(family[kk]=1222)); num2=exp(b2*x[kk,ii]+a002*(family[kk]=1111)+a022*(family[kk]=1122)+ a222*(family[kk]=2222)+a012*(family[kk]=1112)+ a112*(family[kk]=1212)+a122*(family[kk]=1222)); denom=(1+num0+num2); p[1]=num0/denom; p[2]=1/denom; p[3]=num2/denom; dl[1]=(kid[kk,ii]=0)*x[kk,ii]-x[kk,ii]*p[1]; dl[2]=(kid[kk,ii]=2)*x[kk,ii]-x[kk,ii]*p[3]; dl[3]=((kid[kk,ii]=0)-p[1])*(family[kk]=1111); dl[4]=((kid[kk,ii]=2)-p[3])*(family[kk]=1111); dl[5]=((kid[kk,ii]=0)-p[1])*(family[kk]=1112); dl[6]=((kid[kk,ii]=2)-p[3])*(family[kk]=1112); dl[7]=((kid[kk,ii]=0)-p[1])*(family[kk]=1122); dl[8]=((kid[kk,ii]=2)-p[3])*(family[kk]=1122); dl[9]=((kid[kk,ii]=0)-p[1])*(family[kk]=1212); dl[10]=((kid[kk,ii]=2)-p[3])*(family[kk]=1212); dl[11]=((kid[kk,ii]=0)-p[1])*(family[kk]=1222); dl[12]=((kid[kk,ii]=2)-p[3])*(family[kk]=1222); dl[13]=((kid[kk,ii]=0)-p[1])*(family[kk]=2222); dl[14]=((kid[kk,ii]=2)-p[3])*(family[kk]=2222); h[1,1]=(x[kk,ii]*x[kk,ii])*p[1]*(p[1]-1); h[1,2]=(x[kk,ii]*x[kk,ii])*p[1]*p[3]; h[1,3]=(x[kk,ii])*p[1]*(p[1]-1)*(family[kk]=1111); h[1,4]=(x[kk,ii])*p[1]*p[3]*(family[kk]=1111); h[1,5]=(x[kk,ii])*p[1]*(p[1]-1)*(family[kk]=1112); h[1,6]=(x[kk,ii])*p[1]*p[3]*(family[kk]=1112); h[1,7]=(x[kk,ii])*p[1]*(p[1]-1)*(family[kk]=1122); h[1,8]=(x[kk,ii])*p[1]*p[3]*(family[kk]=1122); h[1,9]=(x[kk,ii])*p[1]*(p[1]-1)*(family[kk]=1212); h[1,10]=(x[kk,ii])*p[1]*p[3]*(family[kk]=1212); h[1,11]=(x[kk,ii])*p[1]*(p[1]-1)*(family[kk]=1222); h[1,12]=(x[kk,ii])*p[1]*p[3]*(family[kk]=1222); h[1,13]=(x[kk,ii])*p[1]*(p[1]-1)*(family[kk]=2222); h[1,14]=(x[kk,ii])*p[1]*p[3]*(family[kk]=2222); h[2,2]=(x[kk,ii]*x[kk,ii])*p[3]*(p[3]-1); h[2,3]=(x[kk,ii])*p[1]*p[3]*(family[kk]=1111); h[2,4]=(x[kk,ii])*p[3]*(p[3]-1)*(family[kk]=1111); h[2,5]=(x[kk,ii])*p[1]*p[3]*(family[kk]=1112); h[2,6]=(x[kk,ii])*p[3]*(p[3]-1)*(family[kk]=1112); h[2,7]=(x[kk,ii])*p[1]*p[3]*(family[kk]=1122); h[2,8]=(x[kk,ii])*p[3]*(p[3]-1)*(family[kk]=1122); h[2,9]=(x[kk,ii])*p[1]*p[3]*(family[kk]=1212); h[2,10]=(x[kk,ii])*p[3]*(p[3]-1)*(family[kk]=1212); h[2,11]=(x[kk,ii])*p[1]*p[3]*(family[kk]=1222); h[2,12]=(x[kk,ii])*p[3]*(p[3]-1)*(family[kk]=1222); h[2,13]=(x[kk,ii])*p[1]*p[3]*(family[kk]=2222); h[2,14]=(x[kk,ii])*p[3]*(p[3]-1)*(family[kk]=2222); h[3,3]=p[1]*(p[1]-1)*(family[kk]=1111); h[3,4]=p[1]*p[3]*(family[kk]=1111); h[3,5]=0; h[3,6]=0; h[3,7]=0; h[3,8]=0; h[3,9]=0; h[3,10]=0; h[3,11]=0; h[3,12]=0; h[3,13]=0; h[3,14]=0; h[4,4]=p[3]*(p[3]-1)*(family[kk]=1111); h[4,5]=0; h[4,6]=0; h[4,7]=0; h[4,8]=0; h[4,9]=0; h[4,10]=0; h[4,11]=0; h[4,12]=0; h[4,13]=0; h[4,14]=0; h[5,5]=p[1]*(p[1]-1)*(family[kk]=1112); h[5,6]=p[1]*p[3]*(family[kk]=1112); h[5,7]=0; h[5,8]=0; h[5,9]=0; h[5,10]=0; h[5,11]=0; h[5,12]=0; h[5,13]=0; h[5,14]=0; h[6,6]=p[3]*(p[3]-1)*(family[kk]=1112); h[6,7]=0; h[6,8]=0; h[6,9]=0; h[6,10]=0; h[6,11]=0; h[6,12]=0; h[6,13]=0; h[6,14]=0; h[7,7]=p[1]*(p[1]-1)*(family[kk]=1122); h[7,8]=p[1]*p[3]*(family[kk]=1122); h[7,9]=0; h[7,10]=0; h[7,11]=0; h[7,12]=0; h[7,13]=0; h[7,14]=0; h[8,8]=p[3]*(p[3]-1)*(family[kk]=1122); h[8,9]=0; h[8,10]=0; h[8,11]=0; h[8,12]=0; h[8,13]=0; h[8,14]=0; h[9,9]=p[1]*(p[1]-1)*(family[kk]=1212); h[9,10]=p[1]*p[3]*(family[kk]=1212); h[9,11]=0; h[9,12]=0; h[9,13]=0; h[9,14]=0; h[10,10]=p[3]*(p[3]-1)*(family[kk]=1212); h[10,11]=0; h[10,12]=0; h[10,13]=0; h[10,14]=0; h[11,11]=p[1]*(p[1]-1)*(family[kk]=1222); h[11,12]=p[1]*p[3]*(family[kk]=1222); h[11,13]=0; h[11,14]=0; h[12,12]=p[3]*(p[3]-1)*(family[kk]=1222); h[12,13]=0; h[12,14]=0; h[13,13]=p[1]*(p[1]-1)*(family[kk]=2222); h[13,14]=p[1]*p[3]*(family[kk]=2222); h[14,14]=p[3]*(p[3]-1)*(family[kk]=2222); do i=1 to nparms; do j=1 to nparms; if (j=i) then hfill[j,i]=h[i,j]; end; end; do i=1 to nparms; dl[i]=(1/nkids[kk])*dl[i]; do j=1 to nparms; hfill[i,j]=(1/nkids[kk])*hfill[i,j]; end; end; addparms=addparms+dl; addhs=addhs-hfill; if (ii=1) then do; v=v+addvars; addvars=dl*dl`; end; if (ii^=1) then do; dlbydlt=dl*dl`; addvars=addvars+dlbydlt; end; if (lastone=count) then do; v=v+addvars; covlambda=inv(addhs); pluslambda=covlambda*addparms; lambda=lambda+pluslambda; halfcov=covlambda*v; covwill=halfcov*covlambda; flagit=b0-lambda[1]+ b2-lambda[2]+ a000-lambda[3]+ a002-lambda[4]+ a010-lambda[5]+ a012-lambda[6]+ a020-lambda[7]+ a022-lambda[8]+ a110-lambda[9]+ a112-lambda[10]+ a120-lambda[11]+ a122-lambda[12]+ a220-lambda[13]+ a222-lambda[14]; end; end; end; end; likepre=loglike; loglike=0; oldb0=b0; oldb2=b2; olda000=a000; olda002=a002; olda010=a010; olda012=a012; olda020=a020; olda022=a022; olda110=a110; olda112=a112; olda120=a120; olda122=a122; olda220=a220; olda222=a222; b0=lambda[1]; b2=lambda[2]; a000=lambda[3]; a002=lambda[4]; a010=lambda[5]; a012=lambda[6]; a020=lambda[7]; a022=lambda[8]; a110=lambda[9]; a112=lambda[10]; a120=lambda[11]; a122=lambda[12]; a220=lambda[13]; a222=lambda[14]; do jj=1 to 15; if(jj=1) then do; loglike=0; do kk=1 to sample; do ii=1 to nkids[kk]; if (missing[kk]=0) then do; num0=exp(b0*x[kk,ii]+a000*(family[kk]=1111)+a020*(family[kk]=1122)+ a220*(family[kk]=2222)+a010*(family[kk]=1112)+ a110*(family[kk]=1212)+a120*(family[kk]=1222)); num2=exp(b2*x[kk,ii]+a002*(family[kk]=1111)+a022*(family[kk]=1122)+ a222*(family[kk]=2222)+a012*(family[kk]=1112)+ a112*(family[kk]=1212)+a122*(family[kk]=1222)); denom=(1+num0+num2); p[1]=num0/denom; p[2]=1/denom; p[3]=num2/denom; loglike=loglike + log(((kid[kk,ii]=0)*p[1])+((kid[kk,ii]=1)*p[2])+((kid[kk,ii]=2)*p[3])); end; end; end; end; if (loglike=i) then hfill[j,i]=h[i,j]; end; end; do i=1 to nparms2; dl[i]=(1/nkids[kk])*dl[i]; do j=1 to nparms2; hfill[i,j]=(1/nkids[kk])*hfill[i,j]; end; end; addparms=addparms+dl; addhs=addhs-hfill; if (ii=1) then do; v=v+addvars; addvars=dl*dl`; end; if (ii^=1) then do; dlbydlt=dl*dl`; addvars=addvars+dlbydlt; end; if (lastone=count) then do; v=v+addvars; covlambda=inv(addhs); pluslambda=covlambda*addparms; mlambda=mlambda+pluslambda; halfcov=covlambda*v; covwill=halfcov*covlambda; flagit=c01-mlambda[1]+ c02-mlambda[2]+ c11-mlambda[3]+ c12-mlambda[4]+ c22-mlambda[5]+ b01-mlambda[6]+ b02-mlambda[7]+ b11-mlambda[8]+ b12-mlambda[9]+ b22-mlambda[10]+ a01-mlambda[11]+ a02-mlambda[12]+ a11-mlambda[13]+ a12-mlambda[14]+ a22-mlambda[15]; end; end; end; end; likepre=loglike; loglike=0; oldc01=c01; oldc02=c02; oldc11=c11; oldc12=c12; oldc22=c22; oldb01=b01; oldb02=b02; oldb11=b11; oldb12=b12; oldb22=b22; olda01=a01; olda02=a02; olda11=a11; olda12=a12; olda22=a22; c01=mlambda[1]; c02=mlambda[2]; c11=mlambda[3]; c12=mlambda[4]; c22=mlambda[5]; b01=mlambda[6]; b02=mlambda[7]; b11=mlambda[8]; b12=mlambda[9]; b22=mlambda[10]; a01=mlambda[11]; a02=mlambda[12]; a11=mlambda[13]; a12=mlambda[14]; a22=mlambda[15]; do jj=1 to 15; if(jj=1) then do; loglike=0; do kk=1 to sample; do ii=1 to nkids[kk]; if (missing[kk]=0) then do; num0=exp(c01*x[kk,ii]*x[kk,ii]+b01*x[kk,ii]+a01); num1=exp(c02*x[kk,ii]*x[kk,ii]+b02*x[kk,ii]+a02); num2=exp(c11*x[kk,ii]*x[kk,ii]+b11*x[kk,ii]+a11); num3=exp(c12*x[kk,ii]*x[kk,ii]+b12*x[kk,ii]+a12); num4=exp(c22*x[kk,ii]*x[kk,ii]+b22*x[kk,ii]+a22); denom=(num0+num1+num2+num3+num4+1); p[1]=1/denom; p[2]=num0/denom; p[3]=num1/denom; p[4]=num2/denom; p[5]=num3/denom; p[6]=num4/denom; loglike=loglike+log(((family[kk]=1111)*p[1])+((family[kk]=1112)*p[2])+((family[kk]=1122)*p[3]) +((family[kk]=1212)*p[4])+((family[kk]=1222)*p[5])+((family[kk]=2222)*p[6])); end; end; end; end; if (loglike=rp1) then missparent[kk]=1; notmissing[kk]=0; end; if ((obsparent[kk]=0) & (kid[kk,1]=1)) then do; r=RanUni(0); rp1=.5*p1_01*p01/(.5*p1_01*p01+.5*p1_02*p02); if (r=rp1) then missparent[kk]=2; notmissing[kk]=0; end; if ((obsparent[kk]=1) & (kid[kk,1]=0)) then do; r=RanUni(0); rp1=.5*p0_01*p01/(.5*p0_01*p01+p0_11*p11); if (r=rp1) then missparent[kk]=1; notmissing[kk]=0; end; if ((obsparent[kk]=1) & (kid[kk,1]=1)) then do; r=RanUni(0); rp1=.5*p1_01*p01/(.5*p1_01*p01+p1_11*p11+.5*p1_12*p12); rp2=(.5*p1_01*p01+p1_11*p11)/(.5*p1_01*p01+p1_11*p11+.5*p1_12*p12); if (r=rp1) & (r=rp2) then missparent[kk]=2; notmissing[kk]=0; end; if ((obsparent[kk]=1) & (kid[kk,1]=2)) then do; r=RanUni(0); rp1=p2_11*p11/(p2_11*p11+.5*p2_12*p12); if (r=rp1) then missparent[kk]=2; notmissing[kk]=0; end; if ((obsparent[kk]=2) & (kid[kk,1]=1)) then do; r=RanUni(0); rp1=.5*p1_02*p02/(.5*p1_02*p02+.5*p1_12*p12); if (r=rp1) then missparent[kk]=1; notmissing[kk]=0; end; if ((obsparent[kk]=2) & (kid[kk,1]=2)) then do; r=RanUni(0); rp1=.5*p2_12*p12/(.5*p2_12*p12+p2_22*p22); if (r=rp1) then missparent[kk]=2; notmissing[kk]=0; end; end; if (nkids[kk]=2) then do; if ((obsparent[kk]=0) & (kid[kk,1]=0) & (kid[kk,2]=0)) then do; r=RanUni(0); rp0=0; do jj=1 to nkids[kk]; p0_00=exp(b0*x[kk,jj]+c000)/(exp(b0*x[kk,jj]+c000)+1+exp(b2*x[kk,jj]+c002)); p0_01=exp(b0*x[kk,jj]+c010)/(exp(b0*x[kk,jj]+c010)+1+exp(b2*x[kk,jj]+c012)); total=exp(q00*x[kk,jj]*x[kk,jj]+a00*x[kk,jj]+u00)+ exp(q01*x[kk,jj]*x[kk,jj]+a01*x[kk,jj]+u01)+ exp(q11*x[kk,jj]*x[kk,jj]+a11*x[kk,jj]+u11)+ exp(q12*x[kk,jj]*x[kk,jj]+a12*x[kk,jj]+u12)+ exp(q02*x[kk,jj]*x[kk,jj]+a02*x[kk,jj]+u02)+ exp(q22*x[kk,jj]*x[kk,jj]+a22*x[kk,jj]+u22); p00=exp(q00*x[kk,jj]*x[kk,jj]+a00*x[kk,jj]+u00)/total; p01=exp(q01*x[kk,jj]*x[kk,jj]+a01*x[kk,jj]+u01)/total; rp0=rp0+p0_00*p00/(p0_00*p00+.5*p0_01*p01); end; rp=rp0/2; if (r=rp) then missparent[kk]=1; notmissing[kk]=0; end; if ((obsparent[kk]=1) & (kid[kk,1]=0) & (kid[kk,2]=0)) then do; r=RanUni(0); rp0=0; do jj=1 to nkids[kk]; total=exp(q00*x[kk,jj]*x[kk,jj]+a00*x[kk,jj]+u00)+ exp(q01*x[kk,jj]*x[kk,jj]+a01*x[kk,jj]+u01)+ exp(q11*x[kk,jj]*x[kk,jj]+a11*x[kk,jj]+u11)+ exp(q12*x[kk,jj]*x[kk,jj]+a12*x[kk,jj]+u12)+ exp(q02*x[kk,jj]*x[kk,jj]+a02*x[kk,jj]+u02)+ exp(q22*x[kk,jj]*x[kk,jj]+a22*x[kk,jj]+u22); p01=exp(q01*x[kk,jj]*x[kk,jj]+a01*x[kk,jj]+u01)/total; p11=exp(q11*x[kk,jj]*x[kk,jj]+a11*x[kk,jj]+u11)/total; p0_01=exp(b0*x[kk,jj]+c010)/(exp(b0*x[kk,jj]+c010)+1+exp(b2*x[kk,jj]+c012)); p0_11=exp(b0*x[kk,jj]+c110)/(exp(b0*x[kk,jj]+c110)+1+exp(b2*x[kk,jj]+c112)); rp0=rp0 + .5*p0_01*p01/(.5*p0_01*p01+p0_11*p11); end; rp=rp0/2; if (r=rp) then missparent[kk]=1; notmissing[kk]=0; end; if (((obsparent[kk]=0) & (kid[kk,1]=0) & (kid[kk,2]=1)) | ((obsparent[kk]=0) & (kid[kk,1]=1) & (kid[kk,2]=0))) then do; missparent[kk]=1; notmissing[kk]=0; end; if (((obsparent[kk]=1) & (kid[kk,1]=0) & (kid[kk,2]=1)) | ((obsparent[kk]=1) & (kid[kk,1]=1) & (kid[kk,2]=0))) then do; r=RanUni(0); rp0=0; do jj=1 to nkids[kk]; if (kid[kk,jj]=0) then do; total=exp(q00*x[kk,jj]*x[kk,jj]+a00*x[kk,jj]+u00)+ exp(q01*x[kk,jj]*x[kk,jj]+a01*x[kk,jj]+u01)+ exp(q11*x[kk,jj]*x[kk,jj]+a11*x[kk,jj]+u11)+ exp(q12*x[kk,jj]*x[kk,jj]+a12*x[kk,jj]+u12)+ exp(q02*x[kk,jj]*x[kk,jj]+a02*x[kk,jj]+u02)+ exp(q22*x[kk,jj]*x[kk,jj]+a22*x[kk,jj]+u22); p0_01=exp(b0*x[kk,jj]+c010)/(exp(b0*x[kk,jj]+c010)+1+exp(b2*x[kk,jj]+c012)); p01=exp(q01*x[kk,jj]*x[kk,jj]+a01*x[kk,jj]+u01)/total; p0_11=exp(b0*x[kk,jj]+c110)/(exp(b0*x[kk,jj]+c110)+1+exp(b2*x[kk,jj]+c112)); p11=exp(q11*x[kk,jj]*x[kk,jj]+a11*x[kk,jj]+u11)/total; rp0=rp0+.5*p0_01*p01/(.5*p0_01*p01+p0_11*p11); end; if (kid[kk,jj]=1) then do; total=exp(q00*x[kk,jj]*x[kk,jj]+a00*x[kk,jj]+u00)+ exp(q01*x[kk,jj]*x[kk,jj]+a01*x[kk,jj]+u01)+ exp(q11*x[kk,jj]*x[kk,jj]+a11*x[kk,jj]+u11)+ exp(q12*x[kk,jj]*x[kk,jj]+a12*x[kk,jj]+u12)+ exp(q02*x[kk,jj]*x[kk,jj]+a02*x[kk,jj]+u02)+ exp(q22*x[kk,jj]*x[kk,jj]+a22*x[kk,jj]+u22); p1_01=1/(exp(b0*x[kk,jj]+c010)+1+exp(b2*x[kk,jj]+c012)); p01=exp(q01*x[kk,jj]*x[kk,jj]+a01*x[kk,jj]+u01)/total; p1_11=1/(exp(b0*x[kk,jj]+c110)+1+exp(b2*x[kk,jj]+c112)); p11=exp(q11*x[kk,jj]*x[kk,jj]+a11*x[kk,jj]+u11)/total; rp0=rp0+.5*p1_01*p01/(.5*p1_01*p01+p1_11*p11); end; end; rp=rp0/2; if (r=rp) then missparent[kk]=1; notmissing[kk]=0; end; if (((obsparent[kk]=1) & (kid[kk,1]=0) & (kid[kk,2]=2)) | ((obsparent[kk]=1) & (kid[kk,1]=2) & (kid[kk,2]=0))) then do; missparent[kk]=1; notmissing[kk]=0; end; if ((obsparent[kk]=0) & (kid[kk,1]=1) & (kid[kk,2]=1)) then do; r=RanUni(0); rp0=0; do jj=1 to nkids[kk]; p1_02=1/(exp(b0*x[kk,jj]+c020)+1+exp(b2*x[kk,jj]+c022)); p1_01=1/(exp(b0*x[kk,jj]+c010)+1+exp(b2*x[kk,jj]+c012)); total=exp(q00*x[kk,jj]*x[kk,jj]+a00*x[kk,jj]+u00)+ exp(q01*x[kk,jj]*x[kk,jj]+a01*x[kk,jj]+u01)+ exp(q11*x[kk,jj]*x[kk,jj]+a11*x[kk,jj]+u11)+ exp(q12*x[kk,jj]*x[kk,jj]+a12*x[kk,jj]+u12)+ exp(q02*x[kk,jj]*x[kk,jj]+a02*x[kk,jj]+u02)+ exp(q22*x[kk,jj]*x[kk,jj]+a22*x[kk,jj]+u22); p01=exp(q01*x[kk,jj]*x[kk,jj]+a01*x[kk,jj]+u01)/total; p02=exp(q02*x[kk,jj]*x[kk,jj]+a02*x[kk,jj]+u02)/total; rp0=rp0+.5*p1_01*p01/(.5*p1_01*p01+.5*p1_02*p02); end; rp=rp0/2; if (r=rp) then missparent[kk]=2; notmissing[kk]=0; end; if ((obsparent[kk]=1) & (kid[kk,1]=1) & (kid[kk,2]=1)) then do; r=RanUni(0); rp1=0; rp2=0; do jj=1 to nkids[kk]; p1_01=1/(exp(b0*x[kk,jj]+c010)+1+exp(b2*x[kk,jj]+c012)); p1_11=1/(exp(b0*x[kk,jj]+c110)+1+exp(b2*x[kk,jj]+c112)); p1_12=1/(exp(b0*x[kk,jj]+c120)+1+exp(b2*x[kk,jj]+c122)); total=exp(q00*x[kk,jj]*x[kk,jj]+a00*x[kk,jj]+u00)+ exp(q01*x[kk,jj]*x[kk,jj]+a01*x[kk,jj]+u01)+ exp(q11*x[kk,jj]*x[kk,jj]+a11*x[kk,jj]+u11)+ exp(q12*x[kk,jj]*x[kk,jj]+a12*x[kk,jj]+u12)+ exp(q02*x[kk,jj]*x[kk,jj]+a02*x[kk,jj]+u02)+ exp(q22*x[kk,jj]*x[kk,jj]+a22*x[kk,jj]+u22); p01=exp(q01*x[kk,jj]*x[kk,jj]+a01*x[kk,jj]+u01)/total; p11=exp(q11*x[kk,jj]*x[kk,jj]+a11*x[kk,jj]+u11)/total; p12=exp(q12*x[kk,jj]*x[kk,jj]+a12*x[kk,jj]+u12)/total; rp1= rp1 + .5*p1_01*p01/(.5*p1_01*p01+p1_11*p11+.5*p1_12*p12); rp2= rp2 + (.5*p1_01*p01+p1_11*p11)/(.5*p1_01*p01+p1_11*p11+.5*p1_12*p12); end; rp1=rp1/2; rp2=rp2/2; if (r=rp1) & (r=rp2) then missparent[kk]=2; notmissing[kk]=0; end; if ((obsparent[kk]=2) & (kid[kk,1]=1) & (kid[kk,2]=1)) then do; r=RanUni(0); rp0=0; do jj=1 to nkids[kk]; p1_02=1/(exp(b0*x[kk,jj]+c020)+1+exp(b2*x[kk,jj]+c022)); p1_12=1/(exp(b0*x[kk,jj]+c120)+1+exp(b2*x[kk,jj]+c122)); total=exp(q00*x[kk,jj]*x[kk,jj]+a00*x[kk,jj]+u00)+ exp(q01*x[kk,jj]*x[kk,jj]+a01*x[kk,jj]+u01)+ exp(q11*x[kk,jj]*x[kk,jj]+a11*x[kk,jj]+u11)+ exp(q12*x[kk,jj]*x[kk,jj]+a12*x[kk,jj]+u12)+ exp(q02*x[kk,jj]*x[kk,jj]+a02*x[kk,jj]+u02)+ exp(q22*x[kk,jj]*x[kk,jj]+a22*x[kk,jj]+u22); p12=exp(q12*x[kk,jj]*x[kk,jj]+a12*x[kk,jj]+u12)/total; p02=exp(q02*x[kk,jj]*x[kk,jj]+a02*x[kk,jj]+u02)/total; rp0= rp0 + .5*p1_02*p02/(.5*p1_02*p02+.5*p1_12*p12); end; rp=rp0/2; if (r=rp) then missparent[kk]=1; notmissing[kk]=0; end; if (((obsparent[kk]=1) & (kid[kk,1]=1) & (kid[kk,2]=2)) | ((obsparent[kk]=1) & (kid[kk,1]=2) & (kid[kk,2]=1))) then do; r=RanUni(0); rp0=0; do jj=1 to nkids[kk]; if (kid[kk,jj]=1) then do; p1_11=1/(exp(b0*x[kk,jj]+c110)+1+exp(b2*x[kk,jj]+c112)); p1_12=1/(exp(b0*x[kk,jj]+c120)+1+exp(b2*x[kk,jj]+c122)); total=exp(q00*x[kk,jj]*x[kk,jj]+a00*x[kk,jj]+u00)+ exp(q01*x[kk,jj]*x[kk,jj]+a01*x[kk,jj]+u01)+ exp(q11*x[kk,jj]*x[kk,jj]+a11*x[kk,jj]+u11)+ exp(q12*x[kk,jj]*x[kk,jj]+a12*x[kk,jj]+u12)+ exp(q02*x[kk,jj]*x[kk,jj]+a02*x[kk,jj]+u02)+ exp(q22*x[kk,jj]*x[kk,jj]+a22*x[kk,jj]+u22); p11=exp(q11*x[kk,jj]*x[kk,jj]+a11*x[kk,jj]+u11)/total; p12=exp(q12*x[kk,jj]*x[kk,jj]+a12*x[kk,jj]+u12)/total; rp0=rp0+(p1_11*p11)/(p1_11*p11+.5*p1_12*p12); end; if (kid[kk,jj]=1) then do; p2_11=exp(b2*x[kk,jj]+c112)/(exp(b0*x[kk,jj]+c110)+1+exp(b2*x[kk,jj]+c112)); p2_12=exp(b2*x[kk,jj]+c122)/(exp(b0*x[kk,jj]+c120)+1+exp(b2*x[kk,jj]+c122)); total=exp(q00*x[kk,jj]*x[kk,jj]+a00*x[kk,jj]+u00)+ exp(q01*x[kk,jj]*x[kk,jj]+a01*x[kk,jj]+u01)+ exp(q11*x[kk,jj]*x[kk,jj]+a11*x[kk,jj]+u11)+ exp(q12*x[kk,jj]*x[kk,jj]+a12*x[kk,jj]+u12)+ exp(q02*x[kk,jj]*x[kk,jj]+a02*x[kk,jj]+u02)+ exp(q22*x[kk,jj]*x[kk,jj]+a22*x[kk,jj]+u22); p11=exp(q11*x[kk,jj]*x[kk,jj]+a11*x[kk,jj]+u11)/total; p12=exp(q12*x[kk,jj]*x[kk,jj]+a12*x[kk,jj]+u12)/total; rp0=rp0+p2_11*p11/(p2_11*p11+.5*p2_12*p12); end; end; rp=rp0/2; if (r=rp) then missparent[kk]=2; notmissing[kk]=0; end; if (((obsparent[kk]=2) & (kid[kk,1]=1) & (kid[kk,2]=2)) | ((obsparent[kk]=2) & (kid[kk,1]=2) & (kid[kk,2]=1))) then do; missparent[kk]=1; notmissing[kk]=0; end; if ((obsparent[kk]=1) & (kid[kk,1]=2) & (kid[kk,2]=2)) then do; r=RanUni(0); rp0=0; do jj=1 to nkids[kk]; p2_11=exp(b2*x[kk,jj]+c112)/(exp(b0*x[kk,jj]+c110)+1+exp(b2*x[kk,jj]+c112)); p2_12=exp(b2*x[kk,jj]+c122)/(exp(b0*x[kk,jj]+c120)+1+exp(b2*x[kk,jj]+c122)); total=exp(q00*x[kk,jj]*x[kk,jj]+a00*x[kk,jj]+u00)+ exp(q01*x[kk,jj]*x[kk,jj]+a01*x[kk,jj]+u01)+ exp(q11*x[kk,jj]*x[kk,jj]+a11*x[kk,jj]+u11)+ exp(q12*x[kk,jj]*x[kk,jj]+a12*x[kk,jj]+u12)+ exp(q02*x[kk,jj]*x[kk,jj]+a02*x[kk,jj]+u02)+ exp(q22*x[kk,jj]*x[kk,jj]+a22*x[kk,jj]+u22); p11=exp(q11*x[kk,jj]*x[kk,jj]+a11*x[kk,jj]+u11)/total; p12=exp(q12*x[kk,jj]*x[kk,jj]+a12*x[kk,jj]+u12)/total; rp0=rp0+p2_11*p11/(p2_11*p11+.5*p2_12*p12); end; rp=rp0/2; if (r=rp) then missparent[kk]=2; notmissing[kk]=0; end; if ((obsparent[kk]=2) & (kid[kk,1]=2) & (kid[kk,2]=2)) then do; r=RanUni(0); rp0=0; do jj=1 to nkids[kk]; p2_12=exp(b2*x[kk,jj]+c122)/(exp(b0*x[kk,jj]+c120)+1+exp(b2*x[kk,jj]+c122)); p2_22=exp(b2*x[kk,jj]+c222)/(exp(b0*x[kk,jj]+c220)+1+exp(b2*x[kk,jj]+c222)); total=exp(q00*x[kk,jj]*x[kk,jj]+a00*x[kk,jj]+u00)+ exp(q01*x[kk,jj]*x[kk,jj]+a01*x[kk,jj]+u01)+ exp(q11*x[kk,jj]*x[kk,jj]+a11*x[kk,jj]+u11)+ exp(q12*x[kk,jj]*x[kk,jj]+a12*x[kk,jj]+u12)+ exp(q02*x[kk,jj]*x[kk,jj]+a02*x[kk,jj]+u02)+ exp(q22*x[kk,jj]*x[kk,jj]+a22*x[kk,jj]+u22); p12=exp(q12*x[kk,jj]*x[kk,jj]+a12*x[kk,jj]+u12)/total; p22=exp(q22*x[kk,jj]*x[kk,jj]+a22*x[kk,jj]+u22)/total; rp0=rp0+.5*p2_12*p12/(.5*p2_12*p12+p2_22*p22); end; rp=rp0/2; if (r=rp) then missparent[kk]=2; notmissing[kk]=0; end; end; if (nkids[kk]=3) then do; if ((obsparent[kk]=0) & (kid[kk,1]=0) & (kid[kk,2]=0) & (kid[kk,3]=0)) then do; r=RanUni(0); rp0=0; do jj=1 to nkids[kk]; p0_00=exp(b0*x[kk,jj]+c000)/(exp(b0*x[kk,jj]+c000)+1+exp(b2*x[kk,jj]+c002)); p0_01=exp(b0*x[kk,jj]+c010)/(exp(b0*x[kk,jj]+c010)+1+exp(b2*x[kk,jj]+c012)); total=exp(q00*x[kk,jj]*x[kk,jj]+a00*x[kk,jj]+u00)+ exp(q01*x[kk,jj]*x[kk,jj]+a01*x[kk,jj]+u01)+ exp(q11*x[kk,jj]*x[kk,jj]+a11*x[kk,jj]+u11)+ exp(q12*x[kk,jj]*x[kk,jj]+a12*x[kk,jj]+u12)+ exp(q02*x[kk,jj]*x[kk,jj]+a02*x[kk,jj]+u02)+ exp(q22*x[kk,jj]*x[kk,jj]+a22*x[kk,jj]+u22); p00=exp(q00*x[kk,jj]*x[kk,jj]+a00*x[kk,jj]+u00)/total; p01=exp(q01*x[kk,jj]*x[kk,jj]+a01*x[kk,jj]+u01)/total; rp0=rp0+p0_00*p00/(p0_00*p00+.5*p0_01*p01); end; rp=rp0/3; if (r=rp) then missparent[kk]=1; notmissing[kk]=0; end; if ((obsparent[kk]=1) & (kid[kk,1]=0) & (kid[kk,2]=0) & (kid[kk,3]=0)) then do; r=RanUni(0); rp0=0; do jj=1 to nkids[kk]; total=exp(q00*x[kk,jj]*x[kk,jj]+a00*x[kk,jj]+u00)+ exp(q01*x[kk,jj]*x[kk,jj]+a01*x[kk,jj]+u01)+ exp(q11*x[kk,jj]*x[kk,jj]+a11*x[kk,jj]+u11)+ exp(q12*x[kk,jj]*x[kk,jj]+a12*x[kk,jj]+u12)+ exp(q02*x[kk,jj]*x[kk,jj]+a02*x[kk,jj]+u02)+ exp(q22*x[kk,jj]*x[kk,jj]+a22*x[kk,jj]+u22); p0_01=exp(b0*x[kk,jj]+c010)/(exp(b0*x[kk,jj]+c010)+1+exp(b2*x[kk,jj]+c012)); p0_11=exp(b0*x[kk,jj]+c110)/(exp(b0*x[kk,jj]+c110)+1+exp(b2*x[kk,jj]+c112)); p01=exp(q01*x[kk,jj]*x[kk,jj]+a01*x[kk,jj]+u01)/total; p11=exp(q11*x[kk,jj]*x[kk,jj]+a11*x[kk,jj]+u11)/total; rp0= rp0 + .5*p0_01*p01/(.5*p0_01*p01+p0_11*p11); end; rp=rp0/3; if (r=rp) then missparent[kk]=1; notmissing[kk]=0; end; if (((obsparent[kk]=0) & (kid[kk,1]=0) & (kid[kk,2]=0) & (kid[kk,3]=1)) | ((obsparent[kk]=0) & (kid[kk,1]=0) & (kid[kk,2]=1) & (kid[kk,3]=0)) | ((obsparent[kk]=0) & (kid[kk,1]=1) & (kid[kk,2]=0) & (kid[kk,3]=0)) | ((obsparent[kk]=0) & (kid[kk,1]=0) & (kid[kk,2]=1) & (kid[kk,3]=1)) | ((obsparent[kk]=0) & (kid[kk,1]=1) & (kid[kk,2]=0) & (kid[kk,3]=1)) | ((obsparent[kk]=0) & (kid[kk,1]=1) & (kid[kk,2]=1) & (kid[kk,3]=0))) then do; missparent[kk]=1; notmissing[kk]=0; end; if (((obsparent[kk]=1) & (kid[kk,1]=0) & (kid[kk,2]=0) & (kid[kk,3]=1)) | ((obsparent[kk]=1) & (kid[kk,1]=0) & (kid[kk,2]=1) & (kid[kk,3]=0)) | ((obsparent[kk]=1) & (kid[kk,1]=1) & (kid[kk,2]=0) & (kid[kk,3]=0)) | ((obsparent[kk]=1) & (kid[kk,1]=0) & (kid[kk,2]=1) & (kid[kk,3]=1)) | ((obsparent[kk]=1) & (kid[kk,1]=1) & (kid[kk,2]=0) & (kid[kk,3]=1)) | ((obsparent[kk]=1) & (kid[kk,1]=1) & (kid[kk,2]=1) & (kid[kk,3]=0))) then do; r=RanUni(0); rp0=0; do jj=1 to nkids[kk]; if (kid[kk,jj]=0) then do; total=exp(q00*x[kk,jj]*x[kk,jj]+a00*x[kk,jj]+u00)+ exp(q01*x[kk,jj]*x[kk,jj]+a01*x[kk,jj]+u01)+ exp(q11*x[kk,jj]*x[kk,jj]+a11*x[kk,jj]+u11)+ exp(q12*x[kk,jj]*x[kk,jj]+a12*x[kk,jj]+u12)+ exp(q02*x[kk,jj]*x[kk,jj]+a02*x[kk,jj]+u02)+ exp(q22*x[kk,jj]*x[kk,jj]+a22*x[kk,jj]+u22); p0_01=exp(b0*x[kk,jj]+c010)/(exp(b0*x[kk,jj]+c010)+1+exp(b2*x[kk,jj]+c012)); p01=exp(q01*x[kk,jj]*x[kk,jj]+a01*x[kk,jj]+u01)/total; p0_11=exp(b0*x[kk,jj]+c110)/(exp(b0*x[kk,jj]+c110)+1+exp(b2*x[kk,jj]+c112)); p11=exp(q11*x[kk,jj]*x[kk,jj]+a11*x[kk,jj]+u11)/total; rp0 = rp0 +.5*p0_01*p01/(.5*p0_01*p01+p0_11*p11); end; if (kid[kk,jj]=1) then do; total=exp(q00*x[kk,jj]*x[kk,jj]+a00*x[kk,jj]+u00)+ exp(q01*x[kk,jj]*x[kk,jj]+a01*x[kk,jj]+u01)+ exp(q11*x[kk,jj]*x[kk,jj]+a11*x[kk,jj]+u11)+ exp(q12*x[kk,jj]*x[kk,jj]+a12*x[kk,jj]+u12)+ exp(q02*x[kk,jj]*x[kk,jj]+a02*x[kk,jj]+u02)+ exp(q22*x[kk,jj]*x[kk,jj]+a22*x[kk,jj]+u22); p1_01=1/(exp(b0*x[kk,jj]+c010)+1+exp(b2*x[kk,jj]+c012)); p01=exp(q01*x[kk,jj]*x[kk,jj]+a01*x[kk,jj]+u01)/total; p1_11=1/(exp(b0*x[kk,jj]+c110)+1+exp(b2*x[kk,jj]+c112)); p11=exp(q11*x[kk,jj]*x[kk,jj]+a11*x[kk,jj]+u11)/total; rp0 = rp0 + .5*p1_01*p01/(.5*p1_01*p01+p1_11*p11); end; end; rp=(rp0)/3; if (r=rp) then missparent[kk]=1; notmissing[kk]=0; end; if (((obsparent[kk]=1) & (kid[kk,1]=0) & (kid[kk,2]=0) & (kid[kk,3]=2)) | ((obsparent[kk]=1) & (kid[kk,1]=0) & (kid[kk,2]=2) & (kid[kk,3]=0)) | ((obsparent[kk]=1) & (kid[kk,1]=2) & (kid[kk,2]=0) & (kid[kk,3]=0)) | ((obsparent[kk]=1) & (kid[kk,1]=0) & (kid[kk,2]=2) & (kid[kk,3]=2)) | ((obsparent[kk]=1) & (kid[kk,1]=2) & (kid[kk,2]=0) & (kid[kk,3]=2)) | ((obsparent[kk]=1) & (kid[kk,1]=2) & (kid[kk,2]=2) & (kid[kk,3]=0))) then do; missparent[kk]=1; notmissing[kk]=0; end; if (((obsparent[kk]=1) & (kid[kk,1]=0) & (kid[kk,2]=1) & (kid[kk,3]=2)) | ((obsparent[kk]=1) & (kid[kk,1]=1) & (kid[kk,2]=0) & (kid[kk,3]=2)) | ((obsparent[kk]=1) & (kid[kk,1]=0) & (kid[kk,2]=2) & (kid[kk,3]=1)) | ((obsparent[kk]=1) & (kid[kk,1]=1) & (kid[kk,2]=2) & (kid[kk,3]=0)) | ((obsparent[kk]=1) & (kid[kk,1]=2) & (kid[kk,2]=0) & (kid[kk,3]=1)) | ((obsparent[kk]=1) & (kid[kk,1]=2) & (kid[kk,2]=1) & (kid[kk,3]=0))) then do; missparent[kk]=1; notmissing[kk]=0; end; if ((obsparent[kk]=0) & (kid[kk,1]=1) & (kid[kk,2]=1) & (kid[kk,3]=1)) then do; r=RanUni(0); rp0=0; do jj=1 to nkids[kk]; p1_02=1/(exp(b0*x[kk,jj]+c020)+1+exp(b2*x[kk,jj]+c022)); p1_01=1/(exp(b0*x[kk,jj]+c010)+1+exp(b2*x[kk,jj]+c012)); total=exp(q00*x[kk,jj]*x[kk,jj]+a00*x[kk,jj]+u00)+ exp(q01*x[kk,jj]*x[kk,jj]+a01*x[kk,jj]+u01)+ exp(q11*x[kk,jj]*x[kk,jj]+a11*x[kk,jj]+u11)+ exp(q12*x[kk,jj]*x[kk,jj]+a12*x[kk,jj]+u12)+ exp(q02*x[kk,jj]*x[kk,jj]+a02*x[kk,jj]+u02)+ exp(q22*x[kk,jj]*x[kk,jj]+a22*x[kk,jj]+u22); p01=exp(q01*x[kk,jj]*x[kk,jj]+a01*x[kk,jj]+u01)/total; p02=exp(q02*x[kk,jj]*x[kk,jj]+a02*x[kk,jj]+u02)/total; rp0=rp0 + .5*p1_01*p01/(.5*p1_01*p01+.5*p1_02*p02); end; rp0=rp0/3; if (r=rp0) then missparent[kk]=2; notmissing[kk]=0; end; if ((obsparent[kk]=1) & (kid[kk,1]=1) & (kid[kk,2]=1) & (kid[kk,3]=1)) then do; r=RanUni(0); rp1=0; rp2=0; do jj=1 to nkids[kk]; p1_01=1/(exp(b0*x[kk,jj]+c010)+1+exp(b2*x[kk,jj]+c012)); p1_11=1/(exp(b0*x[kk,jj]+c110)+1+exp(b2*x[kk,jj]+c112)); p1_12=1/(exp(b0*x[kk,jj]+c120)+1+exp(b2*x[kk,jj]+c122)); total=exp(q00*x[kk,jj]*x[kk,jj]+a00*x[kk,jj]+u00)+ exp(q01*x[kk,jj]*x[kk,jj]+a01*x[kk,jj]+u01)+ exp(q11*x[kk,jj]*x[kk,jj]+a11*x[kk,jj]+u11)+ exp(q12*x[kk,jj]*x[kk,jj]+a12*x[kk,jj]+u12)+ exp(q02*x[kk,jj]*x[kk,jj]+a02*x[kk,jj]+u02)+ exp(q22*x[kk,jj]*x[kk,jj]+a22*x[kk,jj]+u22); p01=exp(q01*x[kk,jj]*x[kk,jj]+a01*x[kk,jj]+u01)/total; p11=exp(q11*x[kk,jj]*x[kk,jj]+a11*x[kk,jj]+u11)/total; p12=exp(q12*x[kk,jj]*x[kk,jj]+a12*x[kk,jj]+u12)/total; rp1= rp1 + .5*p1_01*p01/(.5*p1_01*p01+p1_11*p11+.5*p1_12*p12); rp2= rp2 + (.5*p1_01*p01+p1_11*p11)/(.5*p1_01*p01+p1_11*p11+.5*p1_12*p12); end; rp1=rp1/3; rp2=rp2/3; if (r=rp1) & (r=rp2) then missparent[kk]=2; notmissing[kk]=0; end; if ((obsparent[kk]=2) & (kid[kk,1]=1) & (kid[kk,2]=1) & (kid[kk,3]=1)) then do; r=RanUni(0); rp0=0; do jj=1 to nkids[kk]; p1_02=1/(exp(b0*x[kk,jj]+c020)+1+exp(b2*x[kk,jj]+c022)); p1_12=1/(exp(b0*x[kk,jj]+c120)+1+exp(b2*x[kk,jj]+c122)); total=exp(q00*x[kk,jj]*x[kk,jj]+a00*x[kk,jj]+u00)+ exp(q01*x[kk,jj]*x[kk,jj]+a01*x[kk,jj]+u01)+ exp(q11*x[kk,jj]*x[kk,jj]+a11*x[kk,jj]+u11)+ exp(q12*x[kk,jj]*x[kk,jj]+a12*x[kk,jj]+u12)+ exp(q02*x[kk,jj]*x[kk,jj]+a02*x[kk,jj]+u02)+ exp(q22*x[kk,jj]*x[kk,jj]+a22*x[kk,jj]+u22); p02=exp(q02*x[kk,jj]*x[kk,jj]+a02*x[kk,jj]+u02)/total; p12=exp(q12*x[kk,jj]*x[kk,jj]+a12*x[kk,jj]+u12)/total; rp0=rp0 + .5*p1_02*p02/(.5*p1_02*p02+.5*p1_12*p12); end; rp0=rp0/3; if (r=rp0) then missparent[kk]=1; notmissing[kk]=0; end; if (((obsparent[kk]=1) & (kid[kk,1]=1) & (kid[kk,2]=1) & (kid[kk,3]=2)) | ((obsparent[kk]=1) & (kid[kk,1]=1) & (kid[kk,2]=2) & (kid[kk,3]=1)) | ((obsparent[kk]=1) & (kid[kk,1]=2) & (kid[kk,2]=1) & (kid[kk,3]=1)) | ((obsparent[kk]=1) & (kid[kk,1]=1) & (kid[kk,2]=2) & (kid[kk,3]=2)) | ((obsparent[kk]=1) & (kid[kk,1]=2) & (kid[kk,2]=1) & (kid[kk,3]=2)) | ((obsparent[kk]=1) & (kid[kk,1]=2) & (kid[kk,2]=2) & (kid[kk,3]=1))) then do; r=RanUni(0); rp0=0; do jj=1 to nkids[kk]; if (kid[kk,jj]=1) then do; p1_11=1/(exp(b0*x[kk,jj]+c110)+1+exp(b2*x[kk,jj]+c112)); p1_12=1/(exp(b0*x[kk,jj]+c120)+1+exp(b2*x[kk,jj]+c122)); total=exp(q00*x[kk,jj]*x[kk,jj]+a00*x[kk,jj]+u00)+ exp(q01*x[kk,jj]*x[kk,jj]+a01*x[kk,jj]+u01)+ exp(q11*x[kk,jj]*x[kk,jj]+a11*x[kk,jj]+u11)+ exp(q12*x[kk,jj]*x[kk,jj]+a12*x[kk,jj]+u12)+ exp(q02*x[kk,jj]*x[kk,jj]+a02*x[kk,jj]+u02)+ exp(q22*x[kk,jj]*x[kk,jj]+a22*x[kk,jj]+u22); p11=exp(q11*x[kk,jj]*x[kk,jj]+a11*x[kk,jj]+u11)/total; p12=exp(q12*x[kk,jj]*x[kk,jj]+a12*x[kk,jj]+u12)/total; rp0= rp0 + (p1_11*p11)/(p1_11*p11+.5*p1_12*p12); end; if (kid[kk,jj]=2) then do; p2_11=exp(b2*x[kk,jj]+c112)/(exp(b0*x[kk,jj]+c110)+1+exp(b2*x[kk,jj]+c112)); p2_12=exp(b2*x[kk,jj]+c122)/(exp(b0*x[kk,jj]+c120)+1+exp(b2*x[kk,jj]+c122)); total=exp(q00*x[kk,jj]*x[kk,jj]+a00*x[kk,jj]+u00)+ exp(q01*x[kk,jj]*x[kk,jj]+a01*x[kk,jj]+u01)+ exp(q11*x[kk,jj]*x[kk,jj]+a11*x[kk,jj]+u11)+ exp(q12*x[kk,jj]*x[kk,jj]+a12*x[kk,jj]+u12)+ exp(q02*x[kk,jj]*x[kk,jj]+a02*x[kk,jj]+u02)+ exp(q22*x[kk,jj]*x[kk,jj]+a22*x[kk,jj]+u22); p11=exp(q11*x[kk,jj]*x[kk,jj]+a11*x[kk,jj]+u11)/total; p12=exp(q12*x[kk,jj]*x[kk,jj]+a12*x[kk,jj]+u12)/total; rp0= rp0 + p2_11*p11/(p2_11*p11+.5*p2_12*p12); end; end; rp0=rp0/3; if (r=rp0) then missparent[kk]=2; notmissing[kk]=0; end; if (((obsparent[kk]=2) & (kid[kk,1]=1) & (kid[kk,2]=1) & (kid[kk,3]=2)) | ((obsparent[kk]=2) & (kid[kk,1]=1) & (kid[kk,2]=2) & (kid[kk,3]=1)) | ((obsparent[kk]=2) & (kid[kk,1]=2) & (kid[kk,2]=1) & (kid[kk,3]=1)) | ((obsparent[kk]=2) & (kid[kk,1]=1) & (kid[kk,2]=2) & (kid[kk,3]=2)) | ((obsparent[kk]=2) & (kid[kk,1]=2) & (kid[kk,2]=1) & (kid[kk,3]=2)) | ((obsparent[kk]=2) & (kid[kk,1]=2) & (kid[kk,2]=2) & (kid[kk,3]=1))) then do; missparent[kk]=1; notmissing[kk]=0; end; if ((obsparent[kk]=1) & (kid[kk,1]=2) & (kid[kk,2]=2) & (kid[kk,3]=2)) then do; r=RanUni(0); rp0=0; do jj=1 to nkids[kk]; p2_11=exp(b2*x[kk,jj]+c112)/(exp(b0*x[kk,jj]+c110)+1+exp(b2*x[kk,jj]+c112)); p2_12=exp(b2*x[kk,jj]+c122)/(exp(b0*x[kk,jj]+c120)+1+exp(b2*x[kk,jj]+c122)); total=exp(q00*x[kk,jj]*x[kk,jj]+a00*x[kk,jj]+u00)+ exp(q01*x[kk,jj]*x[kk,jj]+a01*x[kk,jj]+u01)+ exp(q11*x[kk,jj]*x[kk,jj]+a11*x[kk,jj]+u11)+ exp(q12*x[kk,jj]*x[kk,jj]+a12*x[kk,jj]+u12)+ exp(q02*x[kk,jj]*x[kk,jj]+a02*x[kk,jj]+u02)+ exp(q22*x[kk,jj]*x[kk,jj]+a22*x[kk,jj]+u22); p11=exp(q11*x[kk,jj]*x[kk,jj]+a11*x[kk,jj]+u11)/total; p12=exp(q12*x[kk,jj]*x[kk,jj]+a12*x[kk,jj]+u12)/total; rp0=rp0+p2_11*p11/(p2_11*p11+.5*p2_12*p12); end; rp=rp0/3; if (r=rp) then missparent[kk]=2; notmissing[kk]=0; end; if ((obsparent[kk]=2) & (kid[kk,1]=2) & (kid[kk,2]=2) & (kid[kk,3]=2)) then do; r=RanUni(0); rp0=0; do jj=1 to nkids[kk]; p2_12=exp(b2*x[kk,jj]+c122)/(exp(b0*x[kk,jj]+c120)+1+exp(b2*x[kk,jj]+c122)); p2_22=exp(b2*x[kk,jj]+c222)/(exp(b0*x[kk,jj]+c220)+1+exp(b2*x[kk,jj]+c222)); total=exp(q00*x[kk,jj]*x[kk,jj]+a00*x[kk,jj]+u00)+ exp(q01*x[kk,jj]*x[kk,jj]+a01*x[kk,jj]+u01)+ exp(q11*x[kk,jj]*x[kk,jj]+a11*x[kk,jj]+u11)+ exp(q12*x[kk,jj]*x[kk,jj]+a12*x[kk,jj]+u12)+ exp(q02*x[kk,jj]*x[kk,jj]+a02*x[kk,jj]+u02)+ exp(q22*x[kk,jj]*x[kk,jj]+a22*x[kk,jj]+u22); p12=exp(q12*x[kk,jj]*x[kk,jj]+a12*x[kk,jj]+u12)/total; p22=exp(q22*x[kk,jj]*x[kk,jj]+a22*x[kk,jj]+u22)/total; rp0=rp0+.5*p2_12*p12/(.5*p2_12*p12+p2_22*p22); end; rp=rp0/3; if (r=rp) then missparent[kk]=2; notmissing[kk]=0; end; end; if ((obsparent[kk]=0) & (missparent[kk]=0)) then mfamily[kk]=1111; if (((obsparent[kk]=0) & (missparent[kk]=2)) | ((obsparent[kk]=2) & (missparent[kk]=0))) then mfamily[kk]=1122; if ((obsparent[kk]=2) & (missparent[kk]=2)) then mfamily[kk]=2222; if (((obsparent[kk]=0) & (missparent[kk]=1)) | ((obsparent[kk]=1) & (missparent[kk]=0))) then mfamily[kk]=1112; if ((obsparent[kk]=1) & (missparent[kk]=1)) then mfamily[kk]=1212; if (((obsparent[kk]=1) & (missparent[kk]=2)) | ((obsparent[kk]=2) & (missparent[kk]=1))) then mfamily[kk]=1222; if (flag_f[kk]=1) then do; missf[kk]=missparent[kk]; missm[kk]=obsparent[kk]; end; if (flag_m[kk]=1) then do; missm[kk]=missparent[kk]; missf[kk]=obsparent[kk]; end; end; end; return; FINISH MISS; /***********************************************************/ START add_will(x,kid,family,nkids,missing,additive,addlambda,addcovwill); count=0; sample=nrow(family); addlambda=j(additive,1,0); addcovwill=j(additive,additive,0); p=j(2,1,0); dl=j(additive,1,0); addvars=j(additive,additive,0); addhs = j(additive,additive,0); v = j(additive,additive,0); h = j(additive,additive,0); hfill = j(additive,additive,0); dlbydlt = j(additive,additive,0); covlambda = j(additive,additive,0); halfcov=j(additive,additive,0); do kk=1 to sample; do ii=1 to nkids[kk]; if (missing[kk]=0) then do; if ((family[kk]=1112)|(family[kk]=1212)|(family[kk]=1222)) then count=count+1; end; end; end; likepre=0.0; loglike=0.0; do ll=1 to 15; addparms=j(additive,1,0); addhs=j(additive,additive,0); addvars=j(additive,additive,0); v=j(additive,additive,0); covlambda=j(additive,additive,0); pluslambda=j(additive,1,0); if (ll=1) then do; b=addlambda[1]; a01=addlambda[2]; a11=addlambda[3]; a12=addlambda[4]; loglike=0; do kk=1 to sample; do ii=1 to nkids[kk]; if (missing[kk]=0) then do; if ((family[kk]=1112)|(family[kk]=1212)|(family[kk]=1222)) then do; num0=1; num1=exp(b*x[kk,ii]+a01*(family[kk]=1112)+a11*(family[kk]=1212)+a12*(family[kk]=1222)); denom=(num0+num1); p[1]=num0/denom; p[2]=num1/denom; loglike= loglike + log(((family[kk]=1112)*(kid[kk,ii]=0)*p[1])+ ((family[kk]=1112)*(kid[kk,ii]=1)*p[2])+ ((family[kk]=1212)*(kid[kk,ii]=0)*p[1]*p[1])+ ((family[kk]=1212)*(kid[kk,ii]=1)*p[1]*p[2])+ ((family[kk]=1212)*(kid[kk,ii]=2)*p[2]*p[2])+ ((family[kk]=1222)*(kid[kk,ii]=1)*p[1])+ ((family[kk]=1222)*(kid[kk,ii]=2)*p[2])); end; end; end; end; end; lastone=0; do kk=1 to sample; do ii=1 to nkids[kk]; if (missing[kk]=0) then do; if ((family[kk]=1112)|(family[kk]=1212)|(family[kk]=1222)) then do; lastone=lastone+1; num0=1; num1=exp(b*x[kk,ii]+a01*(family[kk]=1112)+a11*(family[kk]=1212)+a12*(family[kk]=1222)); denom=(num0+num1); p[1]=num0/denom; p[2]=num1/denom; if (family[kk]=1112) then do; dl[1]=x[kk,ii]*(kid[kk,ii]=1)-x[kk,ii]*p[2]; dl[2]=((kid[kk,ii]=1)-p[2]); dl[3]=0; dl[4]=0; h[1,1]=(x[kk,ii]*x[kk,ii])*p[2]*(p[2]-1); h[1,2]=(x[kk,ii])*p[2]*(p[2]-1); h[1,3]=0; h[1,4]=0; h[2,2]=p[2]*(p[2]-1); h[2,3]=0; h[2,4]=0; h[3,3]=0; h[3,4]=0; h[4,4]=0; end; if (family[kk]=1212) then do; dl[1]=x[kk,ii]*(kid[kk,ii]=1)+2*x[kk,ii]*(kid[kk,ii]=2)-2*x[kk,ii]*p[2]; dl[2]=0; dl[3]=((kid[kk,ii]=1)+2*(kid[kk,ii]=2)-2*p[2]); dl[4]=0; h[1,1]=2*(x[kk,ii]*x[kk,ii])*p[2]*(p[2]-1); h[1,2]=0; h[1,3]=2*(x[kk,ii])*p[2]*(p[2]-1); h[1,4]=0; h[2,2]=0; h[2,3]=0; h[2,4]=0; h[3,3]=2*p[2]*(p[2]-1); h[3,4]=0; h[4,4]=0; end; if (family[kk]=1222) then do; dl[1]=x[kk,ii]*(kid[kk,ii]=2)-x[kk,ii]*p[2]; dl[2]=0; dl[3]=0; dl[4]=((kid[kk,ii]=2)-p[2]); h[1,1]=(x[kk,ii]*x[kk,ii])*p[2]*(p[2]-1); h[1,2]=0; h[1,3]=0; h[1,4]=(x[kk,ii])*p[2]*(p[2]-1); h[2,2]=0; h[2,3]=0; h[2,4]=0; h[3,3]=0; h[3,4]=0; h[4,4]=p[2]*(p[2]-1); end; do i=1 to additive; do j=1 to additive; if (j=i) then hfill[j,i]=h[i,j]; end; end; do i=1 to additive; dl[i]=(1/nkids[kk])*dl[i]; do j=1 to additive; hfill[i,j]=(1/nkids[kk])*hfill[i,j]; end; end; addparms=addparms+dl; addhs=addhs-hfill; if (ii=1) then do; v=v+addvars; addvars=dl*dl`; end; if (ii^=1) then do; dlbydlt=dl*dl`; addvars=addvars+dlbydlt; end; if (lastone=count) then do; v=v+addvars; covlambda=inv(addhs); pluslambda=covlambda*addparms; addlambda=addlambda+pluslambda; halfcov=covlambda*v; addcovwill=halfcov*covlambda; flagit=b-addlambda[1]+ a01-addlambda[2]+ a11-addlambda[3]+ a12-addlambda[4]; end; end; end; end; end; likepre=loglike; loglike=0; oldb=b; olda01=a01; olda11=a11; olda12=a12; b=addlambda[1]; a01=addlambda[2]; a11=addlambda[3]; a12=addlambda[4]; do jj=1 to 15; if(jj=1) then do; loglike=0; do kk=1 to sample; do ii=1 to nkids[kk]; if (missing[kk]=0) then do; if ((family[kk]=1112)|(family[kk]=1212)|(family[kk]=1222)) then do; num0=1; num1=exp(b*x[kk,ii]+a01*(family[kk]=1112)+a11*(family[kk]=1212)+a12*(family[kk]=1222)); denom=(num0+num1); p[1]=num0/denom; p[2]=num1/denom; loglike=loglike + log(((family[kk]=1112)*(kid[kk,ii]=0)*p[1])+ ((family[kk]=1112)*(kid[kk,ii]=1)*p[2])+ ((family[kk]=1212)*(kid[kk,ii]=0)*p[1]*p[1])+ ((family[kk]=1212)*(kid[kk,ii]=1)*p[1]*p[2])+ ((family[kk]=1212)*(kid[kk,ii]=2)*p[2]*p[2])+ ((family[kk]=1222)*(kid[kk,ii]=1)*p[1])+ ((family[kk]=1222)*(kid[kk,ii]=2)*p[2])); end; end; end; end; end; if (loglike=rp1) then missparent[kk]=1; notmissing[kk]=0; end; if ((obsparent[kk]=0) & (kid[kk,1]=1)) then do; r=RanUni(0); rp1=.5*p1_01*p01/(.5*p1_01*p01+.5*p1_02*p02); if (r=rp1) then missparent[kk]=2; notmissing[kk]=0; end; if ((obsparent[kk]=1) & (kid[kk,1]=0)) then do; r=RanUni(0); rp1=.5*p0_01*p01/(.5*p0_01*p01+p0_11*p11); if (r=rp1) then missparent[kk]=1; notmissing[kk]=0; end; if ((obsparent[kk]=1) & (kid[kk,1]=1)) then do; r=RanUni(0); rp1=.5*p1_01*p01/(.5*p1_01*p01+p1_11*p11+.5*p1_12*p12); rp2=(.5*p1_01*p01+p1_11*p11)/(.5*p1_01*p01+p1_11*p11+.5*p1_12*p12); if (r=rp1) & (r=rp2) then missparent[kk]=2; notmissing[kk]=0; end; if ((obsparent[kk]=1) & (kid[kk,1]=2)) then do; r=RanUni(0); rp1=p2_11*p11/(p2_11*p11+.5*p2_12*p12); if (r=rp1) then missparent[kk]=2; notmissing[kk]=0; end; if ((obsparent[kk]=2) & (kid[kk,1]=1)) then do; r=RanUni(0); rp1=.5*p1_02*p02/(.5*p1_02*p02+.5*p1_12*p12); if (r=rp1) then missparent[kk]=1; notmissing[kk]=0; end; if ((obsparent[kk]=2) & (kid[kk,1]=2)) then do; r=RanUni(0); rp1=.5*p2_12*p12/(.5*p2_12*p12+p2_22*p22); if (r=rp1) then missparent[kk]=2; notmissing[kk]=0; end; end; if (nkids[kk]=2) then do; if ((obsparent[kk]=0) & (kid[kk,1]=0) & (kid[kk,2]=0)) then do; r=RanUni(0); rp0=0; do jj=1 to nkids[kk]; p0_00=1; p0_01=1/(exp(b*x[kk,jj]+c01)+1); total=exp(q00*x[kk,jj]*x[kk,jj]+a00*x[kk,jj]+u00)+ exp(q01*x[kk,jj]*x[kk,jj]+a01*x[kk,jj]+u01)+ exp(q11*x[kk,jj]*x[kk,jj]+a11*x[kk,jj]+u11)+ exp(q12*x[kk,jj]*x[kk,jj]+a12*x[kk,jj]+u12)+ exp(q02*x[kk,jj]*x[kk,jj]+a02*x[kk,jj]+u02)+ exp(q22*x[kk,jj]*x[kk,jj]+a22*x[kk,jj]+u22); p00=exp(q00*x[kk,jj]*x[kk,jj]+a00*x[kk,jj]+u00)/total; p01=exp(q01*x[kk,jj]*x[kk,jj]+a01*x[kk,jj]+u01)/total; end; rp=rp0/2; if (r=rp) then missparent[kk]=1; notmissing[kk]=0; end; if ((obsparent[kk]=1) & (kid[kk,1]=0) & (kid[kk,2]=0)) then do; r=RanUni(0); rp0=0; do jj=1 to nkids[kk]; total=exp(q00*x[kk,jj]*x[kk,jj]+a00*x[kk,jj]+u00)+ exp(q01*x[kk,jj]*x[kk,jj]+a01*x[kk,jj]+u01)+ exp(q11*x[kk,jj]*x[kk,jj]+a11*x[kk,jj]+u11)+ exp(q12*x[kk,jj]*x[kk,jj]+a12*x[kk,jj]+u12)+ exp(q02*x[kk,jj]*x[kk,jj]+a02*x[kk,jj]+u02)+ exp(q22*x[kk,jj]*x[kk,jj]+a22*x[kk,jj]+u22); p01=exp(q01*x[kk,jj]*x[kk,jj]+a01*x[kk,jj]+u01)/total; p11=exp(q11*x[kk,jj]*x[kk,jj]+a11*x[kk,jj]+u11)/total; p0_01=1/(exp(b*x[kk,jj]+c01)+1); p0_11=1/((exp(b*x[kk,jj]+c11)+1)*(exp(b*x[kk,jj]+c11)+1)); rp0=rp0 + .5*p0_01*p01/(.5*p0_01*p01+p0_11*p11); end; rp=rp0/2; if (r=rp) then missparent[kk]=1; notmissing[kk]=0; end; if (((obsparent[kk]=0) & (kid[kk,1]=0) & (kid[kk,2]=1)) | ((obsparent[kk]=0) & (kid[kk,1]=1) & (kid[kk,2]=0))) then do; missparent[kk]=1; notmissing[kk]=0; end; if (((obsparent[kk]=1) & (kid[kk,1]=0) & (kid[kk,2]=1)) | ((obsparent[kk]=1) & (kid[kk,1]=1) & (kid[kk,2]=0))) then do; r=RanUni(0); rp0=0; do jj=1 to nkids[kk]; if (kid[kk,jj]=0) then do; total=exp(q00*x[kk,jj]*x[kk,jj]+a00*x[kk,jj]+u00)+ exp(q01*x[kk,jj]*x[kk,jj]+a01*x[kk,jj]+u01)+ exp(q11*x[kk,jj]*x[kk,jj]+a11*x[kk,jj]+u11)+ exp(q12*x[kk,jj]*x[kk,jj]+a12*x[kk,jj]+u12)+ exp(q02*x[kk,jj]*x[kk,jj]+a02*x[kk,jj]+u02)+ exp(q22*x[kk,jj]*x[kk,jj]+a22*x[kk,jj]+u22); p0_01=1/(exp(b*x[kk,jj]+c01)+1); p01=exp(q01*x[kk,jj]*x[kk,jj]+a01*x[kk,jj]+u01)/total; p0_11=1/((exp(b*x[kk,jj]+c11)+1)*(exp(b*x[kk,jj]+c11)+1)); p11=exp(q11*x[kk,jj]*x[kk,jj]+a11*x[kk,jj]+u11)/total; rp0=rp0+.5*p0_01*p01/(.5*p0_01*p01+p0_11*p11); end; if (kid[kk,jj]=1) then do; total=exp(q00*x[kk,jj]*x[kk,jj]+a00*x[kk,jj]+u00)+ exp(q01*x[kk,jj]*x[kk,jj]+a01*x[kk,jj]+u01)+ exp(q11*x[kk,jj]*x[kk,jj]+a11*x[kk,jj]+u11)+ exp(q12*x[kk,jj]*x[kk,jj]+a12*x[kk,jj]+u12)+ exp(q02*x[kk,jj]*x[kk,jj]+a02*x[kk,jj]+u02)+ exp(q22*x[kk,jj]*x[kk,jj]+a22*x[kk,jj]+u22); p1_01=exp(b*x[kk,jj]+c01)/(exp(b*x[kk,jj]+c01)+1); p01=exp(q01*x[kk,jj]*x[kk,jj]+a01*x[kk,jj]+u01)/total; p1_11=2*exp(b*x[kk,jj]+c11)/((exp(b*x[kk,jj]+c11)+1)*(exp(b*x[kk,jj]+c11)+1)); p11=exp(q11*x[kk,jj]*x[kk,jj]+a11*x[kk,jj]+u11)/total; rp0=rp0+.5*p1_01*p01/(.5*p1_01*p01+p1_11*p11); end; end; rp=rp0/2; if (r=rp) then missparent[kk]=1; notmissing[kk]=0; end; if (((obsparent[kk]=1) & (kid[kk,1]=0) & (kid[kk,2]=2)) | ((obsparent[kk]=1) & (kid[kk,1]=2) & (kid[kk,2]=0))) then do; missparent[kk]=1; notmissing[kk]=0; end; if ((obsparent[kk]=0) & (kid[kk,1]=1) & (kid[kk,2]=1)) then do; r=RanUni(0); rp0=0; do jj=1 to nkids[kk]; p1_02=1; p1_01=exp(b*x[kk,jj]+c01)/(exp(b*x[kk,jj]+c01)+1); total=exp(q00*x[kk,jj]*x[kk,jj]+a00*x[kk,jj]+u00)+ exp(q01*x[kk,jj]*x[kk,jj]+a01*x[kk,jj]+u01)+ exp(q11*x[kk,jj]*x[kk,jj]+a11*x[kk,jj]+u11)+ exp(q12*x[kk,jj]*x[kk,jj]+a12*x[kk,jj]+u12)+ exp(q02*x[kk,jj]*x[kk,jj]+a02*x[kk,jj]+u02)+ exp(q22*x[kk,jj]*x[kk,jj]+a22*x[kk,jj]+u22); p01=exp(q01*x[kk,jj]*x[kk,jj]+a01*x[kk,jj]+u01)/total; p02=exp(q02*x[kk,jj]*x[kk,jj]+a02*x[kk,jj]+u02)/total; rp0=rp0+.5*p1_01*p01/(.5*p1_01*p01+.5*p1_02*p02); end; rp=rp0/2; if (r=rp) then missparent[kk]=2; notmissing[kk]=0; end; if ((obsparent[kk]=1) & (kid[kk,1]=1) & (kid[kk,2]=1)) then do; r=RanUni(0); rp1=0; rp2=0; do jj=1 to nkids[kk]; p1_01=exp(b*x[kk,jj]+c01)/(exp(b*x[kk,jj]+c01)+1); p1_11=2*exp(b*x[kk,jj]+c11)/((exp(b*x[kk,jj]+c11)+1)*(exp(b*x[kk,jj]+c11)+1)); p1_12=1/(exp(b*x[kk,jj]+c12)+1); total=exp(q00*x[kk,jj]*x[kk,jj]+a00*x[kk,jj]+u00)+ exp(q01*x[kk,jj]*x[kk,jj]+a01*x[kk,jj]+u01)+ exp(q11*x[kk,jj]*x[kk,jj]+a11*x[kk,jj]+u11)+ exp(q12*x[kk,jj]*x[kk,jj]+a12*x[kk,jj]+u12)+ exp(q02*x[kk,jj]*x[kk,jj]+a02*x[kk,jj]+u02)+ exp(q22*x[kk,jj]*x[kk,jj]+a22*x[kk,jj]+u22); p01=exp(q01*x[kk,jj]*x[kk,jj]+a01*x[kk,jj]+u01)/total; p11=exp(q11*x[kk,jj]*x[kk,jj]+a11*x[kk,jj]+u11)/total; p12=exp(q12*x[kk,jj]*x[kk,jj]+a12*x[kk,jj]+u12)/total; rp1= rp1 + .5*p1_01*p01/(.5*p1_01*p01+p1_11*p11+.5*p1_12*p12); rp2= rp2 + (.5*p1_01*p01+p1_11*p11)/(.5*p1_01*p01+p1_11*p11+.5*p1_12*p12); end; rp1=rp1/2; rp2=rp2/2; if (r=rp1) & (r=rp2) then missparent[kk]=2; notmissing[kk]=0; end; if ((obsparent[kk]=2) & (kid[kk,1]=1) & (kid[kk,2]=1)) then do; r=RanUni(0); rp0=0; do jj=1 to nkids[kk]; p1_02=1; p1_12=1/(exp(b*x[kk,jj]+c12)+1); total=exp(q00*x[kk,jj]*x[kk,jj]+a00*x[kk,jj]+u00)+ exp(q01*x[kk,jj]*x[kk,jj]+a01*x[kk,jj]+u01)+ exp(q11*x[kk,jj]*x[kk,jj]+a11*x[kk,jj]+u11)+ exp(q12*x[kk,jj]*x[kk,jj]+a12*x[kk,jj]+u12)+ exp(q02*x[kk,jj]*x[kk,jj]+a02*x[kk,jj]+u02)+ exp(q22*x[kk,jj]*x[kk,jj]+a22*x[kk,jj]+u22); p12=exp(q12*x[kk,jj]*x[kk,jj]+a12*x[kk,jj]+u12)/total; p02=exp(q02*x[kk,jj]*x[kk,jj]+a02*x[kk,jj]+u02)/total; rp0= rp0 + .5*p1_02*p02/(.5*p1_02*p02+.5*p1_12*p12); end; rp=rp0/2; if (r=rp) then missparent[kk]=1; notmissing[kk]=0; end; if (((obsparent[kk]=1) & (kid[kk,1]=1) & (kid[kk,2]=2)) | ((obsparent[kk]=1) & (kid[kk,1]=2) & (kid[kk,2]=1))) then do; r=RanUni(0); rp0=0; do jj=1 to nkids[kk]; if (kid[kk,jj]=1) then do; p1_11=2*exp(b*x[kk,jj]+c11)/((exp(b*x[kk,jj]+c11)+1)*(exp(b*x[kk,jj]+c11)+1)); p1_12=1/(exp(b*x[kk,jj]+c12)+1); total=exp(q00*x[kk,jj]*x[kk,jj]+a00*x[kk,jj]+u00)+ exp(q01*x[kk,jj]*x[kk,jj]+a01*x[kk,jj]+u01)+ exp(q11*x[kk,jj]*x[kk,jj]+a11*x[kk,jj]+u11)+ exp(q12*x[kk,jj]*x[kk,jj]+a12*x[kk,jj]+u12)+ exp(q02*x[kk,jj]*x[kk,jj]+a02*x[kk,jj]+u02)+ exp(q22*x[kk,jj]*x[kk,jj]+a22*x[kk,jj]+u22); p11=exp(q11*x[kk,jj]*x[kk,jj]+a11*x[kk,jj]+u11)/total; p12=exp(q12*x[kk,jj]*x[kk,jj]+a12*x[kk,jj]+u12)/total; rp0=rp0+(p1_11*p11)/(p1_11*p11+.5*p1_12*p12); end; if (kid[kk,jj]=2) then do; p2_11=(exp(b*x[kk,jj]+c11)*exp(b*x[kk,jj]+c11))/ ((exp(b*x[kk,jj]+c11)+1)*(exp(b*x[kk,jj]+c11)+1)); p2_12=exp(b*x[kk,jj]+c12)/(exp(b*x[kk,jj]+c12)+1); total=exp(q00*x[kk,jj]*x[kk,jj]+a00*x[kk,jj]+u00)+ exp(q01*x[kk,jj]*x[kk,jj]+a01*x[kk,jj]+u01)+ exp(q11*x[kk,jj]*x[kk,jj]+a11*x[kk,jj]+u11)+ exp(q12*x[kk,jj]*x[kk,jj]+a12*x[kk,jj]+u12)+ exp(q02*x[kk,jj]*x[kk,jj]+a02*x[kk,jj]+u02)+ exp(q22*x[kk,jj]*x[kk,jj]+a22*x[kk,jj]+u22); p11=exp(q11*x[kk,jj]*x[kk,jj]+a11*x[kk,jj]+u11)/total; p12=exp(q12*x[kk,jj]*x[kk,jj]+a12*x[kk,jj]+u12)/total; rp0=rp0+p2_11*p11/(p2_11*p11+.5*p2_12*p12); end; end; rp=rp0/2; if (r=rp) then missparent[kk]=2; notmissing[kk]=0; end; if (((obsparent[kk]=2) & (kid[kk,1]=1) & (kid[kk,2]=2)) | ((obsparent[kk]=2) & (kid[kk,1]=2) & (kid[kk,2]=1))) then do; missparent[kk]=1; notmissing[kk]=0; end; if ((obsparent[kk]=1) & (kid[kk,1]=2) & (kid[kk,2]=2)) then do; r=RanUni(0); rp0=0; do jj=1 to nkids[kk]; p2_11=(exp(b*x[kk,jj]+c11)*exp(b*x[kk,jj]+c11))/ ((exp(b*x[kk,jj]+c11)+1)*(exp(b*x[kk,jj]+c11)+1)); p2_12=exp(b*x[kk,jj]+c12)/(exp(b*x[kk,jj]+c12)+1); total=exp(q00*x[kk,jj]*x[kk,jj]+a00*x[kk,jj]+u00)+ exp(q01*x[kk,jj]*x[kk,jj]+a01*x[kk,jj]+u01)+ exp(q11*x[kk,jj]*x[kk,jj]+a11*x[kk,jj]+u11)+ exp(q12*x[kk,jj]*x[kk,jj]+a12*x[kk,jj]+u12)+ exp(q02*x[kk,jj]*x[kk,jj]+a02*x[kk,jj]+u02)+ exp(q22*x[kk,jj]*x[kk,jj]+a22*x[kk,jj]+u22); p11=exp(q11*x[kk,jj]*x[kk,jj]+a11*x[kk,jj]+u11)/total; p12=exp(q12*x[kk,jj]*x[kk,jj]+a12*x[kk,jj]+u12)/total; rp0=rp0+p2_11*p11/(p2_11*p11+.5*p2_12*p12); end; rp=rp0/2; if (r=rp) then missparent[kk]=2; notmissing[kk]=0; end; if ((obsparent[kk]=2) & (kid[kk,1]=2) & (kid[kk,2]=2)) then do; r=RanUni(0); rp0=0; do jj=1 to nkids[kk]; p2_12=exp(b*x[kk,jj]+c12)/(exp(b*x[kk,jj]+c12)+1); p2_22=1; total=exp(q00*x[kk,jj]*x[kk,jj]+a00*x[kk,jj]+u00)+ exp(q01*x[kk,jj]*x[kk,jj]+a01*x[kk,jj]+u01)+ exp(q11*x[kk,jj]*x[kk,jj]+a11*x[kk,jj]+u11)+ exp(q12*x[kk,jj]*x[kk,jj]+a12*x[kk,jj]+u12)+ exp(q02*x[kk,jj]*x[kk,jj]+a02*x[kk,jj]+u02)+ exp(q22*x[kk,jj]*x[kk,jj]+a22*x[kk,jj]+u22); p12=exp(q12*x[kk,jj]*x[kk,jj]+a12*x[kk,jj]+u12)/total; p22=exp(q22*x[kk,jj]*x[kk,jj]+a22*x[kk,jj]+u22)/total; rp0=rp0+.5*p2_12*p12/(.5*p2_12*p12+p2_22*p22); end; rp=rp0/2; if (r=rp) then missparent[kk]=2; notmissing[kk]=0; end; end; if (nkids[kk]=3) then do; if ((obsparent[kk]=0) & (kid[kk,1]=0) & (kid[kk,2]=0) & (kid[kk,3]=0)) then do; r=RanUni(0); rp0=0; do jj=1 to nkids[kk]; p0_00=1; p0_01=1/(exp(b*x[kk,jj]+c01)+1); total=exp(q00*x[kk,jj]*x[kk,jj]+a00*x[kk,jj]+u00)+ exp(q01*x[kk,jj]*x[kk,jj]+a01*x[kk,jj]+u01)+ exp(q11*x[kk,jj]*x[kk,jj]+a11*x[kk,jj]+u11)+ exp(q12*x[kk,jj]*x[kk,jj]+a12*x[kk,jj]+u12)+ exp(q02*x[kk,jj]*x[kk,jj]+a02*x[kk,jj]+u02)+ exp(q22*x[kk,jj]*x[kk,jj]+a22*x[kk,jj]+u22); p00=exp(q00*x[kk,jj]*x[kk,jj]+a00*x[kk,jj]+u00)/total; p01=exp(q01*x[kk,jj]*x[kk,jj]+a01*x[kk,jj]+u01)/total; rp0=rp0+p0_00*p00/(p0_00*p00+.5*p0_01*p01); end; rp=rp0/3; if (r=rp) then missparent[kk]=1; notmissing[kk]=0; end; if ((obsparent[kk]=1) & (kid[kk,1]=0) & (kid[kk,2]=0) & (kid[kk,3]=0)) then do; r=RanUni(0); rp0=0; do jj=1 to nkids[kk]; total=exp(q00*x[kk,jj]*x[kk,jj]+a00*x[kk,jj]+u00)+ exp(q01*x[kk,jj]*x[kk,jj]+a01*x[kk,jj]+u01)+ exp(q11*x[kk,jj]*x[kk,jj]+a11*x[kk,jj]+u11)+ exp(q12*x[kk,jj]*x[kk,jj]+a12*x[kk,jj]+u12)+ exp(q02*x[kk,jj]*x[kk,jj]+a02*x[kk,jj]+u02)+ exp(q22*x[kk,jj]*x[kk,jj]+a22*x[kk,jj]+u22); p0_01=1/(exp(b*x[kk,jj]+c01)+1); p0_11=1/((exp(b*x[kk,jj]+c11)+1)*(exp(b*x[kk,jj]+c11)+1)); p01=exp(q01*x[kk,jj]*x[kk,jj]+a01*x[kk,jj]+u01)/total; p11=exp(q11*x[kk,jj]*x[kk,jj]+a11*x[kk,jj]+u11)/total; rp0= rp0 + .5*p0_01*p01/(.5*p0_01*p01+p0_11*p11); end; rp=rp0/3; if (r=rp) then missparent[kk]=1; notmissing[kk]=0; end; if (((obsparent[kk]=0) & (kid[kk,1]=0) & (kid[kk,2]=0) & (kid[kk,3]=1)) | ((obsparent[kk]=0) & (kid[kk,1]=0) & (kid[kk,2]=1) & (kid[kk,3]=0)) | ((obsparent[kk]=0) & (kid[kk,1]=1) & (kid[kk,2]=0) & (kid[kk,3]=0)) | ((obsparent[kk]=0) & (kid[kk,1]=0) & (kid[kk,2]=1) & (kid[kk,3]=1)) | ((obsparent[kk]=0) & (kid[kk,1]=1) & (kid[kk,2]=0) & (kid[kk,3]=1)) | ((obsparent[kk]=0) & (kid[kk,1]=1) & (kid[kk,2]=1) & (kid[kk,3]=0))) then do; missparent[kk]=1; notmissing[kk]=0; end; if (((obsparent[kk]=1) & (kid[kk,1]=0) & (kid[kk,2]=0) & (kid[kk,3]=1)) | ((obsparent[kk]=1) & (kid[kk,1]=0) & (kid[kk,2]=1) & (kid[kk,3]=0)) | ((obsparent[kk]=1) & (kid[kk,1]=1) & (kid[kk,2]=0) & (kid[kk,3]=0)) | ((obsparent[kk]=1) & (kid[kk,1]=0) & (kid[kk,2]=1) & (kid[kk,3]=1)) | ((obsparent[kk]=1) & (kid[kk,1]=1) & (kid[kk,2]=0) & (kid[kk,3]=1)) | ((obsparent[kk]=1) & (kid[kk,1]=1) & (kid[kk,2]=1) & (kid[kk,3]=0))) then do; r=RanUni(0); rp0=0; do jj=1 to nkids[kk]; if (kid[kk,jj]=0) then do; total=exp(q00*x[kk,jj]*x[kk,jj]+a00*x[kk,jj]+u00)+ exp(q01*x[kk,jj]*x[kk,jj]+a01*x[kk,jj]+u01)+ exp(q11*x[kk,jj]*x[kk,jj]+a11*x[kk,jj]+u11)+ exp(q12*x[kk,jj]*x[kk,jj]+a12*x[kk,jj]+u12)+ exp(q02*x[kk,jj]*x[kk,jj]+a02*x[kk,jj]+u02)+ exp(q22*x[kk,jj]*x[kk,jj]+a22*x[kk,jj]+u22); p0_01=1/(exp(b*x[kk,jj]+c01)+1); p01=exp(q01*x[kk,jj]*x[kk,jj]+a01*x[kk,jj]+u01)/total; p0_11=1/((exp(b*x[kk,jj]+c11)+1)*(exp(b*x[kk,jj]+c11)+1)); p11=exp(q11*x[kk,jj]*x[kk,jj]+a11*x[kk,jj]+u11)/total; rp0 = rp0 +.5*p0_01*p01/(.5*p0_01*p01+p0_11*p11); end; if (kid[kk,jj]=1) then do; total=exp(q00*x[kk,jj]*x[kk,jj]+a00*x[kk,jj]+u00)+ exp(q01*x[kk,jj]*x[kk,jj]+a01*x[kk,jj]+u01)+ exp(q11*x[kk,jj]*x[kk,jj]+a11*x[kk,jj]+u11)+ exp(q12*x[kk,jj]*x[kk,jj]+a12*x[kk,jj]+u12)+ exp(q02*x[kk,jj]*x[kk,jj]+a02*x[kk,jj]+u02)+ exp(q22*x[kk,jj]*x[kk,jj]+a22*x[kk,jj]+u22); p1_01=exp(b*x[kk,jj]+c01)/(exp(b*x[kk,jj]+c01)+1); p1_11=exp(b*x[kk,jj]+c11)/((exp(b*x[kk,jj]+c11)+1)*(exp(b*x[kk,jj]+c11)+1)); p01=exp(q01*x[kk,jj]*x[kk,jj]+a01*x[kk,jj]+u01)/total; p11=exp(q11*x[kk,jj]*x[kk,jj]+a11*x[kk,jj]+u11)/total; rp0 = rp0 + .5*p1_01*p01/(.5*p1_01*p01+p1_11*p11); end; end; rp=(rp0)/3; if (r=rp) then missparent[kk]=1; notmissing[kk]=0; end; if (((obsparent[kk]=1) & (kid[kk,1]=0) & (kid[kk,2]=0) & (kid[kk,3]=2)) | ((obsparent[kk]=1) & (kid[kk,1]=0) & (kid[kk,2]=2) & (kid[kk,3]=0)) | ((obsparent[kk]=1) & (kid[kk,1]=2) & (kid[kk,2]=0) & (kid[kk,3]=0)) | ((obsparent[kk]=1) & (kid[kk,1]=0) & (kid[kk,2]=2) & (kid[kk,3]=2)) | ((obsparent[kk]=1) & (kid[kk,1]=2) & (kid[kk,2]=0) & (kid[kk,3]=2)) | ((obsparent[kk]=1) & (kid[kk,1]=2) & (kid[kk,2]=2) & (kid[kk,3]=0))) then do; missparent[kk]=1; notmissing[kk]=0; end; if (((obsparent[kk]=1) & (kid[kk,1]=0) & (kid[kk,2]=1) & (kid[kk,3]=2)) | ((obsparent[kk]=1) & (kid[kk,1]=1) & (kid[kk,2]=0) & (kid[kk,3]=2)) | ((obsparent[kk]=1) & (kid[kk,1]=0) & (kid[kk,2]=2) & (kid[kk,3]=1)) | ((obsparent[kk]=1) & (kid[kk,1]=1) & (kid[kk,2]=2) & (kid[kk,3]=0)) | ((obsparent[kk]=1) & (kid[kk,1]=2) & (kid[kk,2]=0) & (kid[kk,3]=1)) | ((obsparent[kk]=1) & (kid[kk,1]=2) & (kid[kk,2]=1) & (kid[kk,3]=0))) then do; missparent[kk]=1; notmissing[kk]=0; end; if ((obsparent[kk]=0) & (kid[kk,1]=1) & (kid[kk,2]=1) & (kid[kk,3]=1)) then do; r=RanUni(0); rp0=0; do jj=1 to nkids[kk]; p1_02=1; p1_01=exp(b*x[kk,jj]+c01)/(exp(b*x[kk,jj]+c01)+1); total=exp(q00*x[kk,jj]*x[kk,jj]+a00*x[kk,jj]+u00)+ exp(q01*x[kk,jj]*x[kk,jj]+a01*x[kk,jj]+u01)+ exp(q11*x[kk,jj]*x[kk,jj]+a11*x[kk,jj]+u11)+ exp(q12*x[kk,jj]*x[kk,jj]+a12*x[kk,jj]+u12)+ exp(q02*x[kk,jj]*x[kk,jj]+a02*x[kk,jj]+u02)+ exp(q22*x[kk,jj]*x[kk,jj]+a22*x[kk,jj]+u22); p01=exp(q01*x[kk,jj]*x[kk,jj]+a01*x[kk,jj]+u01)/total; p02=exp(q02*x[kk,jj]*x[kk,jj]+a02*x[kk,jj]+u02)/total; rp0=rp0 + .5*p1_01*p01/(.5*p1_01*p01+.5*p1_02*p02); end; rp0=rp0/3; if (r=rp0) then missparent[kk]=2; notmissing[kk]=0; end; if ((obsparent[kk]=1) & (kid[kk,1]=1) & (kid[kk,2]=1) & (kid[kk,3]=1)) then do; r=RanUni(0); rp1=0; rp2=0; do jj=1 to nkids[kk]; p1_01=exp(b*x[kk,jj]+c01)/(exp(b*x[kk,jj]+c01)+1); p1_11=2*exp(b*x[kk,jj]+c11)/((exp(b*x[kk,jj]+c11)+1)*(exp(b*x[kk,jj]+c11)+1)); p1_12=1/(exp(b*x[kk,jj]+c12)+1); total=exp(q00*x[kk,jj]*x[kk,jj]+a00*x[kk,jj]+u00)+ exp(q01*x[kk,jj]*x[kk,jj]+a01*x[kk,jj]+u01)+ exp(q11*x[kk,jj]*x[kk,jj]+a11*x[kk,jj]+u11)+ exp(q12*x[kk,jj]*x[kk,jj]+a12*x[kk,jj]+u12)+ exp(q02*x[kk,jj]*x[kk,jj]+a02*x[kk,jj]+u02)+ exp(q22*x[kk,jj]*x[kk,jj]+a22*x[kk,jj]+u22); p01=exp(q01*x[kk,jj]*x[kk,jj]+a01*x[kk,jj]+u01)/total; p11=exp(q11*x[kk,jj]*x[kk,jj]+a11*x[kk,jj]+u11)/total; p12=exp(q12*x[kk,jj]*x[kk,jj]+a12*x[kk,jj]+u12)/total; rp1= rp1 + .5*p1_01*p01/(.5*p1_01*p01+p1_11*p11+.5*p1_12*p12); rp2= rp2 + (.5*p1_01*p01+p1_11*p11)/(.5*p1_01*p01+p1_11*p11+.5*p1_12*p12); end; rp1=rp1/3; rp2=rp2/3; if (r=rp1) & (r=rp2) then missparent[kk]=2; notmissing[kk]=0; end; if ((obsparent[kk]=2) & (kid[kk,1]=1) & (kid[kk,2]=1) & (kid[kk,3]=1)) then do; r=RanUni(0); rp0=0; do jj=1 to nkids[kk]; p1_02=1; p1_12=1/(exp(b*x[kk,jj]+c12)+1); total=exp(q00*x[kk,jj]*x[kk,jj]+a00*x[kk,jj]+u00)+ exp(q01*x[kk,jj]*x[kk,jj]+a01*x[kk,jj]+u01)+ exp(q11*x[kk,jj]*x[kk,jj]+a11*x[kk,jj]+u11)+ exp(q12*x[kk,jj]*x[kk,jj]+a12*x[kk,jj]+u12)+ exp(q02*x[kk,jj]*x[kk,jj]+a02*x[kk,jj]+u02)+ exp(q22*x[kk,jj]*x[kk,jj]+a22*x[kk,jj]+u22); p02=exp(q02*x[kk,jj]*x[kk,jj]+a02*x[kk,jj]+u02)/total; p12=exp(q12*x[kk,jj]*x[kk,jj]+a12*x[kk,jj]+u12)/total; rp0=rp0 + .5*p1_02*p02/(.5*p1_02*p02+.5*p1_12*p12); end; rp0=rp0/3; if (r=rp0) then missparent[kk]=1; notmissing[kk]=0; end; if (((obsparent[kk]=1) & (kid[kk,1]=1) & (kid[kk,2]=1) & (kid[kk,3]=2)) | ((obsparent[kk]=1) & (kid[kk,1]=1) & (kid[kk,2]=2) & (kid[kk,3]=1)) | ((obsparent[kk]=1) & (kid[kk,1]=2) & (kid[kk,2]=1) & (kid[kk,3]=1)) | ((obsparent[kk]=1) & (kid[kk,1]=1) & (kid[kk,2]=2) & (kid[kk,3]=2)) | ((obsparent[kk]=1) & (kid[kk,1]=2) & (kid[kk,2]=1) & (kid[kk,3]=2)) | ((obsparent[kk]=1) & (kid[kk,1]=2) & (kid[kk,2]=2) & (kid[kk,3]=1))) then do; r=RanUni(0); rp0=0; do jj=1 to nkids[kk]; if (kid[kk,jj]=1) then do; p1_11=2*exp(b*x[kk,jj]+c11)/((exp(b*x[kk,jj]+c11)+1)*(exp(b*x[kk,jj]+c11)+1)); p1_12=1/(exp(b*x[kk,jj]+c12)+1); total=exp(q00*x[kk,jj]*x[kk,jj]+a00*x[kk,jj]+u00)+ exp(q01*x[kk,jj]*x[kk,jj]+a01*x[kk,jj]+u01)+ exp(q11*x[kk,jj]*x[kk,jj]+a11*x[kk,jj]+u11)+ exp(q12*x[kk,jj]*x[kk,jj]+a12*x[kk,jj]+u12)+ exp(q02*x[kk,jj]*x[kk,jj]+a02*x[kk,jj]+u02)+ exp(q22*x[kk,jj]*x[kk,jj]+a22*x[kk,jj]+u22); p11=exp(q11*x[kk,jj]*x[kk,jj]+a11*x[kk,jj]+u11)/total; p12=exp(q12*x[kk,jj]*x[kk,jj]+a12*x[kk,jj]+u12)/total; rp0= rp0 + (p1_11*p11)/(p1_11*p11+.5*p1_12*p12); end; if (kid[kk,jj]=2) then do; p2_11=(exp(b*x[kk,jj]+c11)*exp(b*x[kk,jj]+c11))/ ((exp(b*x[kk,jj]+c11)+1)*(exp(b*x[kk,jj]+c11)+1)); p2_12=exp(b*x[kk,jj]+c12)/(exp(b*x[kk,jj]+c12)+1); total=exp(q00*x[kk,jj]*x[kk,jj]+a00*x[kk,jj]+u00)+ exp(q01*x[kk,jj]*x[kk,jj]+a01*x[kk,jj]+u01)+ exp(q11*x[kk,jj]*x[kk,jj]+a11*x[kk,jj]+u11)+ exp(q12*x[kk,jj]*x[kk,jj]+a12*x[kk,jj]+u12)+ exp(q02*x[kk,jj]*x[kk,jj]+a02*x[kk,jj]+u02)+ exp(q22*x[kk,jj]*x[kk,jj]+a22*x[kk,jj]+u22); p11=exp(q11*x[kk,jj]*x[kk,jj]+a11*x[kk,jj]+u11)/total; p12=exp(q12*x[kk,jj]*x[kk,jj]+a12*x[kk,jj]+u12)/total; rp0= rp0 + p2_11*p11/(p2_11*p11+.5*p2_12*p12); end; end; rp0=rp0/3; if (r=rp0) then missparent[kk]=2; notmissing[kk]=0; end; if (((obsparent[kk]=2) & (kid[kk,1]=1) & (kid[kk,2]=1) & (kid[kk,3]=2)) | ((obsparent[kk]=2) & (kid[kk,1]=1) & (kid[kk,2]=2) & (kid[kk,3]=1)) | ((obsparent[kk]=2) & (kid[kk,1]=2) & (kid[kk,2]=1) & (kid[kk,3]=1)) | ((obsparent[kk]=2) & (kid[kk,1]=1) & (kid[kk,2]=2) & (kid[kk,3]=2)) | ((obsparent[kk]=2) & (kid[kk,1]=2) & (kid[kk,2]=1) & (kid[kk,3]=2)) | ((obsparent[kk]=2) & (kid[kk,1]=2) & (kid[kk,2]=2) & (kid[kk,3]=1))) then do; missparent[kk]=1; notmissing[kk]=0; end; if ((obsparent[kk]=1) & (kid[kk,1]=2) & (kid[kk,2]=2) & (kid[kk,3]=2)) then do; r=RanUni(0); rp0=0; do jj=1 to nkids[kk]; p2_11=(exp(b*x[kk,jj]+c11)*exp(b*x[kk,jj]+c11))/ ((exp(b*x[kk,jj]+c11)+1)*(exp(b*x[kk,jj]+c11)+1)); p2_12=exp(b*x[kk,jj]+c12)/(exp(b*x[kk,jj]+c12)+1); total=exp(q00*x[kk,jj]*x[kk,jj]+a00*x[kk,jj]+u00)+ exp(q01*x[kk,jj]*x[kk,jj]+a01*x[kk,jj]+u01)+ exp(q11*x[kk,jj]*x[kk,jj]+a11*x[kk,jj]+u11)+ exp(q12*x[kk,jj]*x[kk,jj]+a12*x[kk,jj]+u12)+ exp(q02*x[kk,jj]*x[kk,jj]+a02*x[kk,jj]+u02)+ exp(q22*x[kk,jj]*x[kk,jj]+a22*x[kk,jj]+u22); p11=exp(q11*x[kk,jj]*x[kk,jj]+a11*x[kk,jj]+u11)/total; p12=exp(q12*x[kk,jj]*x[kk,jj]+a12*x[kk,jj]+u12)/total; rp0=rp0+p2_11*p11/(p2_11*p11+.5*p2_12*p12); end; rp=rp0/3; if (r=rp) then missparent[kk]=2; notmissing[kk]=0; end; if ((obsparent[kk]=2) & (kid[kk,1]=2) & (kid[kk,2]=2) & (kid[kk,3]=2)) then do; r=RanUni(0); rp0=0; do jj=1 to nkids[kk]; p2_12=exp(b*x[kk,jj]+c12)/(exp(b*x[kk,jj]+c12)+1); p2_22=1; total=exp(q00*x[kk,jj]*x[kk,jj]+a00*x[kk,jj]+u00)+ exp(q01*x[kk,jj]*x[kk,jj]+a01*x[kk,jj]+u01)+ exp(q11*x[kk,jj]*x[kk,jj]+a11*x[kk,jj]+u11)+ exp(q12*x[kk,jj]*x[kk,jj]+a12*x[kk,jj]+u12)+ exp(q02*x[kk,jj]*x[kk,jj]+a02*x[kk,jj]+u02)+ exp(q22*x[kk,jj]*x[kk,jj]+a22*x[kk,jj]+u22); p12=exp(q12*x[kk,jj]*x[kk,jj]+a12*x[kk,jj]+u12)/total; p22=exp(q22*x[kk,jj]*x[kk,jj]+a22*x[kk,jj]+u22)/total; rp0=rp0+.5*p2_12*p12/(.5*p2_12*p12+p2_22*p22); end; rp=rp0/3; if (r=rp) then missparent[kk]=2; notmissing[kk]=0; end; end; if ((obsparent[kk]=0) & (missparent[kk]=0)) then mfamily[kk]=1111; if (((obsparent[kk]=0) & (missparent[kk]=2)) | ((obsparent[kk]=2) & (missparent[kk]=0))) then mfamily[kk]=1122; if ((obsparent[kk]=2) & (missparent[kk]=2)) then mfamily[kk]=2222; if (((obsparent[kk]=0) & (missparent[kk]=1)) | ((obsparent[kk]=1) & (missparent[kk]=0))) then mfamily[kk]=1112; if ((obsparent[kk]=1) & (missparent[kk]=1)) then mfamily[kk]=1212; if (((obsparent[kk]=1) & (missparent[kk]=2)) | ((obsparent[kk]=2) & (missparent[kk]=1))) then mfamily[kk]=1222; if (flag_f[kk]=1) then do; missf[kk]=missparent[kk]; missm[kk]=obsparent[kk]; end; if (flag_m[kk]=1) then do; missm[kk]=missparent[kk]; missf[kk]=obsparent[kk]; end; end; end; return; FINISH ADDMISS; /**********************************************************************************/