*! traj - May 17, 2013 Bobby L Jones
**********************************************************************************
program traj, eclass
version 9
syntax [if] [in], Model1(str) Var1(varlist numeric) Indep1(varlist numeric) ///
Order1(numlist int >=-1 <=5) ///
[ Min1(real 0.0) Max1(real 1e-38) Iorder1(numlist int >=-1 <=5) ///
Risk1(varlist numeric) Weight1(varlist numeric) Expos1(varlist numeric) ///
Tcov1(varlist numeric) Plottcov1(string) Refgroup1(int 1) ///
Dropout1(numlist int >=-1 <=2) Dcov1(varlist numeric) Obsmar1(varlist numeric) ///
Outcome1(varlist numeric) Omodel1(str) Rorder1(numlist int >=-1 <=3) ///
Model2(str) Var2(varlist numeric) Indep2(varlist numeric) ///
Order2(numlist int >=-1 <=5) Min2(real 0.0) Max2(real 1e-38) ///
Iorder2(numlist int >=-1 <=5) Risk2(varlist numeric) Expos2(varlist numeric) ///
Tcov2(varlist numeric) Refgroup2(int 1) Dropout2(numlist int >=-1 <=2) ///
Model3(str) Var3(varlist numeric) Indep3(varlist numeric) ///
Order3(numlist int >=-1 <=5) Min3(real 0.0) Max3(real 1e-38) ///
Iorder3(numlist int >=-1 <=5) Risk3(varlist numeric) Expos3(varlist numeric) ///
Tcov3(varlist numeric) Dropout3(numlist int >=-1 <=2) ///
Multgroups(integer 0) Start(string) DETAIL TRACE NOVAR OUTOFSAMPLEPROBS]
marksample touse
quietly count if `touse'
local maxmodels = 3
if `r(N)' == 0 {
error 2000
}
local rN = `r(N)'
forval i = 2/`maxmodels' {
if ( length( "`model`i''" ) == 0 ) {
local model`i' = "none"
}
}
local itdetail = 0
if "`detail'" == "detail" {
local itdetail = 1
}
local ittrace = 0
if "`trace'" == "trace" {
local ittrace = 1
}
local itnovar = 0
if "`novar'" == "novar" {
local itnovar = 1
}
local oosprobs = 0
if "`outofsampleprobs'" == "outofsampleprobs" {
local oosprobs = 1
}
local traj_start_len = 0
if "`start'" != "" {
cap confirm matrix `start'
if _rc {
di as err "Error: start must specify a matrix"
exit 198
}
tempname traj_start
local traj_start_len = colsof(`start')
matrix traj_start = `start'
}
local plottcov1_len = 0
if "`plottcov1'"!="" {
cap confirm matrix `plottcov1'
if _rc {
di as err "Error: plottcov1 must specify a matrix"
exit 198
}
tempname traj_plottcov1
local plottcov1_len = colsof(`plottcov1')
matrix traj_plottcov1 = `plottcov1'
}
forval i = 1/`maxmodels' {
if !("`model`i''" == "none") {
if ( "`model`i''" != "cnorm" & "`model`i''" != "zip" &"`model`i''" != "logit" ) {
di as err "Error: model`i' is not cnorm, zip, or logit."
exit 198
}
if "`model`i''"=="cnorm" {
if `max`i'' == 1e-38 {
di as err "Error: max`i' required for model`i' cnorm."
exit 198
}
tokenize "`var`i''"
while "`1'" != "" {
local vr `1'
cap sum `vr'
cap assert r(min) >= `min`i'' & r(max) <= `max`i''
if _rc {
di as err "`vr' is not within min`i' = `min`i'' and max`i' = `max`i''."
exit 198
}
mac shift
}
}
if "`model`i''"=="logit" {
tokenize "`var`i''"
while "`1'" != "" {
local vr `1'
cap assert `vr' == 0 | `vr' == 1 | `vr' >= .
if _rc {
di as err "Error: `vr' is not 0/1 (logit model)."
exit 198
}
mac shift
}
}
if "`model`i''"=="zip" {
tokenize "`var`i''"
while "`1'" != "" {
local vr `1'
cap assert `vr' >= 0
if _rc {
di as err "Error: `vr' has negative values (zip model)."
exit 198
}
mac shift
}
}
}
local numindep`i' : word count `indep`i''
local numvar`i' : word count `var`i''
if ( `numindep`i'' != `numvar`i'' ) {
di as err "Error: number of variables in var`i' and indep`i' must match."
exit 198
}
local numrisk`i' : word count `risk`i''
local numexpos`i' : word count `expos`i''
local numtcov`i' : word count `tcov`i''
local numdcov`i' : word count `dcov`i''
local numoutcome`i' : word count `outcome`i''
local numweight`i' : word count `weight`i''
local numobsmar`i' : word count `obsmar`i''
if ( `numoutcome`i'' > 0 ) {
if ( "`omodel`i''" != "normal" & "`omodel`i''" != "cnorm" & ///
"`omodel`i''" != "poisson" & ///
"`omodel`i''" != "logit" & "`omodel`i''" != "none" ) {
di as err "Error: omodel`i' is not normal, cnorm, poisson, or logit."
exit 198
}
/*
if "`omodel`i''"=="cnorm" {
if `max`i'' == 1e-38 {
di as err "Error: max`i' required for model`i' cnorm."
exit 198
}
cap sum `outcome`i''
cap assert r(min) >= `omin`i'' & r(max) <= `omax`i''
if _rc {
di as err "`outcome`i'' is not within omin`i' = `omin`i'' and omax`i' = `omax`i''."
exit 198
}
}
*/
if "`omodel`i''"=="logit" {
cap assert `outcome`i'' == 0 | `outcome`i'' == 1 | `outcome`i'' >= .
if _rc {
di as err "Error: `outcome`i'' is not 0/1 (logit model)."
exit 198
}
}
if "`omodel`i''"=="poisson" {
cap assert `outcome`i'' >= 0
if _rc {
di as err "Error: `outcome`i'' has negative values (poisson model)."
exit 198
}
}
}
if ( `numexpos`i'' > 1 ) {
tokenize "`expos`i''"
while "`1'" != "" {
local vr `1'
cap assert `vr' > 0
if _rc {
di as err "Error: exposure variable `vr' is not > 0."
exit 198
}
mac shift
}
if ( `numexpos`i'' != `numvar`i'' ) {
di as err "Error: number of variables in var`i' and expos`i' must match."
exit 198
}
if ( "`model`i''" != "zip" ) {
di as err "Error: exposure is only allowed with the zip model."
exit 198
}
}
if ( `numoutcome`i'' > 1 ) {
di as err "Error: only one variable is allowed for outcome`i'."
exit 198
}
if ( `numobsmar`i'' > 1 ) {
di as err "Error: only one variable is allowed for obsmar`i'."
exit 198
}
if ( `numtcov`i'' != 0 & `numtcov`i'' / `numindep`i'' != int( `numtcov`i'' / `numindep`i'' ) ) {
di as err "Error: number of variables in tcov`i' and var`i' must match or be a multiple."
exit 198
}
if ( `numdcov`i'' != 0 & `numdcov`i'' / `numindep`i'' != int( `numdcov`i'' / `numindep`i'' ) ) {
di as err "Error: number of variables in dcov`i' and var`i' must match or be a multiple."
exit 198
}
if (length("`omodel`i''")==0) {
local omodel`i'="none"
}
}
if ( `plottcov1_len' != 0 & `plottcov1_len' / `numindep1' != int( `plottcov1_len' / `numindep1' ) ) {
di as err "Error: there should be `numtcov1' values for plottcov1."
exit 198
}
if ( `numweight1' > 1 ) {
di as err "Error: only one variable is allowed for weight."
exit 198
}
forval i = 1/`maxmodels' {
local numorder`i' 0
tokenize "`order`i''"
while "`1'" != "" {
local numorder`i' = `numorder`i'' + 1
local c`numorder`i'' `1'
local ordr`i' "`ordr`i'' "`1'""
mac shift
}
if ( `numorder`i'' != 0 ) {
matrix traj_mat_order`i' = J( 1, `numorder`i'', 0 )
tokenize "`order`i''"
forval j = 1/`numorder`i'' {
matrix traj_mat_order`i'[1,`j'] = real("`1'")
mac shift
}
}
local numiorder`i' 0
tokenize "`iorder`i''"
while "`1'" != "" {
local numiorder`i' = `numiorder`i'' + 1
local c`numiorder`i'' `1'
local iordr`i' "`iordr`i'' "`1'""
mac shift
}
if ( `numiorder`i'' != 0 ) {
matrix traj_mat_iorder`i' = J( 1, `numiorder`i'', 0 )
tokenize "`iorder`i''"
forval j = 1/`numiorder`i'' {
matrix traj_mat_iorder`i'[1,`j'] = real("`1'")
mac shift
}
}
if ( "`model`i''" != "zip" & `numiorder`i'' != 0 ) {
di as err "Error: iorder is only allowed with the zip model."
exit 198
}
if ( `numiorder`i'' != `numorder`i'' & `numiorder`i'' != 1 & `numiorder`i'' != 0 ) {
di as err "Error: there should be 1 or `numorder`i'' iorders for iorder`i'."
exit 198
}
}
if ( `numrisk1' > 0 ) {
if ( `numorder1' < 2 ) {
di as err "Error: risk is only valid when the number of groups is > 1."
exit 198
}
}
if ( `refgroup1' > 1 ) {
if ( `numrisk1' == 0 ) {
di as err "Error: refgroup is only valid when the risk is specified."
exit 198
}
if ( `refgroup1' > `numorder1' ) {
di as err "Error: refgroup exceeds the number of groups."
exit 198
}
}
if ( `refgroup2' > 1 ) {
if ( `numrisk2' == 0 ) {
di as err "Error: refgroup2 is only valid when the risk2 is specified."
exit 198
}
if ( `refgroup2' > `numorder2' ) {
di as err "Error: refgroup2 exceeds the number of groups."
exit 198
}
}
local numrorder1 0
tokenize "`rorder1'"
while "`1'" != "" {
local numrorder1 = `numrorder1' + 1
local c`numrorder1' `1'
local rordr1 "`rordr1' "`1'""
mac shift
}
if ( `numrorder1' != 0 ) {
matrix traj_mat_rorder1 = J( 1, `numrorder1', 0 )
tokenize "`rorder1'"
forval i = 1/`numrorder1' {
matrix traj_mat_rorder1[1,`i'] = real("`1'")
mac shift
}
}
/*
if ( `numrorder1' > 0 ) {
di as err "rorder is experimental and may be available. Contact bjones@andrew.cmu.edu"
exit 198
}
*/
if ( "`model1'" != "cnorm" & `numrorder1' != 0 ) {
di as err "Error: rorder is only supported for the cnorm model."
exit 198
}
local numdropout1 0
tokenize "`dropout1'"
while "`1'" != "" {
local numdropout1 = `numdropout1' + 1
local c`numdropout1' `1'
local drpout1 "`drpout1' "`1'""
mac shift
}
if ( `numdropout1' != 0 ) {
matrix traj_mat_dropoutorder1 = J( 1, `numdropout1', 0 )
tokenize "`dropout1'"
forval i = 1/`numdropout1' {
matrix traj_mat_dropoutorder1[1,`i'] = real("`1'")
mac shift
}
}
local numdropout2 0
local numdropout3 0
cap qui drop _traj_Group
qui gen _traj_Group = .
cap qui drop _traj_Prob*
local out1 "_traj_Group"
forvalues j = 1/`numorder1' {
qui gen _traj_ProbG`j' = .
local out1 "`out1' _traj_ProbG`j'"
}
local outMat `out1'
if ( `numorder2' != 0 & `multgroups' == 0 ) {
cap qui drop _traj_Model2_Group
qui gen _traj_Model2_Group = .
local out2 "_traj_Model2_Group"
forvalues j = 1/`numorder2' {
cap qui drop _traj_Model2_ProbG`j'
qui gen _traj_Model2_ProbG`j' = .
local out2 "`out2' _traj_Model2_ProbG`j'"
}
local outMat "`outMat' `out2'"
}
/*
if ( `multgroups' == 0 ) {
forval i = 2/`maxmodels' {
if `numorder`i'' != 0 {
cap qui drop _traj_Model`i'_Group
qui gen _traj_Model`i'_Group = .
local out`i' "_traj_Model`i'_Group"
forvalues j = 1/`numorder`i'' {
cap qui drop _traj_Model`i'_ProbG`j'
qui gen _traj_Model`i'_ProbG`j' = .
local out`i' "`out`i'' _traj_Model`i'_ProbG`j'"
}
local outMat "`outMat' `out`i''"
}
}
}
*/
forval i = 1/`maxmodels' {
if `numorder`i'' != 0 {
local plotVarLabels`i' "trajT"
forvalues j = 1/`numorder`i'' {
local plotVarLabels`i' "`plotVarLabels`i'' Avg`j'"
}
forvalues j = 1/`numorder`i'' {
local plotVarLabels`i' "`plotVarLabels`i'' Est`j'"
}
forvalues j = 1/`numorder`i'' {
local plotVarLabels`i' "`plotVarLabels`i'' L95`j'"
local plotVarLabels`i' "`plotVarLabels`i'' U95`j'"
}
}
tokenize "`plotVarLabels`i''"
}
if `numdropout1' != 0 {
if ( `numdropout1' != 1 & `numdropout1' != `numorder1' ) {
di as err "Error: there should be 1 or `numorder`i'' values for dropout."
exit 198
}
forvalues j = 1/`numdropout1' {
local plotVarLabels1 "`plotVarLabels1' Dropout`j'"
}
tokenize "`plotVarLabels1'"
}
tempname b V parmData varData
matrix parmData = J( 1, 300, 0 )
matrix varData = J( 300, 300, 0 )
forval i = 1/`maxmodels' {
if `numorder`i'' != 0 {
tempname outPlot`i' groupPct`i'
if `numindep`i'' == 4 * `numorder`i'' + 1 + `numdropout`i'' {
matrix outPlot`i' = J( 1 + `numindep`i'', 4 * `numorder`i'' + 1 + `numdropout`i'', . )
}
else {
matrix outPlot`i' = J( `numindep`i'', 4 * `numorder`i'' + 1 + `numdropout`i'', . )
}
matrix groupPct`i' = J( 1, `numorder`i'', . )
}
}
capture program _traj, plugin using("traj.plugin")
/*
`risk3' `expos3' `tcov3'
*/
plugin call _traj `var1' `indep1' `risk1' `weight1' `expos1' ///
`tcov1' `outcome1' `dcov1' `obsmar1' ///
`var2' `indep2' `risk2' `expos2' `tcov2' ///
`var3' `indep3' ///
`outMat' if `touse' `in', ///
`model1' ///
`numindep1' ///
`numrisk1' ///
`numweight1' ///
`numexpos1' ///
`numtcov1' ///
`plottcov1_len' ///
`numorder1' ///
"`min1'" ///
"`max1'" ///
`numiorder1' ///
`numrorder1' ///
`numdropout1' ///
`numdcov1' ///
`numobsmar1' ///
"`omodel1'" ///
`refgroup1' ///
"`model2'" ///
`numindep2' ///
`numrisk2' ///
`numexpos2' ///
`numtcov2' ///
`numorder2' ///
"`min2'" ///
"`max2'" ///
`numiorder2' ///
`refgroup2' ///
"`model3'" ///
`numindep3' ///
`numorder3' ///
"`min3'" ///
"`max3'" ///
`numiorder3' ///
`multgroups' ///
`traj_start_len' ///
`itdetail' ///
`ittrace' ///
`itnovar' ///
`oosprobs'
ereturn clear
local np : word count `outvars'
matrix `b' = parmData[1,1..`np']
matrix `V' = varData[1..`np', 1..`np']
mat colnames `b' = `outvars'
mat rownames `V' = `outvars'
mat colnames `V' = `outvars'
ereturn post `b' `V', e(`touse')
ereturn local cmd "traj"
forval i = 1/`maxmodels' {
if `numorder`i'' != 0 {
mat colnames outPlot`i' = `plotVarLabels`i''
ereturn scalar numGroups`i' = `numorder`i''
ereturn matrix plot`i' = outPlot`i'
ereturn matrix groupSize`i' = groupPct`i'
}
}
ereturn local cmdline `"traj `0'"'
end
/* end of traj.ado */