*** LAB3 2008 ***

*First read in the dataset 

*Calculate each individual's mean meansurement
generate mean_wm = (wm1+wm2)/2
summarize mean_wm

* Save the overall mean and call it wm_mean
local wm_mean = r(mean)

* Plot the data and the individual means/ red line = overall mean
twoway (scatter wm1 id, msymbol(circle)) (scatter wm2 id, symbol(circle_hollow)), xtitle(Subject Id) ytitle(Mini Wright Measurements) legend(order(1 "Occasion 1" 2 "Occasion 2")) yline(`wm_mean')

*Change to long-format for analysis
reshape long wm wp, i(id) j(occasion)

*EDA for normality
sort(id)
by id: egen mean_wm=mean(wm)
hist mean_wm, norm

gen wm_resid = wm-mean_wm
hist wm_resid, norm

*Fit random intercept model ( all three methods are equivalent)
xtreg wm, i(id) mle
xtmixed wm || id:, mle
gllamm wm, i(id) nip(12) adapt

*Predict the random intercept deviations
gllapred eb, u

*Predict the random intercept
gllapred eb, linpred

*Plot the prediction and the data
reshape wide wm wp eb ebm1 ebs1, i(id) j(occasion)
twoway (scatter wm1 id, msymbol(circle)) (scatter wm2 id, msymbol(circle_hollow)) (scatter eb1 id, msymbol(X)), xtitle(Subject Id) ytitle(Mini Wright Measurements) legend(order(1 "Occasion 1" 2 "Occasion 2" 3 "EB Subject-Spec Intercept")) yline(`wm_mean')

*Model checking
hist eb1, norm

reshape long wm wp eb ebm1 ebs1, i(id) j(occasion)
gen eb_resid = wm-eb
hist eb_resid, norm