************************************************* * Linear and Logistic Regression * Ph. D course , Aarhus University * copyright: Morten Frydenberg * Department of Biostatistics * * Most of the STATA code used at * DAY 6 ************************************************* cd "E:\KURSER\PostReg\Version6\Lectures\Day6" capture log close log using RegDay6.log,replace use obese,clear ******************************************************* * Hosmer-Lemeshow test ******************************************************* generate age45=age-45 xi: logit obese i.sex age45 lfit, group(10) table xi: logit obese i.sex*age45 lfit, group(10) table ******************************************************* * using lincom in with logistic regression ******************************************************* xi: logit obese i.sex age45 lincom _cons+_Isex lincom _cons+_Isex,or lincom _cons+_Isex-age45*3,or lincom age45*4.5,or *********************************************** * The ROC curve and area under the ROC curve *********************************************** generate over45=(age>45) if age!=. diagt obese over45 roctab obese over45,graph tab de aspect(1) xsize(4) graph export RegDay6_01.wmf,replace hist age ,by(obese,col(1))name(p1,replace) nodraw hist age ,by(obese,col(1))name(p2,replace) nodraw freq graph combine p2 p1,col(2) name(p3,replace) graph export RegDay6_02.wmf,replace cumul age,ge(cu) by(obese ) label var cu "1-cumulative dist" sort age cu replace cu=1-cu separate cu ,by(obese ) line cu? pobese age,co(J J J)name(p4,replace)lpa(solid dash) lco(red blue black) /// legend(label(1 "1-Specificty") label(2 "Sensitivity") label(3 "proportion obese") /// ring(0) pos(2)) graph export RegDay6_03.wmf,replace roctab obese age,graph de aspect(1)xsize(4) graph export RegDay6_04.wmf,replace ************************************************** * obese I.sex*age ************************************************** xi:logit obese i.sex*age45 lroc ,aspect(1) xsize(4) graph export RegDay6_05a.wmf,replace lsens graph export RegDay6_05b.wmf,replace predict predprob,p hist predprob ,by(obese,col(1))name(p1,replace) nodraw hist predprob ,by(obese,col(1))name(p2,replace) nodraw freq graph combine p2 p1,col(2) name(p3,replace) graph export RegDay6_05.wmf,replace ************************************************** * euroscore ************************************************** use euroscore,clear hist euroscore,by(dead,col(1))name(p1,replace) nodraw hist euroscore,by(dead,col(1))name(p2,replace) nodraw freq graph combine p2 p1,col(2) name(p3,replace) graph export RegDay6_08.wmf,replace capture drop cu cu? cumul euroscore,ge(cu) by(dead ) label var cu "1-cumulative dist" sort euroscore cu replace cu=1-cu separate cu ,by(dead ) char euroscore [omit] 5 xi:logit dead i.euroscore predict pdead,p line cu? pdead euroscore,co(J J J)name(p4,replace)lpa(solid dash) lco(red blue black ) /// legend(label(1 "1-Specificty") label(2 "Sensitivity") label(3 "proportion dead") ring(0) pos(2)) graph export RegDay6_09.wmf,replace roctab dead euroscore,graph de aspect(1)xsize(4) graph export RegDay6_10.wmf,replace *********************************************** * cancer data *********************************************** use oralcancer,clear xi:logit cancer textile i.alkcon i.genage xi:clogit cancer textile i.alkcon, group(genage) xi:clogit cancer textile i.alkcon, group(genage) or ******************************************************* * Using glm to fit logit, log and identity model ******************************************************* xi: glm obese i.sex age, fam(bin) link(logit) eform * alternative xi:binreg obese i.sex age,or predict plogit if e(sample),mu xi: glm obese i.sex age, fam(bin) link(log) eform * alternative xi:binreg obese i.sex age,rr predict plog if e(sample),mu xi: glm obese i.sex age, fam(bin) link(iden) * alternative xi:binreg obese i.sex age,rd predict pid if e(sample),mu ******************************************************* * A plot comparing the models PLOT01 ******************************************************* twoway /// (line plogit age if sex==1,clco(blue) clpa(1) clw(thick)) /// (line plogit age if sex==2,clco(red) clpa(1) clw(thick)) /// (line plog age if sex==1,clco(blue) clpa(dash) clw(thick)) /// (line plog age if sex==2,clco(red) clpa(dash) clw(thick)) /// (line pid age if sex==1,clco(blue) clpa(dot) clw(thick)) /// (line pid age if sex==2,clco(red) clpa(dot) clw(thick)) /// , ytitle("prevalence") xlabel(30(5)70) name(p2, replace) /// legend(ring(0) pos(11) rows(3) /// label(1 "logit men" ) label(2 "logit women") /// label(3 "log men" ) label(4 "log women") /// label(5 "ident men" ) label(6 "ident women")) graph export RegDay6_01.wmf, replace ******************************************************* log close