//We used Stata 14.1 on a Windows 10 PC import excel "\\Client\C$\Abortion\NSFGExtract.xls", sheet("Sheet1") firstrow clear //MLE procedure program drop urzip1 program urzip1 version 14 args lnf theta1 theta2 theta3 tempvar lambda psi quietly gen double `lambda' = exp(`theta1')*invlogit(`theta2') quietly gen double `psi' = invlogit(`theta3') quietly replace `lnf' = ln(1-`psi')-`lambda'+$ML_y1*ln(`lambda')-lnfactorial($ML_y1) if $ML_y1>0 quietly replace `lnf' = ln(`psi'+(1-`psi')*(exp(-`lambda'))) if $ML_y1==0 end //Baseline model ml model lf urzip1 (Count:abcnt = age1519 age2024 age2529 age3034 parity econt hispanic black raceotr educat educ13 catholic protestant relotr nvrcont prghlp infserv homo bisex timescoh fmarno nvrwork,exposure(sexposure))(Reporting:=age1519 age2024 age2529 age3034 hispanic black raceotr educat educ13 catholic protestant relotr married pmarried livtogetr inc100 inc200 inc300 inc400 inc500)(Positive:age1519 age2024 age2529 age3034 hispanic black raceotr educat educ13 catholic protestant relotr nvrcont prghlp infserv homo bisex timescoh fmarno nvrwork born80 born90) if sexposure>0 [aw=wgt], vce(robust) ml max //Reporting rates gen edu11=educat<12 gen edu12=educat==12 mata Z=st_data( ., ("age1519", "age2024", "age2529", "age3034", "hispanic","black", "raceotr", "educat", "educ13", "catholic", "protestant", "relotr", "married", "pmarried", "livtogetr", "inc100","inc200","inc300","inc400","inc500")) G=st_data( ., ("age1519", "age2024", "age2529", "age3034", "age3545", "white","black", "raceotr","hispanic", "edu11", "edu12","educ13","relnone", "catholic", "protestant", "relotr", "nmarried", "married", "pmarried", "livtogetr","inc00", "inc100","inc200","inc300","inc400","inc500")),J(15420,1,1) b=st_matrix("e(b)")' gamahat=b[24..43] gamaC=b[44]*J(15420,1,1) Zgama=Z*gamahat+gamaC pihat=invlogit(Zgama) mean(pihat) sumGpihat=colsum(pihat:*G) sumG=colsum(G) meanpi=sumGpihat:/sumG meanpi V=st_matrix("e(V)")' Vhat=V[24..44,24..44] for (i=1; i<=27; i++){ Gi=G[.,i] ZG=Z:*Gi Ggama=((Z'*(pihat:*(J(15420,1,1)-pihat)))\(Gi'*(pihat:*(J(15420,1,1)-pihat))))/15420 Vpi=Ggama'*Vhat*Ggama sqrt(Vpi) } end //Finding the alternative maximum mat init=e(b) mat a=init[1,1..23] mat b=init[1,24..44] mat c=init[1,45..67] mat bm=-b mat badd=b[1,1..4],0,0,b[1,5..12],0,0,0,0,0,0,0,0,b[1,21] mat ab=badd+a ml model lf urzip1 (Count:abcnt = age1519 age2024 age2529 age3034 parity econt hispanic black raceotr educat educ13 catholic protestant relotr nvrcont prghlp infserv homo bisex timescoh fmarno nvrwork,exposure(sexposure))(Reporting:=age1519 age2024 age2529 age3034 hispanic black raceotr educat educ13 catholic protestant relotr married pmarried livtogetr inc100 inc200 inc300 inc400 inc500)(Positive:age1519 age2024 age2529 age3034 hispanic black raceotr educat educ13 catholic protestant relotr nvrcont prghlp infserv homo bisex timescoh fmarno nvrwork born80 born90) if sexposure>0 [aw=wgt], vce(robust) ml init ab bm c ml max