* do-file for review/tutorial 3b of VHM 812/802, Winter 2022
version 17 /* mostly works also with versions 14-16 */
set more off
cd "r:\"
use daisy2red, clear

* 1a)
regress wpc i.herd
margins herd
lincom _cons+2.herd
tabstat wpc, statistics(mean sd semean count) by(herd)

* 1b)
regress wpc milk120
margins , at( milk120=(1200(1000)5600) )
marginsplot
lincom _cons+milk120*3200
sum milk120, d

* 1c)
regress wpc c.milk120##c.milk120
margins , at( milk120=(1200(1000)5600) )
marginsplot
lincom _cons+milk120*3200+c.milk120#c.milk120*10240000

* 1d)
generate lnwpc=ln(wpc)
regress lnwpc milk120
margins , at( milk120=(1200(1000)5600))
margins , at( milk120=(1200(1000)5600)) expression(exp(predict(xb)))
marginsplot
di exp(3.934371) "  ,  " exp(4.167394) /* correct CI for 1200 */

* 2a)
regress wpc i.rp i.vag_disch
margins rp vag_disch
table rp vag_disch
* Stata 16: table rp vag_disch, row col
lincom _cons+1.rp+1.vag_disch*82/1574 /* rp=1 */
margins rp vag_disch, asbalanced
lincom _cons+1.rp+1.vag_disch*0.5 /* rp=1 */
margins rp, over(vag_disch) /* same as: at(vag_disch=(0 1)) */
marginsplot, noci
lincom _cons+1.rp+1.vag_disch /* rp=1, vag_disch=1 */

* 2b)
regress wpc rp##vag_disch
margins rp#vag_disch
marginsplot, noci /* this is the interaction plot! */
marginsplot, noci x(vag_disch) /* controlling variable on x */
marginsplot, noci x(rp) /* same as default */
marginsplot, noci nosimplelabels /* more informative labels */
table rp vag_disch, statistic(mean wpc)
* Stata 16: table rp vag_disch, contents(mean wpc) row col
margins rp
lincom _cons+1.rp+(1.vag_disch+1.rp#1.vag_disch)*82/1574  /* rp=1 */
margins rp, asbalanced
lincom _cons+1.rp+(1.vag_disch+1.rp#1.vag_disch)*0.5  /* rp=1 */

* 2c)
regress wpc i.dyst milk120
margins dyst, at( milk120=(1200 2200 3200 4300 5500))
marginsplot, noci
lincom _cons+1.dyst+milk120*3200 /* dyst=1, milk120=3200 */
margin dyst, atmeans
lincom _cons+1.dyst+milk120*3215.096
* interaction model
regress wpc dyst##c.milk120
margins dyst, at( milk120=(1200 2200 3200 4300 5500))
marginsplot, noci /* this is the interaction plot! */
lincom _cons+1.dyst+(c.milk120+1.dyst#c.milk120)*3200 /* dyst=1, milk120=3200 */

* 2d)
regress wpc parity milk120
tabulate parity
margins , at( parity=(1(1)6) milk120=(1200 2200 3200 4300 5500))
marginsplot, noci
margins , at( milk120=(1200 2200 3200 4300 5500) parity=(1(1)6) )
marginsplot, noci /* changing roles in plot */
margins , at( milk120=(1200 2200 3200 4300 5500) (median)parity)
marginsplot
lincom _cons+milk120*1200+parity*2 /* milk120=1200, parity=2 */
margins, atmeans
lincom _cons+milk120*3215.096+parity*2.736328 /* both at means */
* interaction model
regress wpc c.parity##c.milk120
margins , at( parity=(1(1)6) milk120=(1200 2200 3200 4300 5500))
marginsplot, noci
lincom _cons+milk120*1200+parity*2+milk120#c.parity*2400 /* milk120=1200, parity=2 */

* VER Ex. 14.16: estimates on sqrt-scale and backtransformed
gen calv_mth=month(calv_dt)
gen aut_calv=(calv_mth>=2 & calv_mth<=7) if !missing(calv_mth)
gen rootwpc=sqrt(wpc)
regress rootwpc aut_calv c.herd_size##c.herd_size parity twin dyst rp##vag_disch
*1) correct estimate and CI for rp=1, vag_disch=0, dyst=1
lincom _cons+aut_calv*0+herd_size*251+c.herd_size#c.herd_size*63001+parity*1+twin*0+dyst*1+1.rp+0.vag_disch+1.rp#0.vag_disch
*1) correct estimate and SE using margins command
margins , over(rp vag_disch dyst) at(parity=1 aut_calv=0 herd_size=251 twin=0)
*1) correct backtransformed estimate, but inappropriate CI
margins , over(rp vag_disch dyst) at(parity=1 aut_calv=0 herd_size=251 twin=0) expr(predict(xb)^2)
table (rp) (vag_disch dyst), statistic(count wpc) nototals
* Stata 16: table rp dyst vag_disch, contents(count wpc)
table parity
table twin
sum herd_size
table aut_calv
*3) correct estimate and SE using margins command
margins , at(parity=1 aut_calv=0 herd_size=(125 185 201 235 263 294 333) twin=0 rp=0 dyst=0 vag_disch=0)
table herd_size
*3) correct backtransformed estimate, but inappropriate CI
margins , at(parity=1 aut_calv=0 herd_size=(125 185 201 235 263 294 333) twin=0 rp=0 dyst=0 vag_disch=0) expr(predict(xb)^2)
marginsplot, noci