*! version 1.2.2 Yinghui Wei, Patrick Royston. 20 August 2020. program define ipdfc, sortpreserve /* Reconstruct IPD survival data from published KM survival or failure curve. Based on Guyot, Ades, Ouvwens and Welton (2012) BMC Methodology, 12:9. Translation to do-file by YW of their R program. Do-file converted to ado-file by PR. */ version 12.1 // varlist consists of surv, tstart, nrisk, trisk in that order syntax, surv(varname) TStart(varname) TRisk(varname) NRisk(varname) GENerate(string) SAVing(string) /// [ PROBability FAILure ISOtonic TOTevents(string)] if wordcount("`generate'") != 2 { di as err "generate() must provide names for two new variables" exit 198 } if "`totevents'" != "" confirm integer number `totevents' _prefix_saving `saving' local saving `"`s(filename)'"' local replace `"`s(replace)'"' if "`replace'" == "" { confirm new file `"`saving'"' } /* `generate' will create 2 variables known locally as `t_ipd' and `event_ipd' and save them to file `saving' */ local t_ipd : word 1 of `generate' local event_ipd : word 2 of `generate' confirm new var `t_ipd' confirm new var `event_ipd' confirm numeric variable `nrisk' confirm numeric variable `trisk' tempvar s n_censor n_hat cen d KM_hat last_i cen_t qui drop if `surv'==. & `tstart'==. & `trisk' ==. & `nrisk'==. quietly { // initialize vectors count if !missing(`nrisk') local n_int = r(N) /* Check that nrisk and trisk have consecutively non-missing values in obs 1, 2, ..., `n_int' and that the values in trisk are strictly increasing and that the values in nrisk are strictly non-increasing. */ forvalues i = 1/`n_int' { if missing(`trisk'[`i']) { di as err "`trisk' is missing for interval `i'" di as err "times at risk must be consecutively non-missing" error 198 } if missing(`nrisk'[`i']) { di as err "`nrisk' is missing for interval `i'" di as err "numbers at risk must be consecutively non-missing" exit 198 } local im1 = `i' - 1 if `i' > 1 { if `trisk'[`i'] <= `trisk'[`im1'] { di as err "times at risk must be in increasing order, number of participant at risk should be specified only once for each time point" di as err "times for intervals `im1' and `i' are " `trisk'[`im1'] " and " `trisk'[`i'] error 198 } if `nrisk'[`i'] > `nrisk'[`im1'] { di as err "numbers at risk must be in non-increasing order" di as err "numbers for intervals `im1' and `i' are " `nrisk'[`im1'] " and " `nrisk'[`i'] exit 198 } } else { if `trisk'[1] != 0 { di as err "the first time at risk must be 0" exit 198 } } } foreach thing in lower upper { tempvar `thing' gen long ``thing'' = 0 in 1/`n_int' local n_`thing' 0 } replace `upper' = _N in `n_int' forvalues i = 1/`n_int' { if `i'!=1 { count if `tstart' < `trisk'[`i'] replace `lower' = r(N) +1 in `i' } if `i'==1 { replace `lower' =1 in `i' } if (`i'!=`n_int'){ local pos = `i'+1 count if `tstart' < `trisk'[`pos'] replace `upper' = r(N) in `i' } else{ replace `upper' = _N in `n_int' } } count if !missing(`lower') in 1/`n_int' local n_lower = r(N) if `tstart'[_N] < `trisk'[`n_int'] { local tlast = `trisk'[`n_int'] di as err "Please add data points in surv(), tstart(), trisk(), nrisk() at time `tlast'" exit 198 } // !! Better to work with missing values in `s' via marksample etc. gen `s' = `surv' if ("`probability'" == "") replace `s' = `s' / 100 // convert data from percents if ("`failure'" == "failure") replace `s' = 1 - `s' // convert failure to survival prob replace `s' = 1 if `s'>1 & `s'!=. count if !missing(`s') local len2 = r(N) sort `tstart' if "`isotonic'" == "isotonic" { tempvar iso irax `s' `tstart', generate(`iso') reverse nograph //noi di "[isotonic regression performed to ensure survival probabilities are monotonic]" replace `s' = `iso' drop `iso' `iso'_p } else { // Simple way to make survival probabilities monotone decreasing forvalues i = 2/`len2'{ local pos = `i'-1 local temp = `s'[`pos'] if `s'[`i'] > `temp' { qui replace `s' = `temp' in `i' //noi di "Warning: survival probabilities must be monotone decreasing! Adjusted." } } } count if !missing(`lower') local n_int = r(N) local n_t = `upper'[`n_int'] local len = `nrisk'[1] if `len' > _N { tempvar to_drop gen byte `to_drop' = 0 set obs `len' replace `to_drop' = 1 if missing(`to_drop') } else { local tot_n_obs = _N+1 set obs `tot_n_obs' } local len2 `n_t' // `len2' = number of extracted data points gen long `n_censor' = 0 in 1/`len2' gen long `n_hat' = 0 in 1/`len2' gen long `cen' = 0 in 1/`len2' gen int `d' = 0 in 1/`len2' gen `KM_hat' = 1 in 1/`len2' gen `last_i' = 1 in 1/`len2' local sumdL 0 if `n_int' > 1 { local n_intervals = `n_int'-1 forvalues i = 1/`n_intervals' { /* Adjust tot. no. censored until n.hat = n.risk at start of interval (i+1) First approximation of no. censored on interval i */ local temp2 = `nrisk'[`i'] // di as error "Checking number at risk at `i' interval: `temp2'" local temp = `nrisk'[`i']*`s'[`lower'[`i'+1]]/`s'[`lower'[`i']]-`nrisk'[`i'+1] replace `n_censor' = round(`temp') in `i' if `n_censor'[`i']<0 replace `n_censor' = 0 in `i' while (`n_hat'[`lower'[`i'+1]] > `nrisk'[`i'+1]) | /// ((`n_hat'[`lower'[`i'+1]] < `nrisk'[`i'+1]) & (`n_censor'[`i']>=0)) { if `n_censor'[`i'] <= 0 { local start = `lower'[`i'] local end = `upper'[`i'] replace `cen' = 0 in `start'/`end' replace `n_censor' = 0 in `i' } if `n_censor'[`i'] > 0 & !missing(`n_censor'[`i']) { local len = `n_censor'[`i'] /*need "=" to specify the local variable correctly*/ di `len' gen `cen_t' = 0 in 1/`len' // distribute censored observations evenly over time forvalues j = 1/`len' { local temp =`tstart'[`lower'[`i']] + `j'*(`tstart'[`lower'[`i'+1]]-`tstart'[`lower'[`i']])/(`n_censor'[`i']+1) replace `cen_t' = `temp' in `j' } // Find no. censored on each time interval local start = `lower'[`i'] local end = `upper'[`i'] // one-step in R using hist(,plot=F)$count. not sure in stata forvalues k = `start'/`end' { qui count if `cen_t' >= `tstart'[`k'] & `cen_t' < `tstart'[`k'+1] replace `cen' = r(N) in `k' // di `k' _column(5) `cen'[`k'] } capture drop `cen_t' } // Find no. events and no. at risk on each interval to agree with KM estimates read from curves local pos = `lower'[`i'] replace `n_hat' = `nrisk'[`i'] in `pos' // gen `last' = `last_i'[`i'] !! Surely last does not need to be a variable, a local will do. local last = `last_i'[`i'] local start = `lower'[`i'] local end = `upper'[`i'] forvalues k = `start'/`end' { if `i'==1 & `k'==`start' { replace `d' = 0 in `k' replace `KM_hat' = 1 in `k' } else { replace `d' = round(`n_hat'[`k']*(1-(`s'[`k']/`KM_hat'[`last']))) in `k' replace `KM_hat' = `KM_hat'[`last']*(1-(`d'[`k']/`n_hat'[`k'])) in `k' } local index = `k' + 1 replace `n_hat' = `n_hat'[`k']-`d'[`k']-`cen'[`k'] in `index' if `d'[`k'] > 0 { local last `k' } if `d'[`k'] <0 { di as error "number of events is " `d'[`k'] " in interval `i'. It can't be less than zero!" } } local pos = `lower'[`i'+1] replace `n_censor' = `n_censor'[`i'] +(`n_hat'[`pos'] - `nrisk'[`i'+1]) in `i' } local pos = `lower'[`i'+1] if `n_hat'[`pos']< `nrisk'[`i'+1] { local index = `i' + 1 //di as text "`nrisk' was " `nrisk'[`index'] replace `nrisk' = `n_hat'[`pos'] in `index' //di as text "and now " `n_hat'[`pos'] } local index = `i' + 1 replace `last_i' = `last' in `index' di `last_i'[`index'] } } /* Time interval `n_int'. (last interval of risk table. Step 6 in Guyot et al 2012). YW note 03/June/2013. Guyot et al's method consider the case when the number at risk is not reported in the last interval. I modified it here such that if the total number at risk is not reported at the last interval, we do what they do. However, if the total number at risk is reported at the last interval, the computations are carried out using the same method as that for the other intervals and these are done in the previous loop. An option -fullrisktable- was given to note if the number at risk at the last interval is given; now deleted. */ if `n_int' >= 1 { if `n_int'==1 { local N_INT `n_int' replace `n_censor' = 0 in `N_INT' } if `n_int' > 1 { // Assume same censoring rate as average over previous time intervals local len = `n_int' - 1 qui sum `n_censor' in 1/`len' local sum_n_censor = r(sum) local temp1 = `sum_n_censor' * /// (`tstart'[`upper'[`n_int']]-`tstart'[`lower'[`n_int']])/(`tstart'[`upper'[`n_int'-1]] - `tstart'[`lower'[1]]) local pos `n_int' local temp2 = min(round(`temp1'),`nrisk'[`n_int']) replace `n_censor' = `temp2' in `pos' } /* If number at risk is not reported at intervals but only the total number of patients are available to each arm, in this case `n_int' = 1, and it appears that we only do the calculation in the last interval. */ local N_INT = `n_int' //line 253 - 259 added on 29/July/2015 if `n_censor'[`N_INT'] <= 0 { local start = `lower'[`N_INT'] local end = `upper'[`N_INT'] replace `cen' = 0 in `start'/`end' replace `n_censor' = 0 in `N_INT' local temp = `n_censor'[`N_INT'] } if `n_censor'[`N_INT'] > 0 { local len = `n_censor'[`N_INT'] // need "=" to specify the local variable correctly //di as text "number censored at the last interval is " `len' gen `cen_t' = 0 in 1/`len' // distribute censored observations evenly over time forvalues j = 1/`len'{ local temp =`tstart'[`lower'[`N_INT']] + `j'*(`tstart'[`upper'[`N_INT']]-`tstart'[`lower'[`N_INT']])/(`n_censor'[`N_INT']+1) replace `cen_t' = `temp' in `j' //di as text "censoring time is" `temp' } // Find no. censored on each time interval local start = `lower'[`N_INT'] local end = `upper'[`N_INT']-1 // one-step in R using hist(,plot=F)$count. not sure in stata forvalues k = `start'/`end'{ qui count if `cen_t' >=`tstart'[`k'] & `cen_t'<`tstart'[`k'+1] replace `cen' = r(N) in `k' //di `k' _column(5) `cen'[`k'] } capture drop `cen_t' } // Find no. events and no. at risk on each interval to agree with KM estimates read from curves local pos = `lower'[`N_INT'] replace `n_hat' = `nrisk'[`N_INT'] in `pos' local last = `last_i'[`N_INT'] local start = `lower'[`N_INT'] local end = `upper'[`N_INT'] forvalues k = `start'/`end' { if `KM_hat'[`last'] != 0 { replace `d' = round(`n_hat'[`k']*(1-(`s'[`k']/`KM_hat'[`last']))) in `k' } else { replace `d' = 0 in `k' } replace `KM_hat' = `KM_hat'[`last']*(1-(`d'[`k']/`n_hat'[`k'])) in `k' local index = `k' + 1 replace `n_hat' = `n_hat'[`k']-`d'[`k']-`cen'[`k'] in `index' *No. at risk cannot be negative if `n_hat'[`k'+1] < 0 { local index = `k' + 1 replace `n_hat' = 0 in `index' replace `cen' = `n_hat'[`k'] - `d'[`k'] in `k' } if `d'[`k'] !=0 { local last `k' } } /* Comment out the code relevant to totevents , 23 July2014. if "`totevents'" != "" will be expanded if remove comment. */ local N_INT = `n_int' local len = `n_censor'[`N_INT'] // if total no. of events reported, adjust no. censored so that total no. of events agrees. if "`totevents'" != "" { if `N_INT' > 1 { local endp = `upper'[`N_INT'-1] qui sum(`d') in 1/`endp' local sumdL = r(sum) // If total no. events already too big, then set events and censoring = 0 on all further time intervals if `sumdL' >= `totevents' { local startp = `lower'[`N_INT'] local endp = `upper'[`N_INT'] replace `d' = 0 in `startp'/`endp' local endp = `upper'[`N_INT'] -1 replace `cen' = 0 in `startp'/`endp' local startp = `lower'[`N_INT'] +1 local endp = `upper'[`N_INT'] +1 replace `n_hat' = `nrisk'[`N_INT'] in `startp'/`endp' // noi di "[Stop condition 1] Estimated tot events `sumdL' prior to the last interval is greater than reported tot events `totevents'" } } /* Otherwise adjust no. censored to give correct total no. events YW: note. for the ICON7 data, the following code in the condition if(`sumdL'<`tot_events') was not tested . */ if `sumdL' < `totevents' | `n_int'==1 { local endp = `upper'[`N_INT'] qui sum `d' in 1/`endp' local sumd = r(sum) local abc = `n_censor'[`N_INT'] //di as err "estimated tot events (include the last interval) is `sumd'; reported tot events `totevents'; estimated number censor `abc'" if `sumd' == `totevents' { // noi di "[Stop condition 2] Estimated tot events `sumd' is equal to the reported tot events `totevents'" } if `sumd' < `totevents' &`abc' == 0 { // noi di "[Stop condition 3] Estimated censor at the last interval is 0" } if ((`sumd'>`totevents') | (`sumd' < `totevents' & `abc' > 0)) { // noi di "redistributing the estimated censoring over the last interval" } while (`sumd'>`totevents') | (`sumd' < `totevents' & `abc' > 0) { noi di as text "estimated tot events is `sumd'; reported tot events `totevents';estimated number censor `abc'" di `abc' if `n_censor'[`N_INT'] <= 0 { local startp = `lower'[`N_INT'] local endp = `upper'[`N_INT'] replace `cen' = 0 in `startp'/`endp' replace `n_censor' = 0 in `N_INT' } if `n_censor'[`N_INT'] > 0 { local len = `n_censor'[`N_INT'] gen `cen_t' = 0 in 1/`len' forvalues j = 1/`len' { local temp =`tstart'[`lower'[`N_INT']] + `j'*(`tstart'[`upper'[`N_INT']]-`tstart'[`lower'[`N_INT']])/(`n_censor'[`N_INT']+1) replace `cen_t' = `temp' in `j' } local start = `lower'[`N_INT'] local end = `upper'[`N_INT']-1 forvalues k = `start'/`end'{ qui count if `cen_t' >=`tstart'[`k'] & `cen_t'<`tstart'[`k'+1] replace `cen' = r(N) in `k' //di `k' _column(5) `cen'[`k'] } capture drop `cen_t' } local index = `lower'[`N_INT'] replace `n_hat' = `nrisk'[`N_INT'] in `index' local last = `last_i'[`N_INT'] local start = `lower'[`N_INT'] local end = `upper'[`N_INT'] forvalues k = `start'/`end' { replace `d' = round(`n_hat'[`k']*(1-(`s'[`k']/`KM_hat'[`last']))) in `k' replace `KM_hat' = `KM_hat'[`last']*(1-(`d'[`k']/`n_hat'[`k'])) in `k' if `k' != `upper'[`n_int'] { local index = `k' + 1 replace `n_hat' = `n_hat'[`k'] - `d'[`k'] - `cen'[`k'] in `index' // No. at risk cannot be negative if `n_hat'[`k'+1] < 0 { replace `n_hat' = 0 in `index' replace `cen' = `n_hat'[`k'] - `d'[`k'] in `k' } } if `d'[`k'] !=0 { local last `k' } } local endp = `upper'[`N_INT'] qui sum `d' in 1/`endp' local sumd = r(sum) replace `n_censor' = `n_censor'[`N_INT'] + (`sumd' - `totevents') in `N_INT' local abc = `n_censor'[`N_INT'] if !((`sumd'>`totevents') | (`sumd' < `totevents' & `abc' > 0)) { // noi di "[Stop condition 4] Complete after redistributions of censor at the last interval!" } } } } else { // di as text "total number of events was not reported in the publication" } } /* Form IPD (create reconstructed survival data at patient level) */ // initialize vectors local endp = `nrisk'[1] gen `t_ipd' = `tstart'[`n_t'] in 1/`endp' gen byte `event_ipd' = 0 in 1/`endp' // write event time and event indicator (=1) for each event, as separate row in t_IPD and event_IPD local k 1 local N_t=`n_t' //di as error "Checkign N_t = `N_t'" forvalues j = 1/`N_t' { // di as error "check `j'" if `d'[`j'] != 0&`d'[`j']!=. { local startp `k' local endp = `k' + `d'[`j'] - 1 /*if `j'==75{ local temp2 = `d'[`j'] di as error "`temp2'" di as error "`startp'" di as error "`endp'" }*/ replace `t_ipd' = `tstart'[`j'] in `startp'/`endp' replace `event_ipd' = 1 in `startp'/`endp' local k = `startp' + `d'[`j'] } } // write censor time and event indicator (=0) for each censor, as separate row in t_IPD and event_IPd local N_t_1 = `n_t' - 1 forvalues j = 1/`N_t_1' { if `cen'[`j'] != 0 { local startp = `k' local endp = `k' + `cen'[`j'] - 1 local temp = (`tstart'[`j'] +`tstart'[`j'+1])/2 //di `temp' replace `t_ipd' = `temp' in `startp'/`endp' replace `event_ipd' = 0 in `startp'/`endp' local k = `startp' + `cen'[`j'] } } // Save results to file and tidy up preserve keep `t_ipd' `event_ipd' drop if missing(`t_ipd') save `"`saving'"', replace restore if ("`to_drop'" != "") drop if `to_drop' } end * v 1.0.0 PR/WvP 24may2013 /* Simplified version of WvP ira.ado, last updated 26/10/2012. Temporary variables: PART to define the current partition in merged sets. VIOL to detect a violation of the required ordering and the need for an additional merge. M mean of response variable. */ program define irax, sortpreserve version 11.0 quietly { syntax varlist(min=2 max=2 numeric) [if] [in] /// [, COMbine GENerate(string) noPTs REVerse ci noGRaph *] if "`generate'" != "" { gettoken generate rest : generate, parse(" ,") tokenize `"`rest'"', parse(" ,") while "`1'" != "" { if "`1'" == "replace" { local gen_replace replace continue, break } mac shift } } marksample touse gettoken y xvar : varlist sum `y' if `touse', meanonly if r(min) == r(max) { di as err "`y' does not vary" exit 198 } local mean_y = r(mean) count if (`touse'==1) & !inlist(`y', 0, 1) if (r(N)==0) local type LR // Logistic regression else local type LINR // Linear regression noi di _n as txt "type = `type'" tempvar PART VIOL HR M x if "`reverse'"=="" { gen `x' = -`xvar' if `touse' } else gen `x' = `xvar' if `touse' gen int `PART' = 1 if `touse'==1 sort `PART' `x' `y' replace `PART' = `PART'[_n-1] + (`x' > `x'[_n-1]) if !missing(`PART') & _n > 1 if "`reverse'"=="" { replace `x' = -`x' } sort `PART' by `PART': gen `M' = sum(`y') if !missing(`PART') by `PART': replace `M' = `M'[_N] / _N gen int `VIOL' = sum((`M' >= `M'[_n-1]) & (`PART'>`PART'[_n-1]) & !missing(`PART')) while `VIOL'[_N] != 0 { replace `PART' = `PART'[_n-1] + (`M' < `M'[_n-1]) if !missing(`PART') & _n > 1 sort `PART' by `PART': replace `M' = sum(`y') if !missing(`PART') by `PART': replace `M' = `M'[_N] / _N replace `VIOL' = sum((`M' >= `M'[_n-1]) & (`PART' > `PART'[_n-1]) & !missing(`PART')) } if ("`type'"=="LR") label variable `M' "Event probability" else if ("`type'"=="LINR") label variable `M' "Mean of `y'" if "`combine'" != "" { // For the -combine- option, need to create the "other" curve if ("`reverse'" == "") local Reverse reverse tempvar tempfit irax `y' `xvar' if `touse', `Reverse' nograph generate(`tempfit', replace) replace `M' = `M' + `tempfit' - `mean_y' drop `tempfit' `tempfit'_p // Update partition to one based on the new mean-y groups drop `PART' egen int `PART' = group(`M') if `touse' /* // Not clear if we should best re-estimate mean for new partition. Can easily be done: drop `M' sort `PART' by `PART': gen `M' = sum(`y') if !missing(`PART') by `PART': replace `M' = `M'[_N] / _N */ } if "`ci'" != "" { // Crude pointwise CI from partition model tempvar se lci uci tempname z scalar `z' = invnormal((100 + c(level))/200) if "`type'"=="LR" { local cmd logit local backtransf invlogit local transf logit } else { local cmd regress } `cmd' `y' ibn.`PART' if `touse' predict `se' if `touse', stdp gen `lci' = `backtransf'(`transf'(`M') - `z' * `se') gen `uci' = `backtransf'(`transf'(`M') + `z' * `se') drop `se' } /* Create graph. */ if "`graph'" != "nograph" { local title Isotonic Regression Analysis if ("`combine'" != "") local title `title' (extended) if "`type'" == "LINR" { local varlab: var lab `y' if ("`varlab'"=="" | index("`varlab'", ")") > 0) local varlab `y' local yplot = cond("`pts'" == "nopts", "", "`y'") if "`ci'" == "" { scatter `M' `yplot' `xvar' if `touse', sort connect(J .) msymbol(i oh) /// title(`title') ytitle("`varlab'") `options' } else { scatter `M' `lci' `uci' `yplot' `xvar' if `touse', sort connect(J J J .) lpattern(l shortdash ..) /// cmissing(y n n) msymbol(i i i oh) title(`title') ytitle("`varlab'") legend(off) `options' } } else { if "`ci'" == "" { line `M' `xvar' if `touse', sort connect(J) title(`title') `options' } else { line `M' `lci' `uci' `xvar' if `touse', sort connect(J J J) lpattern(l shortdash ..) /// cmissing(y n n) title(`title') legend(off) `options' } } } if "`generate'" != "" { local lab : var lab `M' compute `generate' = `M', `gen_replace' label("`lab'") noi di as txt "[variable `generate' created]" compute `generate'_p = `PART', `gen_replace' label("partition of `xvar'") noi di as txt "[variable `generate'_p created]" if "`ci'" != "" { compute `generate'_lci = `lci', `gen_replace' label("lower CL for `lab'") compute `generate'_uci = `uci', `gen_replace' label("upper CL for `lab'") noi di as txt "[variables `generate'_lci and `generate'_uci created]" } } } end * v 1.0.1 PR 03mar2013 program define compute version 8.0 gettoken type 0 : 0, parse("= ") bind // Process putative type local ok 0 foreach t in byte int long float double { if "`t'" == "`type'" { local ok 1 continue, break } } if !`ok' { if substr("`type'", 1, 3) == "str" { local strn = substr("`type'", 4, .) confirm integer number `strn' local ok 1 } } if !`ok' { local newvar `type' local type } else gettoken newvar 0 : 0, parse("= ") bind gettoken eqs 0 : 0, parse("= ") if "`eqs'" != "=" { di "{p}{err}syntax is {cmd:compute [{it:type}] {it:existing_var}|{it:newvar} = {it:exp}}" /// " [, {cmd:replace} {cmd:force} {cmd:label(}{it:label}{cmd:)}]{p_end}" exit 198 } gettoken Exp 0 : 0, parse(", ") bind syntax [if] [in] [, replace LABel(string) force ] if "`type'" != "" { local ok 0 foreach t in byte int long float double { if "`t'" == "`type'" { local ok 1 continue, break } } if !`ok' { if substr("`type'", 1, 3) == "str" { local strn = substr("`type'", 4, .) confirm integer number `strn' local ok 1 } } if !`ok' { di as err "type `type' not recognized" exit 198 } } capture confirm var `newvar', exact local rc = c(rc) if `rc' != 0 { // `newvar' does not exist; safe to create it from `Exp' generate `type' `newvar' = `Exp' `if' `in' } else { // `newvar' exists if "`replace'" != "" { // safe to recreate `newvar' replace `newvar' = `Exp' `if' `in' if "`type'" != "" { recast `type' `newvar', `force' } } else { di as err "`newvar' already defined" exit 110 } } if "`label'" != "" { label var `newvar' "`label'" } end exit History v 1.2.1 PR 31/July/2015. Tidy up code for lower and upper. Includes YW's recent changes in int() -> round(). v 1.1.0 YW 26/July/2015. Remove lower() and uppwer() as options as they are not needed. v 1.0.9 YW 25/July/2014. Output error message if the data information on the last interval of `trisk' are not given. Include messages on the conditions of completing reconstruction. v 1.0.8 YW 23/July/2014. Delete the option "fullrisktable" as it makes little difference. Comment out the relevant code. Keep the option "totevents". v 1.0.7 YW 23/July/2014. Minor update. Set the while condition in the redistribution of total events. v 1.0.7 YW 12/June/2013. Minor update. Set survival probability to be 1 if recorded as greater than 1. v 1.0.7 YW 03/June/2013. Add -fullrisktable- option, default ("`fullrisktable'" == "") means that the total number at risk is not given for the last interval. The reconstruction for last interval are following those specified in Guyot. But if the total number at risks is given at the last interval, the reconstruction for last interval are following those specified for the other intervals. v 1.0.6 PR 29may2013. Using irax not ira3 for isotonic regression. Minor tidying up. Renamed to reconstruct_ipd. Changed the -probability- option, default ("`probability'" == "") means percent. v 1.0.5 YW 24may2013 minor updates. correct a bracket matching error in "if notal no. of events reported...." v 1.0.4 PR/YW 24may2013 add isotonic regression as an alternative approach to correct survival/failure// probabilities if they are not monotonic v 1.0.3 YW 22may2013 modify the first while condition, to allow the computation when number of censor (n_censor) at interval i is equal to zero. prior to the first while condition, set the number of censor to be zero if negative. PR 22may2013 modify the code for lower() and upper() v 1.0.2 YW 20may2013 allow index of bigger interval (lower, upper) to be missing v 1.0.1 YW 17may2013 correct survival/failure probabilities if they are not monotonic v 1.0.0 YW/PR 16may2013