* do-file for lecture 10 of VHM 802, Winter 2014 (somewhat expanded for wetland ex)
version 13
set more off
cd "h:\vhm\vhm802\data_csv"

* Example 16.5
insheet using ch16ta1.csv, clear
anova errors anxiety##tension / subject|anxiety#tension anxiety##tension##type 
regress
* residual analysis
predict pred, xb
predict stdres, rstandard
scatter stdres pred
qnorm stdres
swilk stdres

* random effect checks
preserve /* run next 7 lines together with this one */
collapse (mean) errors, by(anxiety tension subject)
anova errors anxiety##tension
predict pred2, xb
predict stdres2, rstandard
scatter stdres2 pred2
qnorm stdres2
swilk stdres2

* likelihood-based analysis
egen id=group(anxiety tension subject)
mixed errors anxiety##tension##type || id:, reml
* residual analysis
predict pred2, fitted
predict stdres2, rstandard
scatter stdres2 pred2
qnorm stdres2
predict pred_ref, reffects
bysort id: gen within_n=_n
qnorm pred_ref if within_n==1
swilk pred_ref if within_n==1

* Example 16.7
insheet using ch16ta2.csv, clear
anova biomass n table / n#table n##weed / n#table#weed n##weed##clip
predict pred, xb
predict stdres, rstandard
scatter stdres pred
qnorm stdres
* computation of means for model checking:
* collapse (mean) biomass, by(table n weed)
* collapse (mean) biomass, by(table n)

* likelihood-based analysis
mixed biomass table n##weed##clip || tray: || wetland:, reml
testparm n#weed#clip
margins n#weed clip
* residuals analysis
predict pred, xb
predict stdres, rstandard
scatter stdres pred
qnorm stdres
swilk stdres
predict pred_ref*, reffects
bysort tray: gen within1_n=_n
bysort wetland: gen within2_n=_n
qnorm pred_ref1 if within1_n==1
swilk pred_ref1 if within1_n==1
qnorm pred_ref2 if within2_n==1
swilk pred_ref2 if within2_n==1