*! survlsl v1.2.0 Long Hong 15Oct2017 * This version only contains left censoring and left truncation * Arguments: income, threshold, model, censor/trunc, percentage (only for censor) capture program drop survlsl program define survlsl, rclass version 12 syntax varlist(max=1 numeric), THREShold(real) CENSORPCT(real) MODEL(string asis) marksample touse qui count if `touse' if `r(N)' == 0 { error 2000 } *** Assert capture assert `varlist' > 0 // Varlist if c(rc) { display in red "Observations must be > 0" exit 9 } capture assert `varlist' != . if c(rc) { display in red "No missing observation is allowed" exit 9 } capture assert `threshold' > 0 // Threshold if c(rc) { display in red "Observations must be > 0" exit 9 } qui sum `varlist' local min = r(min) capture assert `threshold' <= `min' if c(rc) { display in red "Threshold must be no larger than the smallest obversation" exit 9 } capture assert "`model'" == "lognormal" | "`model'" == "weibull" | "`model'" == "loglogistic" // Models if c(rc) { display in red "Please select a proper model using the exact names:" display in blue "(1) lognormal (2) weibull (3) loglogistic " exit 9 } capture assert `censorpct' >= 0 & `censorpct' < 1 // Censor Percent if c(rc) { display in red "Percentage of censoring must be in (0, 1)" display in blue "Value 0 implies 'truncated data' " exit 9 } *** Programme tempname k censor if `censorpct' == 0 { qui gen `k' = `threshold' qui ml model lf `model'_lefttrunc (alpha: `varlist' = ) (beta: `k' = ) ml maximize display " " display in yellow "Left Truncated Model" } else { preserve local newobs = round(_N/(1-`censorpct'), 1) qui set obs `newobs' qui gen `censor' = (`varlist' != .) qui replace `varlist' = `threshold' if `varlist' == . qui ml model lf `model'_leftcensor (alpha: `varlist' = ) (beta: `censor' = ) ml maximize restore display " " display in yellow "Left Censored Model" } *** Gini local alpha_mle = [alpha]_cons local beta_mle = [beta]_cons local gini_lognormal "2*normal(`beta_mle'/(sqrt(2)))-1" local gini_weibull "1 - 2^(-1/`beta_mle')" local gini_loglogistic "1/`beta_mle'" *** Confidence Intervals tempname table mat `table' = r(table) local beta_sd = `table'[2,2] // S.D. for beta local lognormal_fd "sqrt(2)*normalden(`beta_mle'/(sqrt(2)))" local weibull_fd "`beta_mle'^(-2) * 2^(-(1/`beta_mle')) * ln(2)" local loglogistic_fd "-`beta_mle'^(-2)" local IC1_low = `gini_`model'' - invnormal(0.975) * ``model'_fd' * `beta_sd' local IC1_high = `gini_`model'' + invnormal(0.975) * ``model'_fd' * `beta_sd' local IC2_high_lognormal = 2*normal((`beta_mle' + invnormal(0.975)*`beta_sd')/sqrt(2)) - 1 local IC2_low_lognormal = 2*normal((`beta_mle' - invnormal(0.975)*`beta_sd')/sqrt(2)) - 1 local IC2_low_weibull = 1 - 2^(-(1/(`beta_mle' - invnormal(0.975)*`beta_sd'))) local IC2_high_weibull = 1 - 2^(-(1/(`beta_mle' + invnormal(0.975)*`beta_sd'))) local IC2_low_loglogistic = 1/(`beta_mle' - invnormal(0.975)*`beta_sd') local IC2_high_loglogistic = 1/(`beta_mle' + invnormal(0.975)*`beta_sd') tempname ic mat `ic' = (round(`IC1_low', 0.00001), round(`IC1_high', 0.00001) \ round(`IC2_low_`model'', 0.00001), round(`IC2_high_`model'', 0.00001)) matrix rownames `ic' = "Conf Interval 1" "Conf Interval 2" matrix colnames `ic' = "Lower" "Upper" *** Display display " " display in blue "Estimated Parameters: " if "`model'" == "lognormal" { display in blue " MLE location = " round(`alpha_mle', 0.00001) display in blue " MLE scale = " round(`beta_mle', 0.00001) } else { display in blue " MLE scale = " round(`alpha_mle', 0.00001) display in blue " MLE shape = " round(`beta_mle', 0.00001) } display "" display in yellow "Parametric Gini = " round(`gini_`model'', 0.00001) display " " display in blue "Parametric Gini 95% Confidence Interval: " display in blue "C.I. 1 is derived from the delta method;" display in blue "C.I. 2 is derived from a direct approach." matlist `ic', twidth(15) border(rows) rowtitle("") *** Saved Estimates tempname b V return scalar beta = [beta]_cons return scalar alpha = [alpha]_cons return scalar gini = `gini_`model'' mat `b' = (`alpha_mle', `beta_mle') matrix rownames `b' = "MLE" matrix colnames `b' = alpha beta local alpha_sd = `table'[2,1] local beta_sd = `table'[2,2] mat `V' = (`alpha_sd', `beta_sd') matrix rownames `V' = Std_Dev matrix colnames `V' = alpha beta return matrix conf_interval = `ic' return matrix variances = `V' return matrix estimates = `b' end