*! version 1.5.7  28nov2000 / 24sep2004
program define pkexamine_7, rclass byable(recall) sortpreserve
	version 7, missing

	syntax varlist(numeric min=2 max=2) [if] [in] /*
		*/ [, Trapezoid fit(integer 3) line log exp(real -1) noZero /*
	        */ Graph Connect(string) Symbol(string) T1title(string) * ]
		
	if "`graph'" == "" {
		syntax varlist(numeric min=2 max=2) [if] [in] /*
		*/ [, Trapezoid fit(integer 3) line log exp(real -1) noZero /*
	        */  Graph Connect(string) Symbol(string) T1title(string) ]
	}

	marksample touse
	tokenize `varlist'
	local time "`1'"
	local conc "`2'"

	tempvar lnc toest
	tempname auc_ln auc0_tm slope cmax auc_ex auc_lg /* 
		*/ lpt b0 b1 lastconc lc_cons

	if `fit' < 2 | `fit' > _N { 
		di as err /*
		*/ "fit() must be at least 2 and less than or equal to _N=" _N
		exit 198
	}

	qui count if `touse' 
	if `fit'>r(N) { 
		if "`if'"~="" {
			local myif=`"`if'"'
			gettoken myif2 myif: myif  
			local ifmsg="with`myif'"
		}

		di as err "fit() is greater than the number of obs. `ifmsg'"
		error cond(r(N)==0,2000,2001) 
	}

	if "`exp'" != "-1" { 
		qui su `time', meanonly
		capture assert `exp' >= r(max) | `exp' == 0
		if _rc { 
			di as err "exp must be 0 or greater than or equal to tmax= " r(max)
			exit 198
		}
	}
	else local exp = ""	

	if ("`log'"!="") + ("`exp'"!="") + ("`line'"!="") > 1 { 
		di as err /*
		*/ "only one of options log, exp, and line may be specified"
		exit 198
	}

	sort `touse' `time'
	capture by `touse' `time': assert _N==1 if `touse'
	if _rc { 
		di as err "time variable takes on repeated values"
		exit 459
	}
	capture assert `time'>=0 if `touse'
	if _rc { 
		di as err "time variable takes on negative values"
		exit 459
	}

	if "`zero'" == "" {	
		qui sort `touse' `time'
		capture assert `time'[1] == 0 if `touse'
		if _rc {
			di
			di as txt _col(5) /*
*/ "Warning: the point (0, 0) is not in your data.  It will be added."
			preserve
			local obs = _N + 1
			qui set obs `obs'
			qui replace `time' = 0 in l
			qui replace `conc' = 0 in l
			qui replace `touse' = 1 in l
		
		}
		capture assert `time' == `conc' if `time' == 0 & `touse'
		if _rc {
			di as err _col(5) "Warning: the point (0, 0) is not in your data. It will be added."
			preserve
			local obs = _N + 1
			qui set obs `obs'
			qui replace `time' = 0 in l
			qui replace `conc' = 0 in l
			qui replace `touse' = 1 in l
		}
		sort `touse' `time'
	}

				/* parsing ends */

	qui {
		scalar `lastconc' = `conc'[_N] 
		return scalar tmax = `time'[_N]
		gen byte `toest' = 0 
		replace `toest' = 1 in -`fit'/l


		gen double `lnc' = ln(`conc') if `touse'
		
		integ `conc' `time' if `touse', `trapezoid' 
		scalar `auc0_tm' = r(integral)

		regress `lnc' `time' if `toest'
		scalar `slope' = _b[`time']
		scalar `lc_cons' = _b[_cons]
		capture assert `slope' < 0
		if _rc {
			scalar `slope' = .
		}
		
		sum `conc' if `touse', meanonly
		scalar `cmax' = r(max)
		sum `time' if `touse' & `conc'==`cmax', meanonly
		return scalar tomc = r(max)

		MYLine `time' `conc' `toest' `touse'
		scalar `auc_ln' = `auc0_tm' + r(area)
		scalar `lpt' = r(l_point)

		MYExp `time' `conc' `toest' `touse'
		scalar `auc_ex' = `auc0_tm' + r(area)
		scalar `b0' = r(b0)
		scalar `b1' = r(b1)

		scalar `auc_lg' = `auc0_tm' + (`lastconc' / (-`slope'))
	}



	return scalar cmax = `cmax'
	/* return scalar tmax */ 
	return scalar ke = -`slope'
	return scalar half = ln(2) / -`slope'
	return scalar auc = `auc0_tm'
	return scalar auc_ln = `auc_lg'
	return scalar auc_exp = `auc_ex'
	return scalar auc_line = `auc_ln'



#delimit ;

	di ;
	di as txt _col(5) "                                      Maximum concentration = " 
		as res %9.0g return(cmax) ;
	di as txt _col(5) "                              Time of maximum concentration = " 
		as res %9.0g return(tomc) ;

	di as txt _col(5) "                            Time of last observation (Tmax) = "
		as res %9.0g return(tmax) ;
	di as txt _col(5) "                                           Elimination rate = " 
		as res %9.4f return(ke) ;
	di as txt _col(5) "                                                  Half life = " 
		as res %9.4f return(half) ; 

	di _n as txt _col(5)
	"Area under the curve" ;
	di as txt _col(5) /*
	*/ "{hline 16}{c TT}{hline 21}{c TT}{hline 15}{c TT}{hline 16}" ;
	di as txt _col(21) 
	"{c |}    AUC [0, inf.)    {c |} AUC [0, inf.) {c |}  AUC [0, inf.)" ;

	di as txt _col(5) 
	" AUC [0, Tmax]  {c |} Linear of log conc. {c |}   Linear fit  {c |} Exponential fit" ;

	di as txt _col(5) /*
	*/ "{hline 16}{c +}{hline 21}{c +}{hline 15}{c +}{hline 16}" ;

	di _col(6) as res %9.2f return(auc) 
			as txt "      {c |}      " as res %9.3f return(auc_ln)
			as txt "      {c |}  " as res %9.3f return(auc_line)
			as txt "    {c |}   " as res %9.3f return(auc_exp) ;
	di as txt _col(5) /*
	*/ "{hline 16}{c BT}{hline 21}{c BT}{hline 15}{c BT}{hline 16}" ;
	di as txt	_col(5) "Fit based on last " as res `fit' " {gr}points." ;

#delimit cr

	if "`graph'"=="" { 
		exit
	}

	local tlab: variable label `time'
	if "`tlab'" == "" {
		local tlab "Analysis Time"
	}
	label var `time' "`tlab'"

	local clab: variable label `conc'
	if "`clab'" == "" {
		local clab "Concentration"
	}
	label var `conc' "`clab'"

	if "`log'" != "" {
		local lnlab "Log Concentration"
	}
	label var `lnc' "`lnlab'"

	if "`connect'" != "" {
		local conopt "connect(`connect') sort"
	}
	else 	local conopt "connect(ls) sort"

	if "`symbol'" != "" {
		local symopt "symbol(`symbol') "
	}
	else 	local symopt "symbol(o.)"

	if "`t1title'" != "" {
		local t1opt "t1title(`t1title')"
	}
	else 	local t1opt   "t1title(" ")"

	
	if "`log'" == "" & "`line'" == "" & "`exp'" == "" {
		gr7 `conc' `time', `t1opt' `conopt' `options' `symopt', if `touse'
	}

	else if "`log'" != "" {
		if `auc_lg'>=. {
			noi di as err /*
		*/ "unable to compute area under extension, graph not produced"
			exit 459
		}
		capture preserve
		local obs = _N+1
		qui set obs `obs'
		qui replace `lnc' = 0 if _n == _N
		qui replace `time' = -`lc_cons' / `slope' if _n == _N
		gr7 `lnc' `time', `conopt' `options' `symopt', if `touse'
	}

	else if "`line'" != "" {
		if `auc_ln'>=. {
			noi di as err /*
		*/ "unable to compute area under extension, graph not produced"
			exit 459
		}
		qui {
			capture preserve
			local obs = _N+1
			set obs `obs'
			replace `conc' = 0 if _n == _N
			replace `time' = `lpt' if _n == _N
			gr7 `conc' `time', `conopt' `options' `symopt', if `touse'
		}
	}


	else if "`exp'" != "" {
		if `auc_ex'>=. {
			noi di as err /*
		*/ "unable to compute area under extension, graph not produced"
			exit 459
		}
		qui {
			if "`exp'" == "0" {
				local exp = `lpt'
			}

			capture preserve
			local org = _N
			local obs = _N+10
			set obs `obs'
			local inc = (`exp' - return(tmax)) / 10

			replace `time' = `time'[_n-1] + `inc' if _n > `org'

			tempvar conc1
			gen double `conc1'  = exp((`b0' + `b1' * `time')) if _n > `org'
			replace `conc1' = `conc' if _n == `org'

			gr7 `conc' `conc1' `time', `conopt' `options' `symopt' `t1opt', if `touse'
		}
	}

	else {
		di as err "Impossible"
		exit 198
	}
end


program define MYLine, rclass
	args time conc toest touse
	tempvar time1 conc1

	regress `conc' `time' if `toest'
	capture assert _b[`time'] < 0
	if _rc {
		return scalar area = .
		exit
	}
	return scalar l_point = -_b[_cons]/_b[`time']
	sort `touse' `time'
	gen double `time1' = `time'[_N] in 1
	replace `time1' = return(l_point) in 2
	gen double `conc1' = `conc'[_N] in 1
	replace `conc1' = 0 in 2
	return scalar area = .5 * `conc1'[1] * (`time1'[2] - `time1'[1])
end


program define MYExp, rclass
	args time conc toest touse

	ereg `conc' `time' if `toest'
	capture assert _b[`time'] < 0
	if _rc {
		return scalar area = .
		exit
	}
	ret scalar b0 = _b[_cons]
	ret scalar b1 = _b[`time']
	tempvar time1
	sort `touse' `time' /* in case ereg changed sort order */
	gen double `time1' = `time'[_N] in 1
	ret scalar area = (-1 / _b[`time']) * exp((_b[_cons] + _b[`time']*`time1'[1]))
end
