set more off capture log close * do "~/Documents/teaching/c2017/biostat208/homework/homework3/hwk3.do" cd "~/Documents/teaching/c2017/biostat208/homework/homework3" log using hwk3, replace set scheme s1mono use hwk3, clear gen lntgl = log(tgl) gen whr10 = whr/0.1 gen lnwhr = log(whr) gen lnbmi = log(bmi) recode bmi min/25 = 1 25.0001/30 = 2 30.0001/35 = 3 35.0001/max = 4, gen(bmicat) recode whr min/.8 = 1 .80001/.85 = 2 .85001/.9 = 3 .90001/max = 4, gen(whrcat) reg tgl age i.raceth exercise drinkany whr10 i.bmicat predict resid, resid reg lntgl age i.raceth exercise drinkany whr10 i.bmicat, eform("exp(beta)") * percent increase in TGL for a 0.1 unit increase in WHR nlcom 100*(exp(_b[whr10])-1) * corresponding absolute increase at mean TGL nlcom 164*(exp(_b[whr10])-1) predict lnresid, resid * normality checks set graphics off qnorm resid, msize(small) name(qn1, replace) ytitle("residuals") qnorm lnresid, msize(small) name(qn2, replace) ytitle("residuals after log transformation") set graphics on graph combine qn1 qn2, name(qqplots, replace) graph export qqplots.pdf, replace * sensitivity analysis dropping smallest TGL value qnorm lnresid if tgl>36, /// msize(small) name(qn2, replace) ytitle("residuals after log transformation") reg lntgl age i.raceth exercise drinkany whr10 i.bmicat if tgl>36, eform("exp(beta)") reg lntgl age i.raceth exercise drinkany whr10 i.bmicat, eform("exp(beta)") * linearity checks reg tgl age i.raceth exercise drinkany whr10 i.bmicat set graphics off cprplot whr10, /// rlopts(clpatt(solid)) lsopts(bw(.5) clpat(longdash) clw(thick)) /// plotregion(style(none)) msize(vtiny) /// ytitle("component plus residual") xtitle("WHR") /// name(cpr1, replace) reg tgl age i.raceth exercise drinkany lnwhr i.bmicat cprplot lnwhr, /// rlopts(clpat(solid)) lsopts(bw(.5) clpat(longdash) clw(thick)) /// plotregion(style(none)) msize(vtiny) /// ytitle("component plus residual") xtitle("ln WHR") /// name(cpr2, replace) set graphics on graph combine cpr1 cpr2, name(cprplots, replace) graph export cprplots.pdf, replace reg tgl age i.raceth exercise drinkany i.whrcat i.bmicat testparm i.whrcat testparm i.whrcat, equal lincom 3.whrcat - 2.whrcat lincom 4.whrcat - 3.whrcat margins whrcat matrix list r(b) gen catfit = . forvalues i = 1/4 { qui replace catfit = el(r(b), 1, `i') if whrcat==`i' } twoway (line catfit whr, sort c(J) lw(thick)), /// ytitle("Mean Triglyceride Level (mg/dl)") /// xtitle("Waist/Hip Ratio") legend(off) /// name(catfit, replace) graph export catfit.pdf, replace qui drop catfit qui xi: reg tgl age i.raceth exercise drinkany i.whrcat i.bmicat adjust age _Irac* exercise drinkany _Ibmi*, gen(catfit) twoway (line catfit whr, sort c(J) lw(thick)), /// ytitle("Mean Triglyceride Level (mg/dl)") /// xtitle("Waist/Hip Ratio") legend(off) /// name(catfit, replace) graph export catfit.pdf, replace * repeat analysis for log TGL reg lntgl age i.raceth exercise drinkany i.whrcat i.bmicat testparm i.whrcat testparm i.whrcat, equal lincom 3.whrcat - 2.whrcat lincom 4.whrcat - 3.whrcat margins whrcat replace catfit = . forvalues i = 1/4 { qui replace catfit = el(r(b), 1,`i') if whrcat==`i' } twoway (line catfit whr, sort c(J) lw(thick)), /// ytitle("Mean Triglyceride Level (mg/dl)") /// xtitle("Waist/Hip Ratio") legend(off) /// name(catfit, replace) graph export catfit.pdf, replace * influence checks capture drop _dfbeta_* reg tgl age i.raceth exercise drinkany i.whrcat i.bmicat dfbeta graph box _dfbeta_*, name(dfbetaboxplot, replace) graph export dfbetaboxplot.pdf, replace graph export dfbetaboxplot.eps, replace capture program drop comparebeta program define comparebeta tempname i nvars nx range omit betaall betaomit diff pctdiff nomit output syntax anything [, cutoff(real 0.2) minx(integer 1) maxx(integer 100)] preserve qui reg `anything' matrix `betaall' = e(b) scalar `nx' = colsof(`betaall') quietly dfbeta quietly dfbeta gen `omit' = 0 forvalues i = 1/100 { if `i' >= `minx' & `i' <= min(`maxx', `nx'-1) { quietly replace `omit' = 1 if abs(_dfbeta_`i') >= `cutoff' & _dfbeta_`i' ~= . } } egen `nomit' = total(`omit') scalar `nomit' = `nomit' qui reg `anything' if `omit' == 0 matrix `betaomit' = e(b) matrix `diff' = `betaomit'[1, 1..`nx'] - `betaall'[1, 1..`nx'] matrix `pctdiff' = `diff' forvalues i = 1/100 { if `i' <= `nx' { matrix `pctdiff'[1, `i'] = `diff'[1, `i'] / `betaall'[1, `i'] * 100 } } matrix `output' = (`betaall'[1, 1..`nx'] \ `betaomit'[1, 1..`nx'] \ `diff' \ `pctdiff')' matrix colnames `output' = "all data" "omitting IPs" "difference" "% difference" display "Differences in beta coefficients after omitting influential points (IPs)" display " with DFbetas > `cutoff' in absolute value" display "Number of omitted IPs:" scalar list `nomit' matrix list `output', format(%6.1f) restore end capture drop maxDF gen maxDF = 0 forvalues i = 1/11 { replace maxDF = abs(_dfbeta_`i') if abs(_dfbeta_`i') > maxDF & _dfbeta_`i' ~= . } list age raceth exercise drinkany whr bmi if maxDF >= 0.35 & maxDF ~= . reg tgl age i.raceth exercise drinkany i.whrcat i.bmicat if maxDF < 0.35 comparebeta tgl age i.raceth exercise drinkany i.whrcat i.bmicat, cutoff(0.35) capture drop maxDF gen maxDF = 0 forvalues i = 6/11 { replace maxDF = abs(_dfbeta_`i') if abs(_dfbeta_`i') > maxDF & _dfbeta_`i' ~= . } list age raceth exercise drinkany whr bmi if maxDF >= 0.25 & maxDF ~= . reg tgl age i.raceth exercise drinkany i.whrcat i.bmicat if maxDF < 0.25 comparebeta tgl age i.raceth exercise drinkany i.whrcat i.bmicat, /// cutoff(0.25) minx(6) maxx(11) * overlap check foreach x in age whr bmi { table exercise, c(mean `x' min `x' p75 `x' p25 `x' max `x') } foreach x in drinkany smoking { tab exercise `x', row nokey } logistic exercise age i.raceth smoking drinkany i.whrcat i.bmicat predict logit_pscore, xb twoway (kdensity logit_pscore if exercise==1, lp(solid) lw(thick)) /// (kdensity logit_pscore if exercise==0, lp(longdash) lw(thick)), /// ytitle("Density") xtitle("Logit Propensity Score") /// legend(order(2 "Exercise: <3 times/week" 1 ">=3 times/week") textfirst) /// name(pscores, replace) graph export pscores.pdf, replace * extra credit mkspline whr1 8 whr2 8.5 whr3 9 whr4 = whr10 mkspline bmi1 20 bmi2 25 bmi3 30 bmi4 35 bmi5 = bmi xi: reg tgl age i.raceth exercise drinkany whr1-whr4 bmi1-bmi5 testparm whr1-whr4, equal testparm whr2-whr4, equal capture drop lsfit quietly adjust age _Iraceth* exercise drinkany bmi1-bmi5, gen(lsfit) twoway (scatter tgl whr, msize(vsmall)) /// (line lsfit whr, sort lpattern(dash_dot) lwidth(thick)) /// (lowess tgl whr, sort lpattern(longdash) lwidth(thick)), /// ytitle("Triglyceride levels (mg/dL)") /// legend(order(2 "Linear spline" 3 "LOWESS") rows(1)) graph export tglfits.pdf, replace log close