************************************************* * Linear and Logistic Regression * Ph. D course , Aarhus University * copyright: Morten Frydenberg * Department of Biostatistics * * Most of the STATA code used at * DAY 3 ************************************************* cd "E:\KURSER\PostReg\Lectures\" capture log close log using PostReg3.log,replace ******************************************************* * Plotting the relationship between logit and probability * PLOT01 ******************************************************* twoway /// ( function y=exp(x)/(1+exp(x) ),range(-5 5) clco(black)) /// , xlabel( -5(1)5) xtitle( "logit=ln(odds)") /// ylabel( 0(.1)1) ytitle( "Probability") graph export Reg3_plot01.wmf, replace ******************************************************* * Obsesity and sex ******************************************************* use obese,clear sort age tabu obese sex,col char sex[omit]1 xi: logit obese i.sex xi: logit obese i.sex,or ******************************************************* * Obsesity and age ******************************************************* generate age45=age-45 logit obese age45 logit obese age45,or ******************************************************* * Obsesity and age plot of fitted * PLOT02 and PLOT03 ******************************************************* predict fit1 if e(sample),xb predict p1 if e(sample),p line fit1 age ,clco(black) clw(thick) ytitle("log odds")xlabel(30(5)70) xline(45) graph export Reg3_plot02.wmf, replace line p1 age ,clco(black) clw(thick) ytitle("prevalence")xlabel(30(5)70) xline(45) graph export Reg3_plot03.wmf, replace ******************************************************* * Obsesity and age in groups ******************************************************* egen agegrp7=cut(age), at(0,35,40,45,50,55,60,120) label table agegrp7 ,c(min age max age count obese sum obese)row char agegrp7[omit]0 xi: logit obese i.agegrp7 xi: logit obese i.agegrp7,or testparm _Iage* char agegrp7[omit]3 xi: logit obese i.agegrp7 xi: logit obese i.agegrp7,or ******************************************************* * Obsesity and age in groups plot of fitted * PLOT04 ******************************************************* predict fit2 if e(sample),xb predict p2 if e(sample),p line fit2 age ,clco(black) connect(J) clpa(1) clw(thick) /// , ytitle("log odds")xlabel(30(5)70) saving(dum1, replace) line p2 age ,clco(black) connect(J) clpa(1) clw(thick) /// , ytitle("prevalence") xlabel(30(5)70) saving(dum2, replace) graph combine dum1.gph dum2.gph graph export Reg3_plot04.wmf, replace ******************************************************* * Obsesity and age comparing the two modeæs * PLOT05 ******************************************************* twoway /// (line fit1 age ,clco(black) clw(thick) clpa(1) ) /// (line fit2 age ,clco(blue) connect(J) clw(thick) clpa(dash)) /// , ytitle("log odds")xlabel(30(5)70) saving(dum1, replace) /// legend(ring(0) pos(11) label(1 "model1" )label(2 "model2") ) twoway /// (line p1 age ,clco(black) clpa(1) clw(thick)) /// (line p2 age ,clco(blue) connect(J) clw(thick) clpa(dash)) /// , ytitle("prevalence") xlabel(30(5)70) saving(dum2, replace) /// legend(ring(0) pos(11) label(1 "model1" )label(2 "model2") ) graph combine dum1.gph dum2.gph graph export Reg3_plot05.wmf, replace ******************************************************* * Obsesity and sex and age no interaction ******************************************************* xi: logit obese i.sex age45 xi: logit obese i.sex age45,or ******************************************************* * PLOT06 ******************************************************* predict fit3 if e(sample),xb predict p3 if e(sample),p twoway /// (line fit3 age if sex==1 & e(sample) ,clco(blue) clpa(1) clw(thick) ) /// (line fit3 age if sex==2 & e(sample) ,clco(red) clpa(dash) clw(thick) ) /// , ytitle("log odds")xlabel(30(5)70) saving(dum1, replace)xline(45) /// legend(ring(0) pos(11) label(1 "men" )label(2 "women") ) twoway /// (line p3 age if sex==1 & e(sample) ,clco(blue) clpa(1) clw(thick) ) /// (line p3 age if sex==2 & e(sample) ,clco(red) clpa(dash) clw(thick) ) /// , ytitle("prevalence")xlabel(30(5)70) saving(dum2, replace) xline(45) /// legend(ring(0) pos(11) label(1 "men" )label(2 "women") ) graph combine dum1.gph dum2.gph graph export Reg3_plot06.wmf, replace ******************************************************* * Obsesity and sex and age with interaction ******************************************************* xi: logit obese i.sex*age45 xi: logit obese i.sex*age45,or ******************************************************* * PLOT07 ******************************************************* predict fit4 if e(sample),xb predict p4 if e(sample) twoway /// (line fit4 age if sex==1 & e(sample) ,clco(blue) clpa(1) clw(thick) ) /// (line fit4 age if sex==2 & e(sample) ,clco(red) clpa(dash) clw(thick) ) /// , ytitle("log odds")xlabel(30(5)70) xline(45) saving(dum1, replace) /// legend(ring(0) pos(11) label(1 "men" )label(2 "women")) twoway /// (line p4 age if sex==1 & e(sample) ,clco(blue) clpa(1) clw(thick) ) /// (line p4 age if sex==2 & e(sample) ,clco(red) clpa(dash) clw(thick) ) /// , ytitle("prevalence")xlabel(30(5)70) xline(45)saving(dum2, replace) /// legend(ring(0) pos(11) label(1 "men" )label(2 "women")) graph combine dum1.gph dum2.gph graph export Reg3_plot07.wmf, replace *********************************************** * cancer data *********************************************** use case_control,clear tabodds cancer age tabodds cancer age, or tabodds cancer age, or base(3) char age [omit]1 xi:logit cancer i.smoker i.age,or char age [omit]3 xi:logit cancer i.smoker i.age,or estimates store model1 xi:logit cancer i.smoker estimates store model2 lrtest model1 model2 *********************************************** log close