*! version 1.0.2 June 2000				 (STB-56: sg148)
/* Mark Pearce								*/
/* Logistic regression - profile likelihood confidence intervals*/
/* Works with a logistic model from logistic or logit			*/
program define logprof, eclass
	version 6.0
	if "`e(cmd)'" == "" {
		di in red "Results of last estimate not found"
		exit 301
	}
	if e(cmd)!= "logistic" & e(cmd) != "logit" {
		di in red "Can only be used after a logistic model"
		error 301
	}
	local options "LEVel(int 95) NOGRaph OR INC(real 98989) From(real 98989) To(real 98989) INCFRAC(int -1) SAVing(str)"
	loc varlist "req ex min(1) max(1)"
	tempname vce names
	matrix `vce' = e(V)
	parse "`*'"
	if `"`saving'"' != "" {
		local saving `"saving(`saving')"'
	}
	parse "`varlist'", parse(" ")
	local varlblx : variable label `1'
	local xvar ="`1'"
	local yvar=e(depvar)
	local obll=e(ll)
	local obb=_b[`xvar']
	local seb=_se[`xvar']
	local fspec=1
	local tspec=1
	local ispec=1

	if `inc'==98989 {
		if `incfrac'!=-1 {
			local inc=`obb'/100
		}
		else local inc=0.01
	}

/*Get weights		*/
	local wtype=e(wtype)
	local wtexp=e(wexp)
	if "`e(wtype)'" == "" {
		if "`e(wexp)'" != "" {
			local weight "[`e(wtype)'`e(wexp)']"
		}
	}
	else local weight " "

/* Omit constant */
	local ncol = colsof(`vce')-1
	if `ncol'>1{
		matrix `names' = `vce'[1..`ncol',1]
		local all : rownames(`names')

/*Omit variable of interest */
		local exclude `1'
		tokenize `all'
		while "`1'" != "" {
        			if "`1'" != "`exclude'" {
                			local explvar "`explvar' `1'"
        			}
       			 mac shift
		}
	}
	else local explvar=" "
	preserve
	tempvar off   logl  off2 touse base voff

qui gen `touse'=e(sample)

/*Get offset	*/
local offset="`e(offset)'"
quietly {
	if "`e(offset)'" == "" {
		gen `off'=0
	}
	else gen `off'=`offset'

	tempname prof
	tempfile profile
	estimates hold `base'
	gen `off2'=0
	gen `voff'=0

	local levp=1-(`level'/100)
	local z=0.5*(invchi(1,`levp'))
	local oblow=`obb'-(`z'*`seb')
	local obupp=`obb'+(`z'*`seb')

	if `from'==98989 {
		local from=`oblow'
		local fspec=0
	}

	if `to'==98989 {
		local to=`obupp'
		local tspec=0
	}

/*Check range is valid	*/
/*Check from is low enough	*/
	replace `voff'=(`from'*`xvar')+`off'
	logistic `yvar' `explvar' if `touse' `weight', off(`voff')
	local dev2lik=e(ll)
	if `dev2lik'<`obll'-`z' {
		local beta=`from'
	}
	else {
		local j=1
		while `j'<=10 {
			replace `voff'=((`from'-`j')*`xvar')+`off'
			logistic `yvar' `explvar' if `touse' `weight',off(`voff')
			local dev2lik=e(ll)
			if `dev2lik'<`obll'-`z' {
				local beta=`from'-`j'
				local j=11
			}
			if `j'==10 { 
				local i=1
				while `i'<=1000 {
					replace `voff'=((`from'-(`i'*`seb'))*`xvar')+`off'
					logistic `yvar' `explvar' if `touse' `weight',off(`voff')
					local dev2lik=e(ll)
					if `dev2lik'<`obll'-`z' {
						local beta=`from'-(`i'*`seb')
						local i=1001
						local j=11	
					}
					else local i=`i'+1
				}
				if `i'==1000 {
					if `fspec'==1 {
						di in red "Range not Valid 
						exit 499
					}
					else {
			 di in red "Lower Confidence Limit not reached after 1000 iterations"
						exit 499
					}
				}
			}
			local j=`j'+1
		}
	}


/*Check to is high enough	*/
	replace `voff'=(`to'*`xvar')+`off'
	logistic `yvar' `explvar' if `touse' `weight', off(`voff')
	local dev2lik=e(ll)
	if `dev2lik'<`obll'-`z' {
		local to=`to'+`inc'
	}
	else {
		local j=1
		while `j'<=10 {
			replace `voff'=((`to'+`j')*`xvar')+`off'
			logistic `yvar' `explvar' if `touse' `weight',off(`voff')
			local dev2lik=e(ll)
			if `dev2lik'<`obll'-`z' {
				local to=`to'+`j'+`inc'
				local j=11
			}
			if `j'==10 { 
				local i=1
				while `i'<=1000 {
					replace `voff'=((`to'+(`i'*`seb'))*`xvar')+`off'
					logistic `yvar' `explvar' if `touse' `weight',off(`voff')
					local dev2lik=e(ll)
					if `dev2lik'<`obll'-`z' {
						local to=`to'+(`i'*`seb')
						local i=1001
						local j=11	
					}
					else local i=`i'+1
				}
				if `i'==1000 {
					if `tspec'==1 {
						di in red "Range not Valid 
						exit 499
					}
					else {
			 di in red "Upper Confidence Limit not reached after 1000 iterations"
						exit 499
					}
				}
			}
			local j=`j'+1
		}
	}
	postfile `prof' beta loglik using `profile', `replace' double

	while `beta'<= `to' {
		replace `off2'=(`beta'*`xvar')+`off'
		logistic `yvar' `explvar' if `touse' `weight', off(`off2')
		local dev2lik=e(ll)
		post `prof' (`beta') (`dev2lik')
		local beta=`beta'+`inc'
	}
	postclose `prof'
	use `profile', clear
	label data
	label var beta "Coefficient for `xvar'"
	label var loglik "Log-Likelihood"
}
	local cp=`obll'-`z'
	if "`nograph'"=="" {
graph loglik beta, xlab ylab t1(Profile log-likelihood function for `xvar') c(l) s(i) l2(Profile log-likelihood) yline(`cp') `saving'
	}
	estimates unhold `base'

quietly {
	drop if beta > `obb'
	summ beta if loglik==`cp'
	if r(N)>=1 {
		loc lower=(r(min)+r(max))/2
	}
	else {
		local obs=_N+1
		set obs `obs'
		replace beta=-999999 if beta==.
		replace loglik=`cp' if beta==-999999
		sort loglik
		tempvar where
		gen `where'=_n
		summ `where' if  beta==-999999
		local pos=r(min)
		summ beta if `where'==`pos'+1
		local beta1=r(min)
		summ loglik if beta==`beta1'
		local loglik1=r(min)
		 summ beta if `where'==`pos'-1
		local beta2=r(min)
		summ loglik if beta==`beta2'
		local loglik2=r(min)
		local lower = `beta1'-(`inc'*(`loglik1'-`cp')/(`loglik1'-`loglik2'))
	}
	use `profile', clear
	drop if beta < `obb'

	summ beta if loglik==`cp'
	if r(N)>=1 {
		loc upper=(r(min)+r(max))/2
	}
	else {
		local obs=_N+1
		set obs `obs'
		replace beta=-999999 if beta==.
		replace loglik=`cp' if beta==-999999
		sort loglik
		tempvar where
		gen `where'=_n
		summ `where' if beta==-999999
		local pos=r(min)
		summ beta if `where'==`pos'+1
		local beta1=r(min)
		summ loglik if beta==`beta1'
		local loglik1=r(min)
		summ beta if `where'==`pos'-1
		local beta2=r(min)
		summ loglik if beta==`beta2'
		local loglik2=r(min)
		local upper = `beta1'+(`inc'*(`loglik1'-`cp')/(`loglik1'-`loglik2'))
	}
}
	restore
	di _n in g "Maximum Profile Log-Likelihood = " `obll'
	di _n
	di in g _n _dup(79) "-"
	if "`or'"!="" {
		local betamax=exp(`obb')
		local lower=exp(`lower')
		local upper=exp(`upper')
di in g "`yvar'" _col(16) "|"  _col(20) "Odds Ratio" _col(40) "[" `level' "% Conf. Interval ]"
di in g _dup(15) "-" "+" _dup(63) "-"
di in g "`xvar'" _col(16) "|" in y _col(20)`betamax' _col(40) `lower' _col(52) `upper'
	}
else {
	di in g "`yvar'" _col(16) "|" _col(20) "Coef." _col(40) "[" `level' "% Conf. Interval]"
	di in g _dup(15) "-" "+" _dup(63) "-"
	di in g "`xvar'" _col(16) "|" in y _col(20)`obb' _col(40) `lower' _col(52) `upper'
}
di in g _dup(79) "-"
end