*****************************************************************************
* Stata do-file to replicate results given in the following presentation:
*
* TITLE:	Running MIXREGLS from within Stata: The runmixregls command
* AUTHORS:	George Leckie
* LOCATION:	Bristol
* DATE:		27-06-2013
*
* George Leckie
* Centre for Multilevel Modelling, 2013
*****************************************************************************
* This do-file differs slightly from that released with the JSS article.
* We have improved runmixregls over time and that has required are to make 
* a few some edits to this do-file to make sure it still replicates the article
* We document these changes below.
*
*****************************************************************************



********************************************************************************
* Installing runmixregls
********************************************************************************

* Install the runmixregls command from the Statistical Software Components
* (SSC) archive. You only need to run this command once.
//ssc install runmixregls

* If you have already installed runmixregls from the SSC, you can check that
* you are using the latest version by typing the following command. 
* Uncomment this command if you wish to run it (i.e. remove the //).
//. adoupdate runmixregls

* View the runmixregls help file
help runmixregls



****************************************************************************
* Example - Riesby depression data
****************************************************************************

* Load the Riesby data into memory
use "http://www.bristol.ac.uk/cmm/media/runmixregls/reisby", clear

* Describe the data contents
codebook, compact

* Drop cases with missing response
recode hamdep (-9=.)

* Declare the data to be panel data
xtset id

* Plot a Spaghetti plot of the observed data
twoway (line hamdep week, connect(ascending)), xlabel(0(1)5) by(endog)

* Fit the model (cf Hedeker and Nordgren, 2013, page 16)
runmixregls hamdep week endog endweek, between(endog) within(week endog)

* Store the estimation results
estimates store ex1m1

* Refit the model specifying the reffects() and residuals() options to
* retrieve the standardized EB random effects and residuals, respecitvely
* Also specify the iterate() option to specify the maximum number of
* iterations to be 100. Specify the noheader option to suppress the
* table header.
runmixregls hamdep week endog endweek, between(endog) within(week endog) ///
  association(none) reffects(theta) residuals(estd) ///
  iterate(100) noheader

* Store the estimation results
estimates store ex1m2

* Figure 1
scatter theta1 theta0

* Figure 2
histogram estd, width(0.5) start(-3) frequency

* LR test comparing Model 1 and 2
lrtest ex1m1 ex1m2

* Wald test for whether the two patient groups differ significantly or not
test endog endweek

* Calculate the Scale variance
nlcom (sigma2_v: [Scale]sigma^2)



****************************************************************************
exit