cd "~/Documents/teaching/c2016/biostat208/labs/lab6" * do "~/Documents/teaching/c2016/biostat208/labs/lab6/lab6.do" capture log close set more off log using "lab6.smcl", replace use lab6, clear * univariate linearity lowess bmd age, name(linearity, replace) graph export linearity.pdf, replace * multivariate linearity reg bmd age weight cprplot weight, lowess name(cprplot1, replace) graph export cprplot1.pdf, replace quietly adjust age, by(weight) gen(fitted1) regress bmd age lweight cprplot lweight, lowess name(cprplot2, replace) nlcom _b[lweight]*log(1.1) graph export cprplot2.pdf, replace quietly adjust age, by(weight) gen(fitted2) regress bmd age weight weight2 quietly adjust age, gen(fitted3) predict residuals2, resid rvpplot weight, yline(0) addplot(lowess residuals2 weight) name(rvpplot3, replace) graph export rvpplot3.pdf, replace mkspline wt = weight, cubic regress bmd age wt1-wt4 quietly adjust age, by(weight) gen(fitted4) twoway (scatter bmd weight, msize(small)) /// (line fitted1 weight, sort lpattern(solid) lwidth(medthick)) /// (line fitted2 weight, sort lpattern(dash_dot) lwidth(thick)) /// (line fitted3 weight, sort lpattern(longdash) lwidth(thick)) /// (line fitted4 weight, sort lpattern(shortdash) lwidth(thick)), /// legend(order(2 "Linear" 3 "Log" 4 "Quadratic" 5 "Cubic spline") rows(2)) /// name(bmdfits, replace) graph export bmdfits.pdf, replace drop fitted* * normality of residuals reg bmd age weight weight2 predict bmdresid, resid qnorm bmdresid, name(qqplot1, replace) graph export qqplot1.pdf, replace kdensity bmdresid, normal graph export kdplot1.pdf, replace reg eeu age poorhlth gaitspd predict eeuresid, res qnorm eeuresid, name(qqplot2, replace) graph export qqplot2.pdf, replace replace eeu = eeu+0.01 qladder eeu, name(qladder, replace) graph export qladder.pdf, replace gen ln_eeu = log(eeu) reg ln_eeu age poorhlth gaitspd predict ln_eeuresid, res qnorm ln_eeuresid, name(qqplot3, replace) graph export qqplot3.pdf, replace * constant variance rvfplot, name(rvfplot1, replace) graph export rvfplot1.pdf, replace rvpplot age, name(rvpplot1, replace) graph export rvpplot1.pdf, replace rvpplot gaitspd, name(rvpplot2, replace) graph export rvpplot2.pdf, replace table poorhlth, c(n ln_eeuresid sd ln_eeuresid) * leverage and influence reg bmd age lweight dfbeta graph box _dfbeta_1 _dfbeta_2, name(boxplot1, replace) graph export boxplot1.pdf, replace list bmd age lweight _dfbeta_1 if abs(_dfbeta_1) > 0.2 & ~missing(_dfbeta_1) list bmd age lweight _dfbeta_2 if abs(_dfbeta_2) > 0.2 & ~missing(_dfbeta_2) reg bmd age lweight if abs(_dfbeta_1) <= .2 & abs(_dfbeta_2) <= 0.2 * checking overlap recode estrogen 1=0 2=1, gen(ht_current) label values ht_current yesno tab estrogen ht_current foreach x in age bmi gaitspd has10 { table ht_current, c(mean `x' p5 `x' min `x' p95 `x' max `x') } foreach x in usearms poorhlth calsupp etid nsfx2 { tab ht_current `x', row } mkspline agesp = age, cubic mkspline bmisp = bmi, cubic mkspline gaitspdsp = gaitspd, cubic mkspline has10sp = has10, cubic sw, pr(.2): logistic ht_current (agesp*) (bmisp*) (has10*) nsfx2 calsupp etid (gaitspdsp*) logistic ht_current agesp* calsupp testparm agesp* predict logit_pscore, xb tab ht_current, sum(logit_pscore) twoway (kdensity logit_pscore if estrogen==1, area(1) lpattern(solid)) /// (kdensity logit_pscore if estrogen==0, area(1) lpattern(longdash)), /// ytitle("Density") xtitle("Logit Propensity Score") /// legend(order(1 "Current estrogen users" 2 "Past and never users")) /// name(pscores, replace) graph export pscores.pdf, replace * optional section on splines capture drop bmicat capture drop fitted* recode bmi min/18.5=1 18.5001/25=2 25.0001/30=3 30.00001/35=4 35.00001/max=5, gen(bmicat) xi: reg bmd i.bmicat age i.estrogen qui adjust age _Iestrogen_*, gen(fitted1) xb * rerun using Version 13 notation to test for trend, departure from linearity reg bmd i.bmicat age i.estrogen * test for linear trend across categories contrast q(1).bmicat * test for departure from linearity contrast q(2/4).bmicat * linear splines capture drop bmi1 bmi2 bmi3 bmi4 bmi5 mkspline bmi1 18.5 bmi2 25 bmi3 30 bmi4 35 bmi5 = bmi xi: reg bmd bmi1-bmi5 age i.estrogen qui adjust age _Iestrogen_*, gen(fitted2) xb * test for departure from linear trend testparm bmi*, equal * cubic splines capture drop bmisp* mkspline bmisp = bmi, cubic xi: reg bmd bmisp* age i.estrogen qui adjust age _Iestrogen_*, gen(fitted3) xb * test for overall BMI effect testparm bmisp* * test for departure from linearity testparm bmisp2 bmisp3 bmisp4 twoway /// (scatter bmd bmi, msize(vtiny)) /// (line fitted1 bmi, sort c(J) lp(longdash) lcol(red) lw(medthick)) /// (line fitted2 bmi, sort lp(shortdash) lcol(green) lw(medthick)) /// (line fitted3 bmi, sort lp(shortdash) lcol(black) lw(medthick)) /// (lowess bmd bmi, lp(solid) lcol(blue) lw(medthick)), /// plotregion(style(none)) scheme(s1color) ytitle("BMD (gm/cm^2)") /// legend(order(2 "Categorical" 3 "Linear Spline" 4 "Cubic Spline" 5 "Lowess") /// rows(2)) caption(Adjusted for age and estrogen use) name(splinefit, replace) graph export splinefit.pdf, replace log close