*** do file for recreating models for behaviour paper
* want to change analysis to include baseline as covariate and no variable selection

set more off 

* data set

cd "G:\My Drive\Research\Active\Grant\Thesis Data sets and Do files\Chapter 2 AFI project\Final Data set"
use "stoddard_cramer_data_stata.dta", clear

* to get baseline as a covariate
keep if time==1
* round baseline values and generated baseline variable
gen time_0_tyield = round(tyield)
gen time_0_tactivity = round(tactivity)
gen time_0_trest = round(trest)
gen time_0_trestbouts = round(trestbouts)

* table 1 data
ci means dim lact tyield tactivity trest trestbouts milk

* of cows with data from each farm and average herdsize
tab farmid, sum(herdsize)


* keep just baseline variables
drop original_cowid time_of_trim lact lactno season dim /// 
tyield trestbouts trimmed herdsize milk time2 dim2 dim3 trestbouts2 trest tactivity

save "JDS_time_0", replace


* merge back into main dataset
use "AFI project final data set 2 7 18.dta", clear
merge m:1 farmid new_cowid using "JDS_time_0.dta"

* drop time 0 from repeated measures and merge variable
drop _merge
*keep if time>1

save "JDS_model", replace

use "JDS_model", clear

* label variables to get better graphs
label variable farmid "Farm"
label variable time2 "Days relative to hoof trimming"
label define Farm 1 "Farm 1" 2 "Farm 2" 3 "Farm 3" 4 "Farm 4"
label values farmid Farm
label define Time 1 "Baseline" 2 "Day of" 3 "+1" 4 "+2" 5 "+3" 6 "+4" 7 "+5" 8 "+6" 9 "+7"
label values time2 Time

* set color scheme
set scheme s2mono


* milk model
mixed tyield i.time2##i.farmid i.lact i.season dim time_0_tyield|| new_cowid:, var noconst residuals(exchangeable,) 
margins time2#farmid
marginsplot, by(farmid) xlabel(, angle(45)) ytitle("Predicted Milk Production") byopts(title("")) plotopts(lcolor("red")) 
estat ic 
contrast i.time2 i.farmid i.lact i.season 

* alternative plotting difference Have to keep baseline in model as it's referent
mixed tyield i.time2##i.farmid i.lact i.season dim time_0_tyield|| ///
new_cowid:, var noconst residuals(exchangeable,) 
margins r.time2@farmid
marginsplot, by(farmid) allsimplelabels ///
 xlabel(2 "0" 3 "+1" 4 "+2" 5 "+3" 6 "+4" 7 "+5" 8 "+6" 9 "+7") ///
 ytitle("Predicted difference from baseline") byopts(title("")) ///
 plotopts(lcolor("red")) yline(0, lpattern(dash))

* resting time
mixed trest i.time2##i.farmid i.lact i.season dim time_0_trest|| ///
new_cowid:, var noconst residuals(exchangeable,) 
margins r.time2@farmid
marginsplot, by(farmid) allsimplelabels ///
 xlabel(2 "0" 3 "+1" 4 "+2" 5 "+3" 6 "+4" 7 "+5" 8 "+6" 9 "+7") ///
 ytitle("Predicted difference from baseline (min)") byopts(title("")) ///
 plotopts(lcolor("red")) yline(0, lpattern(dash))
 
 * activity
mixed tactivity i.time2##i.farmid i.lact i.season dim time_0_tactivity|| ///
new_cowid:, var noconst residuals(exchangeable,) 
margins r.time2@farmid
marginsplot, by(farmid) allsimplelabels ///
 xlabel(2 "0" 3 "+1" 4 "+2" 5 "+3" 6 "+4" 7 "+5" 8 "+6" 9 "+7") ///
 ytitle("Predicted difference from baseline(steps)") byopts(title("")) ///
 plotopts(lcolor("red")) yline(0, lpattern(dash))
 
 * bouts
mixed trestbouts i.time2##i.farmid i.lact i.season dim time_0_trestbouts|| ///
new_cowid:, var noconst residuals(exchangeable,) 
margins r.time2@farmid
marginsplot, by(farmid) allsimplelabels ///
 xlabel(2 "0" 3 "+1" 4 "+2" 5 "+3" 6 "+4" 7 "+5" 8 "+6" 9 "+7") ///
 ytitle("Predicted difference from baseline") byopts(title("")) ///
 plotopts(lcolor("red")) yline(0, lpattern(dash))