*! version 2.0.1 JMGarrett 05May98 /* Program to graph predicted values for logistic regression models */ /* 02/20/95 Joanne M. Garrett */ /* Form: logpred y x, from(#) to(#) inc(#) adj(cov_list) options */ /* Options required: from(#), to(#), increment(#) */ /* Options allowed: class(var), poly, adjust, model, graph, nolist */ /* Note: X variable must be continuous (interval or ordinal) */ program define logpred version 4.0 #delimit ; local options "From(real 1) To(real 1) Inc(real 1) Poly(real 0) CLass(string) Model Adjust(string) NOList GRaph T1title(string) T2title(string) L1title(string) L2title(string) Level(real 95) *" ; #delimit cr local varlist "req ex min(2) max(2)" local if "opt" parse "`*'" parse "`varlist'", parse(" ") preserve if ("`class'"~="" & `poly'~=0) { disp in red "Polynomial terms will not work with class() option" exit } capture keep `if' local varlbly : variable label `1' local yvar="`1'" local varlblx : variable label `2' local xvar="`2'" quietly drop if `xvar'==. | `yvar'==. if "`class'"~="" { local clvar="`class'" quietly drop if `clvar'==. local varlblc: variable label `class' } * If there are covariates, drop missing values, calculate means parse "`adjust'", parse(" ") 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'=_result(3) } local i=`i'+1 } keep `yvar' `xvar' `clvar' `covlist' * 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' { quietly gen Xxclss`i'=`xvar' * clss`i' local clist `clist' clss`i' local ilist `ilist' Xxclss`i' local i=`i'+1 } } * 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" } * Run logistic models and test for interaction if class specified if "`model'"~="" { logistic `yvar' `xvar' `polylst' `covlist' `clist' `ilist' more } if "`model'"=="" { quietly logistic `yvar' `xvar' `polylst' `covlist' `clist' `ilist' } local newn=_result(1) * LR test for interaction if present if "`class'"~="" { estimates hold logest local l1=-2*_result(2) quietly logistic `yvar' `xvar' `covlist' `clist' local l2=-2*_result(2) local df=`numcat'-1 local chisq=`l2'-`l1' local probchi=chiprob(`df',`chisq') estimates unhold logest } * 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=_result(2) local i=2 while `i'<=`numcat' { quietly gen Xxclss`i'=`xvar' * clss`i' local clist `clist' clss`i' local ilist `ilist' Xxclss`i' local i=`i'+1 } } * Calculate the predicted values and 95% confidence intervals tempvar linpred predict pred_y predict se, stdp predict `linpred', xb local z=invnorm((1-`level'/100)/2) gen lower=1/(1+exp(-`linpred'-`z'*se)) gen upper=1/(1+exp(-`linpred'+`z'*se)) * List results if "`class'"~="" {sort `clvar' `xvar'} if "`nolist'"~="nolist" { display " " display in green "Probabilities and `level'% Confidence Intervals" display " " display in gr" Outcome:" in yel " `varlbly' -- $S_E_depv" display in gr" X Variable:" in yel" `varlblx' -- `xvar'" if "`class'"~="" { display in gr " Class:" in yel " `varlblc' -- `clvar'" display in gr " Interaction:" in yel " `xvar' by `clvar'" } if `poly'==2 | `poly'==3 { display in gr " Polynomials:" in yel " `polylst'" } if `numcov'>0 {display in gr " Covariates:" in yel " `covdisp'"} if `numcov'==0 {display in gr " Covariates:" in yel " (none)"} display in gr " Observations:" in yel " `newn'" if "`class'"=="" {list `xvar' pred lower upper} if "`class'"~="" { by `clvar': list `xvar' pred lower upper disp " " #delimit ; disp in green " Likelihood ratio test for interaction of" in yellow " `xvar' * `class'" in green ":"; #delimit cr disp " " disp in green " LR Chi2(`df') = " in yellow %6.2f `chisq' if `probchi'>=.0001 { disp in green " Prob > Chi2 = " in yellow %7.4f `probchi' } if `probchi'<.0001 { disp in green " Prob > Chi2 < " in yellow "0.0001" } } if "`graph'"~="" {more} } * Graph results if requested if "`graph'"~="" { if "`l1title'"=="" & "`class'"~="" { if "`varlbly'"=="" { local l1title "Predicted Probabilities for $S_E_depv" } if "`varlbly'"~="" { local l1title "Predicted Probabilities for `varlbly'" } } if "`l1title'"=="" & "`class'"=="" { local l1title "Predicted Probabilities and 95% CI" if "`l2title'"=="" { if "`varlbly'"=="" {local l2title "for $S_E_depv" } if "`varlbly'"~="" {local l2title "for `varlbly'" } } } if "`class'"~="" & "`l2title'"=="" { if "`varlblc'"=="" {local l2title "By Categories of `clvar'"} if "`varlblc'"~="" {local l2title "By Categories of `varlblc'"} } if "`class'"=="" { #delimit ; graph pred_y upper lower `xvar', sort c(sss) s(Oii) l1("`l1title'") l2("`l2title'") `options' ; #delimit cr } if "`class'"~="" { rename pred_y _P local clval : value label `clvar' 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 } reshape groups `clvar' `cval' reshape vars _P reshape cons `xvar' reshape wide 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 } #delimit ; graph _P* `xvar', l1("`l1title'") l2("`l2title'") `options' ; #delimit cr } } end