***Exercise 6.2  Postnatal-depression data

*These data are from a 1996 study (Gregoire, Kumar Everitt, Henderson & Studd) on the efficacy of estrogen patches in treating postnatal depression. 
*Women were randomly assigned to either a placebo control group (group=0, n=27) or estrogen patch group (group=1, n=34). 
*Prior to the first treatment all patients took the Edinburgh Postnatal Depression Scale (EPDS).
*EPDS data was collected monthly for six months once the treatment began. Higher scores on the EDPS are indicative of higher levels of depression.

* Read data set
use http://www.stata-press.com/data/mlmus3/postnatal, clear

* 1.
reshape long dep, i(subj) j(month)

* Recode missing values
mvdecode dep pre, mv(-9)

* Describe structure of missing data
xtset subj month
xtdescribe if dep<.

generate time = month - 1

* Extra: 
* Simple linear regression ignoring subject effects and correleted within-subject error terms.
reg dep pre group time
reg dep pre group time,vce(robust)
reg dep pre group time,vce(cluster subj)


* 2.
* Unstructured
xtmixed dep pre group time || subj:, noconstant residuals(unstructured, t(month)) mle
estimates store un

* 3.
* Exchangeable
xtmixed dep pre group time || subj:, noconstant residuals(exchangeable) mle
estimates store exch
lrtest exch un

* 4.
* RI
xtmixed dep pre group time || subj:, mle variance
estimates store ri
display sqrt(14.48409 +11.20199)
display 14.48409/(14.48409 +11.20199)
 
* 5.
* RI & AR(1)
xtmixed dep pre group time || subj:, residuals(ar 1, t(month)) mle
estimates store ri_ar1
lrtest ri_ar1 ri

* 6.
* Toeplitz
xtmixed dep pre group time || subj:, noconstant residuals(toeplitz 5, t(month)) mle
estimates store toep
lrtest toep ri_ar1, force
lrtest toep ri, force  /* or exchangeable */
lrtest toep un

* 7.
* RC
xtmixed dep pre group time || subj: time, covariance(unstructured) mle
estimates store rc
lrtest rc ri

* 8.
* RC & AR(1)
xtmixed dep pre group time || subj: time,  covariance(unstructured) ///
   residuals(ar 1, t(time)) mle
estimates store rc_ar1
lrtest rc_ar1 rc 
lrtest rc_ar1 ri_ar1
lrtest rc_ar1 ri
 
 * 9.
* Comparison
estimates stats un exch ri ri_ar1 toep rc rc_ar1 

* Extra: 

* RI with time as factor: 
xtmixed dep pre group i.time || subj:, mle variance
estimates store rifac
display sqrt(14.48409 +11.20199)
display 14.48409/(14.48409 +11.20199)

* Display correlation matrix: 
ssc install xtmixed_corr, replace /* installs program */

* corr.matrix for first subject:
xtmixed_corr


* Residual diagnostics: 

predict res, residuals
histogram res, normal xtitle("Predicted level-1 residuals")
swilk res

* qq-plott: 
qnorm res

* Fitted values: 
predict pred, fitted


* Contrasts:
lincom _cons - 1.time
lincom 1.time - 2.time
lincom 2.time - 3.time
lincom 3.time - 4.time
lincom 4.time - 5.time