*! version 1.1.1 JMGarrett 16May06
/*  Graphs predicted values for linear or logistic models for survey data  */
/*  Form:  svypxcon y, xvar(xvar) f(#) t(#) i(#) options                   */
/*  Options required: xvar(), from(#), to(#), inc(#)                       */
/*  Options allowed: poly, adjust, model, graph, nolist, linear, subpop    */
/*  Note:  X variable should be continuous (interval or ordinal)           */
/*  Updated for version 9.0                                                */

program define svypxcon
  version 9.0
  #delimit ;
    syntax varlist (min=1 max=1) [if] [in], Xvar(varlist)
        [From(real 0) To(real 0) Inc(real 1) Poly(real 0) CLass(string)
         MODel Adjust(string) NOList GRaph LINear Level(real 95) 
         SAVepred(string) SUBpop(string) * ] ;
  #delimit cr
  marksample touse
  markout `touse' `xvar'
  tokenize "`varlist'"
  preserve
  if `from'==0 & `to'==0 {
    disp _n(1) as error "Both " as result "from() " as error "and "  ///
      as result "to() " as error "must be specified"
      exit
      }
  quietly keep if `touse'
  local varlbly : variable label `1'
  local yvar "`1'"
  local varlblx : variable label `xvar'
  if "`poly'"~="0" {
    capture assert `poly'==2 | `poly'==3
    if _rc~=0 {
       disp as error "Error: Poly() option can only contain 2 or 3"
       exit
       }
    }
  capture assert `yvar'==1 | `yvar'==0
  if _rc==0 & "`linear'"=="linear" {
     local reg2cat "01"
     }
  if _rc==0 {
     local regtype="log"
     }
  if _rc~=0 | "`linear'"=="linear" {
     local regtype="lin"
     }
  if "`class'"~="" {
     local clvar="`class'"
     local clval : value label `clvar'
     quietly drop if `clvar'==.
     local varlblc: variable label `class'
     }

* Read in subpop variable if it is an option
  if "`subpop'"~=""  {
    local sub "subpop(`subpop')"
    }

* If there are covariates, drop missing values, calculate means
  tokenize "`adjust'" 
  local numcov 0
  local i 1
  while "`1'"~="" {
    local equal=index("`1'","=")
    if `equal'==0  {
       local cov`i'="`1'"
       local mcov`i'="mean"
       }
    if `equal'~=0  {
       local cov`i'=substr("`1'",1,`equal'-1)
       local mcov`i'=substr("`1'",`equal'+1,length("`1'"))
       }
    quietly drop if `cov`i''==.
    local covlist `covlist' `cov`i''
    local covdisp `covdisp' `1'
    local i=`i'+1
    macro shift
    local numcov=`i'-1
    }
  local i 1
  while `i'<=`numcov' {
    if "`mcov`i''"=="mean" {
      quietly sum `cov`i''
      local mcov`i'=r(mean)
      }
    local i=`i'+1
    }
  * keep `yvar' `xvar' `clvar' `covlist'
  local newn=_N

* If polynomial terms are requested, create them
  if `poly'==2  {
     gen x_sq=`xvar'^2
     local polylst="x_sq"
     }
  if `poly'==3  {
     gen x_sq=`xvar'^2
     gen x_cube=`xvar'^3
     local polylst="x_sq x_cube"
     }
 
* If there is a class variable, set up dummy variables and interactions
  if "`class'"~="" {
     quietly tab `clvar', gen(clss)
     local numcat=_result(2)
     local i 2
     while `i'<=`numcat' {
      *** New section allows interaction terms with polynomials
       quietly gen Xxclss`i'=`xvar' * clss`i'
       if `poly'==2 | `poly'==3 {
          quietly gen Xxclss`i'sq=x_sq * clss`i'
          }
       if `poly'==3 {
          quietly gen Xxclss`i'cb=x_cube * clss`i'
          }
       local clist `clist' clss`i'
       if `poly'==0 {
          local ilist `ilist' Xxclss`i'
          }
       if `poly'==2 {
          local ilist `ilist' Xxclss`i' Xxclss`i'sq
          }
       if `poly'==3 {
          local ilist `ilist' Xxclss`i' Xxclss`i'sq Xxclss`i'cb
          }
       local i=`i'+1
       }
     }

* Run regression models and test for interaction if class specified
  svyset
  if "`model'"~=""  {
    if "`regtype'"=="lin" {
      svy, `sub': regress `yvar' `xvar' `polylst' `covlist' `clist' `ilist' 
      more
      }
    if "`regtype'"=="log" {
      svy, `sub': logistic `yvar' `xvar' `polylst' `covlist' `clist' `ilist'
      more
      }
    }
  if "`model'"==""  {
    if "`regtype'"=="lin" {
      qui svy, `sub': regress `yvar' `xvar' `polylst' `covlist' `clist' `ilist'
      }
    if "`regtype'"=="log" {
      qui svy, `sub': logistic `yvar' `xvar' `polylst' `covlist' `clist' `ilist'
      }
    }
  if "`class'"~="" {
    quietly test `ilist'
    local df1=r(df)
    local df2=r(df_r)
    local f=r(F)
    local probf=r(p)
    }
  else {  
    * tests for x, x_sq, or x_cube
    if `poly'==0 {
       quietly test `xvar'
       }
    if `poly'==2 {
       quietly test x_sq
       }
    if `poly'==3 {
       quietly test x_cube
       }
    local df1=r(df)
    local df2=r(df_r)
    local f=r(F)
    local probf=r(p)
    }

* Save the sample n, pop n, and subpop n's
    if "`sub'"=="" {
       local sampn=e(N)
       local popn=round(e(N_pop),.1)
       }
    if "`sub'"~="" {
       local sampn=e(N_sub)
       local popn=round(e(N_subpop),.1)
       }

* If there is a class variable, retain values for later
  if "`class'"~="" {
     tempvar count
     quietly gen `count'=1
     sort `clvar'
     collapse `count', by(`clvar')
     local class1=`clvar'
     local i 2
       while `i'<=`numcat' {
         local class`i'=`clvar'[_n+`i'-1]
         local i=`i'+1
         }
     }

* Generate the values of x to calculate the predicted values
  drop _all
  local i `from'
  while `i'<`to'  {
    local i=`i'+`inc'
    }
  if `i'>`to'  {
    local to=`i'-`inc'
    }
  local newobs=((`to'-`from')/`inc')+1
  local newobs=round(`newobs',1)
  quietly range `xvar' `from' `to' `newobs'
  label var `xvar' "`varlblx'"
  if `poly'==2  {
     gen x_sq=`xvar'^2
     }
  if `poly'==3  {
     gen x_sq=`xvar'^2
     gen x_cube=`xvar'^3
     }
  local i=1
  while `i'<=`numcov'  {
    quietly gen `cov`i''=`mcov`i''
    local i=`i'+1
    }

* If interaction, expand data, create dummy and interaction variables
  if "`class'"~=""  {
     quietly expand `numcat'
     sort `xvar'
     quietly by `xvar': gen `clvar'=_n
     local i 1
     while `i'<=`numcat'  {
       quietly replace `clvar'=`class`i'' if `clvar'==`i'
       local i=`i'+1
       }
     quietly tab `clvar', gen(clss)
     local numcat=r(r)
     local i=2
     while `i'<=`numcat' {
    *** New section allows interaction terms with polynomials
      quietly gen Xxclss`i'=`xvar' * clss`i'
       if `poly'==2 | `poly'==3 {
          quietly gen Xxclss`i'sq=x_sq * clss`i'
          }
       if `poly'==3 {
          quietly gen Xxclss`i'cb=x_cube * clss`i'
          }
       local clist `clist' clss`i'
       if `poly'==0 {
          local ilist `ilist' Xxclss`i'
          }
       if `poly'==2 {
          local ilist `ilist' Xxclss`i' Xxclss`i'sq
          }
       if `poly'==3 {
          local ilist `ilist' Xxclss`i' Xxclss`i'sq Xxclss`i'cb
          }
       local i=`i'+1
       } 
     }

* Calculate the predicted values and confidence intervals
  if "`regtype'"=="log" {
     predict pred_y, p
     }
     else {
       predict pred_y, xb
       }
  predict se, stdp
  local z=invnorm((1-`level'/100)/2)
  if "`regtype'"=="lin" | "`regtype'"=="med" {
    gen lower=pred_y+`z'*se
    gen upper=pred_y-`z'*se
    }
  if "`regtype'"=="log"  {
    tempvar linpred
    predict `linpred', xb
    gen lower=1/(1+exp(-`linpred'-`z'*se))
    gen upper=1/(1+exp(-`linpred'+`z'*se))
    }

* Print results
  if "`class'"~="" {
     sort `clvar' `xvar'
     }
  display "  "
  if "`regtype'"=="lin" {
     disp as text "Predicted Values and `level'% Confidence Intervals"
     }
  if "`regtype'"=="log" {
     disp as text "Predicted Probabilities and `level'% Confidence Intervals"
     }
  display "  "

  if "`regtype'"=="lin" & "`reg2cat'"=="01" {
      disp "  "
       #delimit  ;
       disp as result "    Warning:" as text "  Y variable is coded as " 
          as result "0" as text " and " as result "1" as text  ". Perhaps" ;
       display as text "               you didn't mean to use the "
          as result "linear" as text " option" ;
       #delimit cr
      disp "  "
      }
     
  if "`regtype'"=="lin" {
    display as text "  Model Type:" as result "    Survey Linear Regression"
    }
  if "`regtype'"=="log" {
    display as text "  Model Type:" as result "    Survey Logistic Regression"
    }
  display as text "  Outcome:" as result "       `varlbly' -- $S_E_depv"
  display as text "  X Variable:" as result"    `varlblx' -- `xvar'"
  if "`class'"~="" {
     display as text "  Class:" as result "         `varlblc' -- `clvar'"
     display as text "  Interaction:" as result "   `xvar' by `clvar'"
     }
  if `poly'==2 | `poly'==3  {
     display as text "  Polynomials:" as result "   `polylst'"
     }
  if `numcov'>0 {
     display as text "  Covariates:" as result "    `covdisp'"
     }
  if `numcov'==0 {
     display as text "  Covariates:" as result "    (none)"
     }
  if "`sub'"=="" {
     display as text "  Sample N:      " as result `sampn'
     display as text "  Population N:  " as result `popn'
     }
  if "`sub'"~="" {
     display as text "  Sub Sample N:  " as result `sampn'
     display as text "  Sub Pop N:     " as result `popn'
     }
  disp "  "
  if "`class'"=="" & "`nolist'"=="" {
     list `xvar' pred lower upper, noob separator(0)
     }
  if "`class'"~="" & "`nolist'"=="" {
    by `clvar': list `xvar' pred lower upper, noob separator(0)
    }

  if "`class'"=="" {
     disp "  "
     if `poly'==0 {
         disp as text "  Wald test for" as result " `xvar'" ///
            as text ":"
         }
     if `poly'==2 {
         disp as text "  Wald test for" as result " x_sq" ///
            as text ":"
         }
     if `poly'==3 {
         disp as text "  Wald test for" as result " x_cube" ///
            as text ":"
         }
     disp "  "
     disp as text "    F(`df1', `df2') =  " as result %6.2f `f'
     if `probf'>=.0001 {
        disp as text "    Prob > F  =  " as result %7.4f `probf'
        }
     if `probf'<.0001 {
        disp as text "    Prob > F   < " as result "0.0001"
        }
    }   

  if "`class'"~="" {
     disp "  "
     #delimit ;
       disp as text "  Wald test for interaction of"
        as result " `xvar' * `class'" as text ":";
     #delimit cr
     disp "  "
     disp as text "    F(`df1', `df2') =  " as result %6.2f `f'
     if `probf'>=.0001 {
       disp as text "    Prob > F   = " as result %7.4f `probf'
       }
     if `probf'<.0001 {
       disp as text "    Prob > F   < " as result "0.0001"
       }
    }

* Save results if requested
  if "`savepred'"~="" {
    keep `clvar' `xvar' pred_y se lower upper
    tempfile tempprd
    quietly save `tempprd'
    }

* Graph results if requested
  
  if "`graph'"~=""  {
    if "`varlbly'"=="" {
       local ytitle "for $S_E_depv"
       }
    if "`varlbly'"~="" {
       local ytitle "for `varlbly'"
       }
    if "`varlblx'"=="" {
       local xtitle "`xvar'"
       }
    if "`varlblx'"~="" {
       local xtitle "`varlblx'"
       }
      if "`poly'"=="2" {
         local xtitle "`xtitle'  (quadratic)"
         }
      if "`poly'"=="3" {
        local xtitle "`xtitle'  (cubic)"
        }
    if "`class'"=="" {
       if "`varlblc'"=="" {
          local leg 'clvar'
          }
       if "`varlblc'"~="" {
          local leg `varlblc'
          }
       if "`regtype'"=="lin" {
          local l2title "Predicted Values and `level'% CI"
          }
       if "`regtype'"=="log" {
          local l2title "Predicted Probabilities and `level'% CI"
          }
  
      twoway (connected pred_y `xvar', sort)                              /// 
          (connected upper `xvar', sort msymbol(none) clcolor(cranberry)  ///
             clpat(dash))                                                 ///                                              
          (line lower `xvar', sort clcolor(cranberry) clpat(dash)),       ///
          ytitle("`ytitle'") l2("`l2title'") xtitle("`xtitle'")           ///
          legend(order(1 "Predicted Value" 2 "`level'% CI"))              ///
          ylabel(, angle(horizontal)) `options'                       
       }

    if "`class'"~="" {
      if "`regtype'"=="lin" {
         local l2title "Predicted Values"
         }
      if "`regtype'"=="log" {
         local l2title "Predicted Probabilities"
         }
      if "`varlblc'"=="" {
         local legtitle="`clvar'"
         }
      if "`varlblc'"~="" {
         local legtitle="`varlblc'"
         }
      if "`clval'"~="" {
        local i 1
        while `i'<=`numcat' {
          local clbl`i' : label `clval' `class`i''
          local i=`i'+1
          }
        }
      local i 2
      local cval `class1'
      while `i'<=`numcat' {
        local cval `cval' `class`i''
        local i=`i'+1
        }
      local i 1
      while `i'<=`numcat' {
        if "`clval'"~=""  {
          local leg`i' `i' "`clbl`i''"
          }
        if "`clval'"==""  {
         local leg`i' `i' "`clvar'=`class`i''"
          }
        local longleg `longleg' `leg`i''
        local i=`i'+1
        }
      rename pred_y _P
      keep _P `xvar' `clvar' upper lower
      quietly reshape wide _P upper lower, i(`xvar') j(`clvar')
      local i 1
      while `i'<=`numcat' {
        if "`clval'"~="" {
           label var _P`class`i'' "`clvar' = `clbl`i''"
           }
           else          {
             label var _P`class`i'' "`clvar' = `class`i''"
             }
        local i=`i'+1
        }
        twoway connect _P* `xvar', ytitle("`ytitle'") xtitle("`xtitle'")  ///
          legend(order("`longleg'") title("`legtitle'", size(default)))   ///
          l2("`l2title'") ylabel(, angle(horizontal)) `options'
      }
    }

if "`savepred'"~="" {
  disp "   "
  use `tempprd'
  save "`savepred'"
  }
end