!* Estimation of the size of a population with 3 Source Capture-Recapture !* including Goodness of fit-based confidence intervals !* written by Matthias an der Heiden, RKI Berlin, AnderHeidenM@rki.de program define recap version 9 capture drop Model qui gen Model = "" capture drop N_hat qui gen N_hat =. capture drop N_lb qui gen N_lb =. capture drop N_ub qui gen N_ub =. capture drop x_hat qui gen x_hat =. capture drop x_lb qui gen x_lb =. capture drop x_ub qui gen x_ub =. capture drop DoF qui gen DoF =. capture drop Loglikelihood qui gen Loglikelihood=. capture drop Gp2 qui gen Gp2 =. capture drop p_value qui gen p_value=. capture drop AIC qui gen AIC =. capture drop BIC qui gen BIC =. tempvar Gpower2 qui gen `Gpower2' =. local sig = invchi2(1,.95) syntax varlist(min=4 max=4) [if] [in] tokenize `varlist' local fre `1' local a `2' local b `3' local c `4' local vars `0' gsort -`a' -`b' -`c' qui replace `fre' = .1 if `fre'<.1 tempvar help qui gen `help' = sum(`fre') local N_obs = floor(`help'[7]) tempvar AB gen `AB' = `a' * `b' tempvar AC gen `AC' = `a' * `c' tempvar BC gen `BC' = `b' * `c' cons 1 `AB' = 0 cons 2 `AC' = 0 cons 3 `BC' = 0 local num=0 local i=1 while (`i' >= 0) { local j=2 while (`j' >= 0) { local k=3 while (`k' >= 0) { if (`i'==1) { local ab } else { local ab = 1 } if (`j'==2) { local ac } else { local ac = 2 } if (`k'==3) { local bc } else { local bc = 3 } // Berechnung der erwarteten Werte // ------------------------------- local ++num qui replace Model = "Model`i'`j'`k'" in `num' qui poisson `vars' `AB' `AC' `BC' `if' `in', const(`ab' `ac' `bc') qui replace Loglikelihood = e(ll) in `num' qui replace DoF = e(df_m) in `num' local deg = 6 - e(df_m) local Observations = e(N) tempvar pred nfreq hlp1 hlp2 qui predict double `pred' qui replace x_hat = floor(`pred'[8]) in `num' qui replace N_hat = `N_obs' + x_hat in `num' qui gen double `nfreq' = `fre' qui replace `nfreq' = `pred' in 8 qui gen double `hlp1' = 2*`nfreq' * log(`nfreq' / `pred') if `nfreq'>0 qui gen double `hlp2' = sum(`hlp1') qui replace `Gpower2' = `hlp2'[8] in `num' drop `hlp1' `hlp2' qui replace p_value = round(1 - chi2(`deg',`Gpower2'),.01) in `num' qui replace p_value = 1 if p_value==. in `num' // Bestimmung der Konfidenz-Intervalle // 1. Bestimmung der oberen Grenze // ------------------------------- tempvar tfreq qui gen double `tfreq' = `nfreq' if (`tfreq'[8] < 1) { qui replace `tfreq' = 1 in 8 } local Goodness = 0 local v=0 local z=0 while (abs(`Goodness' - `sig') > .01){ while (`Goodness' < `sig'){ qui replace `tfreq' = `tfreq'*(1+10^(-`z')) in 8 qui poisson `tfreq' `a' `b' `c' `AB' `AC' `BC' `if' `in', const(`ab' `ac' `bc') tempvar pred hlp1 hlp2 qui predict double `pred' qui gen double `hlp1' = 2*`tfreq' * log(`tfreq' / `pred') if `tfreq'>0 qui gen double `hlp2' = sum(`hlp1') local Goodness = `hlp2'[8] - `Gpower2'[`num'] drop `pred' `hlp1' `hlp2' local ++v } // di as result `v' as text " Iterationen und x= " as result `tfreq'[8] while (`Goodness' > `sig'){ qui replace `tfreq' = `tfreq'*(1-10^(-`z'-.8)) in 8 qui poisson `tfreq' `a' `b' `c' `AB' `AC' `BC' `if' `in', const(`ab' `ac' `bc') tempvar pred hlp1 hlp2 qui predict double `pred' qui gen double `hlp1' = 2*`tfreq' * log(`tfreq' / `pred') if `tfreq'>0 qui gen double `hlp2' = sum(`hlp1') local Goodness = `hlp2'[8]- `Gpower2'[`num'] drop `pred' `hlp1' `hlp2' local ++v } // di as result `v' as text " Iterationen und x= " as result `tfreq'[8] local ++z } qui replace x_ub = ceil(`tfreq'[8]) in `num' drop `tfreq' // di . // di as text "x_ub `i'`j'`k' = " as result x_ub[`num'] qui replace N_ub = `N_obs' + x_ub in `num' // 2. Bestimmung der unteren Grenze // -------------------------------- tempvar tfreq qui gen double `tfreq' = `nfreq' qui replace `tfreq' = 0 in 8 qui poisson `tfreq' `a' `b' `c' `AB' `AC' `BC' `if' `in', const(`ab' `ac' `bc') tempvar pred hlp1 hlp2 qui predict double `pred' qui gen double `hlp1' = 2*`tfreq' * log(`tfreq' / `pred') if `tfreq'>0 qui gen double `hlp2' = sum(`hlp1') local Goodness = `hlp2'[8] - `Gpower2'[`num'] if (`Goodness' > `sig' & `nfreq'[8] > 1) { qui replace `tfreq' = `nfreq' local Goodness = 0 local v=0 local z=0 while (abs(`Goodness' - `sig') > .01){ while (`Goodness' < `sig'){ qui replace `tfreq' = `tfreq'*(1-10^(-max(`z',log(2)/log(100)))) in 8 qui poisson `tfreq' `a' `b' `c' `AB' `AC' `BC' `if' `in', const(`ab' `ac' `bc') tempvar pred hlp1 hlp2 qui predict double `pred' qui gen double `hlp1' = 2*`tfreq' * log(`tfreq' / `pred') if `tfreq'>0 qui gen double `hlp2' = sum(`hlp1') local Goodness = `hlp2'[8] - `Gpower2'[`num'] // di as result `v' as text " Iterationen und x= " as result `tfreq'[8] as text " und Goodness = " /// // as res `Goodness' as text " wobei " as result `tpred'[8] drop `pred' `hlp1' `hlp2' local ++v } while (`Goodness' > `sig'){ qui replace `tfreq' = `tfreq'*(1+10^(-`z'-.8)) in 8 qui poisson `tfreq' `a' `b' `c' `AB' `AC' `BC' `if' `in', const(`ab' `ac' `bc') tempvar pred hlp1 hlp2 qui predict double `pred' qui gen double `hlp1' = 2*`tfreq' * log(`tfreq' / `pred') if `tfreq'>0 qui gen double `hlp2' = sum(`hlp1') local Goodness = `hlp2'[8]- `Gpower2'[`num'] drop `pred' `hlp1' `hlp2' local ++v } // di as result `v' as text " Iterationen und x= " as result `tfreq'[8] local ++z } } qui replace x_lb = floor(`tfreq'[8]) in `num' // di as text "x_lb `i'`j'`k' = " as result x_lb[`num'] qui replace N_lb = `N_obs' + x_lb in `num' di as text "." _c local k = `k'-3 } local j = `j'-2 } local i = `i'-1 } // Ergebnisse // qui replace sAIC= 2* (Loglikelihood+DoF) // qui replace sBIC= 2* Loglikelihood+DoF * log(`Observations') // qui replace BIC= round(sBIC[1]-sBIC,.01) // qui replace AIC= round(sAIC[1]-sAIC,.01) qui replace AIC= -round( 2*(-Loglikelihood[1]+DoF[1]) - 2* (-Loglikelihood+DoF),.01) qui replace BIC= -round(-2* Loglikelihood[1]+ DoF[1] * log(`Observations') /// + 2* Loglikelihood - DoF * log(`Observations'),.01) qui replace Gp2= round(`Gpower2',.01) qui replace p_value=1 in 1 gsort DoF -Model l Model DoF Gp2 p_value AIC BIC x_hat N_hat N_lb N_ub in 1/8, sep(0) noobs gsort -`a' -`b' -`c' end