* Prevent "more" messages from appearing
set more off 
* Control line length
set linesize 150 
 
********************************************************************************
*******            BEGIN DATA MANIPULATION OF CHAPTER 6 EXAMPLE          *******
*******               CHANGE "filesave" to your directory                *******
********************************************************************************
 
* Defining global variable for file location to be replaced in code below
global filesave "C:\Dropbox\PilesOfVariance\Chapter6\STATA" 
 
* Import chapter 6 stacked data and create time predictor variables
use "$filesave\STATA_Chapter6.dta", clear 
 
* Predictors for polynomial time models (need to make quadratic time)
gen time1 = session - 1 
gen time6 = session - 6 
label variable time1  "time1: Session (0=1)" 
label variable time6  "time6: Session (0=6)" 
gen time1sq = time1 * time1 
gen time6sq = time6 * time6 
label variable time1sq "time1sq: Quadratic Session (0=1)" 
label variable time6sq "time6sq: Quadratic Session (0=6)" 
* Predictors for piecewise time models
gen slope12 = session 
gen slope26 = session 
recode slope12 (1=0) if session==1 
recode slope12 (2=1) if session==2 
recode slope12 (3=1) if session==3 
recode slope12 (4=1) if session==4 
recode slope12 (5=1) if session==5 
recode slope12 (6=1) if session==6 
recode slope26 (1=0) if session==1 
recode slope26 (2=0) if session==2 
recode slope26 (3=1) if session==3 
recode slope26 (4=2) if session==4 
recode slope26 (5=3) if session==5 
recode slope26 (6=4) if session==6 
label variable slope12 "slope12: Early Practice Slope (Session 1-2)" 
label variable slope26 "slope26: Later Practice Slope (Session 2-6)" 
 
********************************************************************************
*******                      BEGIN CHAPTER 6 MODELS                      *******
*******    NOTE: NEGATIVE EXPONENTIAL MODELS ARE NOT POSSIBLE IN STATA   *******
********************************************************************************
 
* Save results to separate file
log using $filesave\STATA_Chapter6_Output, replace name(STATA_Chapter6) 
 
display as result "Chapter 6 Example: Means by session for RT outcome" 
tabulate session, summarize(rt) 
 
display as result "Ch 6: 0: Saturated Means, Unstructured Variance Model" 
display as result "TOTAL ANSWER KEY" 
mixed rt i.session, ///
           || personid: , noconstant variance reml covariance(unstructured) ///
           residuals(unstructured,t(session)), 
      estat ic, n(101), 
      estat wcorrelation, covariance, 
      estat wcorrelation, 
      contrast i.session,  
      margins  i.session,  
      margins  i.session, pwcompare(pveffects) 
 
display as result "Ch 6: 1a: Empty Means, E-Only Variance Model" 
mixed rt , ///
           || personid: , noconstant variance reml covariance(unstructured) ///
           residuals(independent,t(session)), 
      estat ic, n(101), 
      estat wcorrelation, covariance, 
      estat wcorrelation, 
 
display as result "Ch 6: 1b: Empty Means, Random Intercept Model" 
mixed rt , ///
           || personid: , variance reml covariance(unstructured) ///
           residuals(independent,t(session)), 
      estat ic, n(101), 
      estat icc, 
      estat recovariance, relevel(personid), 
      estat recovariance, relevel(personid) correlation, 
      estat wcorrelation, covariance, 
      estat wcorrelation, 
 
display as result "Ch 6: 2a: Fixed Linear Time, Random Intercept Model" 
mixed rt c.time1, ///
           || personid: , variance reml covariance(unstructured) ///
           residuals(independent,t(session)), 
      estat ic, n(101), 
      estat recovariance, relevel(personid), 
      estat recovariance, relevel(personid) correlation, 
      estat wcorrelation, covariance, 
      estat wcorrelation, 
      * Intercept at Session=1 Time=0
      lincom _cons*1 + c.time1*0 
      * Intercept at Session=2 Time=1
      lincom _cons*1 + c.time1*1 
      * Intercept at Session=3 Time=2
      lincom _cons*1 + c.time1*2 
      * Intercept at Session=4 Time=3
      lincom _cons*1 + c.time1*3 
      * Intercept at Session=5 Time=4
      lincom _cons*1 + c.time1*4 
      * Intercept at Session=6 Time=5
      lincom _cons*1 + c.time1*5 
      estimates store FitFixLin, 
 
display as result "Ch 6: 2b: Random Linear Time Model" 
mixed rt c.time1, ///
           || personid: time1, variance reml covariance(unstructured) ///
           residuals(independent,t(session)), 
      estat ic, n(101), 
      estat recovariance, relevel(personid), 
      estat recovariance, relevel(personid) correlation, 
      estat wcorrelation, covariance, 
      estat wcorrelation, 
      * Intercept at Session=1 Time=0
      lincom _cons*1 + c.time1*0 
      * Intercept at Session=2 Time=1
      lincom _cons*1 + c.time1*1 
      * Intercept at Session=3 Time=2
      lincom _cons*1 + c.time1*2 
      * Intercept at Session=4 Time=3
      lincom _cons*1 + c.time1*3 
      * Intercept at Session=5 Time=4
      lincom _cons*1 + c.time1*4 
      * Intercept at Session=6 Time=5
      lincom _cons*1 + c.time1*5 
      estimates store FitRandLin, 
      lrtest FitRandLin FitFixLin,  
 
display as result "Ch 6: 3a: Fixed Quadratic, Random Linear Time Model" 
mixed rt c.time1 c.time1#c.time1, ///
           || personid: time1, variance reml covariance(unstructured) ///
           residuals(independent,t(session)), 
      estat ic, n(101), 
      estat recovariance, relevel(personid), 
      estat recovariance, relevel(personid) correlation, 
      estat wcorrelation, covariance, 
      estat wcorrelation, 
      * Intercept at Session=1 Time=0
      lincom _cons*1 + c.time1*0 + c.time1#c.time1*0 
      * Intercept at Session=2 Time=1
      lincom _cons*1 + c.time1*1 + c.time1#c.time1*1 
      * Intercept at Session=3 Time=2
      lincom _cons*1 + c.time1*2 + c.time1#c.time1*4 
      * Intercept at Session=4 Time=3
      lincom _cons*1 + c.time1*3 + c.time1#c.time1*9 
      * Intercept at Session=5 Time=4
      lincom _cons*1 + c.time1*4 + c.time1#c.time1*16 
      * Intercept at Session=6 Time=5
      lincom _cons*1 + c.time1*5 + c.time1#c.time1*25 
      * Linear Slope at Session=1 Time=0
      lincom c.time1*1 + c.time1#c.time1*0 
      * Linear Slope at Session=2 Time=1
      lincom c.time1*1 + c.time1#c.time1*2 
      * Linear Slope at Session=3 Time=2
      lincom c.time1*1 + c.time1#c.time1*4 
      * Linear Slope at Session=4 Time=3
      lincom c.time1*1 + c.time1#c.time1*6 
      * Linear Slope at Session=5 Time=4
      lincom c.time1*1 + c.time1#c.time1*8 
      * Linear Slope at Session=6 Time=5
      lincom c.time1*1 + c.time1#c.time1*10 
      estimates store FitFixQuad, 
 
display as result "Ch 6: 3b: Random Quadratic Time Model" 
mixed rt c.time1 c.time1#c.time1, ///
           || personid: time1 time1sq, variance reml covariance(unstructured) ///
           residuals(independent,t(session)), 
      estat ic, n(101), 
      estat recovariance, relevel(personid), 
      estat recovariance, relevel(personid) correlation, 
      estat wcorrelation, covariance, 
      estat wcorrelation, 
      * Intercept at Session=1 Time=0
      lincom _cons*1 + c.time1*0 + c.time1#c.time1*0 
      * Intercept at Session=2 Time=1
      lincom _cons*1 + c.time1*1 + c.time1#c.time1*1 
      * Intercept at Session=3 Time=2
      lincom _cons*1 + c.time1*2 + c.time1#c.time1*4 
      * Intercept at Session=4 Time=3
      lincom _cons*1 + c.time1*3 + c.time1#c.time1*9 
      * Intercept at Session=5 Time=4
      lincom _cons*1 + c.time1*4 + c.time1#c.time1*16 
      * Intercept at Session=6 Time=5
      lincom _cons*1 + c.time1*5 + c.time1#c.time1*25 
      * Linear Slope at Session=1 Time=0
      lincom c.time1*1 + c.time1#c.time1*0 
      * Linear Slope at Session=2 Time=1
      lincom c.time1*1 + c.time1#c.time1*2 
      * Linear Slope at Session=3 Time=2
      lincom c.time1*1 + c.time1#c.time1*4 
      * Linear Slope at Session=4 Time=3
      lincom c.time1*1 + c.time1#c.time1*6 
      * Linear Slope at Session=5 Time=4
      lincom c.time1*1 + c.time1#c.time1*8 
      * Linear Slope at Session=6 Time=5
      lincom c.time1*1 + c.time1#c.time1*10 
      estimates store FitRandQuad, 
      lrtest FitRandQuad FitFixQuad,  
 
display as result "Ch 6: 3b: Random Quadratic Time Model (0=Session 6)" 
mixed rt c.time6 c.time6#c.time6, ///
           || personid: time6 time6sq, variance reml covariance(unstructured) ///
           residuals(independent,t(session)), 
      estat ic, n(101), 
      estat recovariance, relevel(personid), 
      estat recovariance, relevel(personid) correlation, 
      estat wcorrelation, covariance, 
      estat wcorrelation, 
      * Intercept at Session=1 Time=-5
      lincom _cons*1 + c.time6*-5 + c.time6#c.time6*25 
      * Intercept at Session=2 Time=-4
      lincom _cons*1 + c.time6*-4 + c.time6#c.time6*16 
      * Intercept at Session=3 Time=-3
      lincom _cons*1 + c.time6*-3 + c.time6#c.time6*9 
      * Intercept at Session=4 Time=-2
      lincom _cons*1 + c.time6*-2 + c.time6#c.time6*4 
      * Intercept at Session=5 Time=-1
      lincom _cons*1 + c.time6*-1 + c.time6#c.time6*1 
      * Intercept at Session=6 Time= 0
      lincom _cons*1 + c.time6*0 + c.time6#c.time6*0 
      * Linear Slope at Session=1 Time=-5
      lincom c.time6*1 + c.time6#c.time6*-10 
      * Linear Slope at Session=2 Time=-4
      lincom c.time6*1 + c.time6#c.time6*-8 
      * Linear Slope at Session=3 Time=-3
      lincom c.time6*1 + c.time6#c.time6*-6 
      * Linear Slope at Session=4 Time=-2
      lincom c.time6*1 + c.time6#c.time6*-4 
      * Linear Slope at Session=5 Time=-1
      lincom c.time6*1 + c.time6#c.time6*-2 
      * Linear Slope at Session=6 Time= 0
      lincom c.time6*1 + c.time6#c.time6*0 
 
display as result "Ch 6: 4a: Fixed Slope12, Fixed Slope26, Random Intercept Model" 
mixed rt c.slope12 c.slope26, ///
           || personid: , variance reml covariance(unstructured) ///
           residuals(independent,t(session)), 
      estat ic, n(101), 
      estat recovariance, relevel(personid), 
      estat recovariance, relevel(personid) correlation, 
      estat wcorrelation, covariance, 
      estat wcorrelation, 
      * Intercept at Session=1 Time=0
      lincom _cons*1 + c.slope12*0 + c.slope26*0 
      * Intercept at Session=2 Time=1
      lincom _cons*1 + c.slope12*1 + c.slope26*0 
      * Intercept at Session=3 Time=2
      lincom _cons*1 + c.slope12*1 + c.slope26*1 
      * Intercept at Session=4 Time=3
      lincom _cons*1 + c.slope12*1 + c.slope26*2 
      * Intercept at Session=5 Time=4
      lincom _cons*1 + c.slope12*1 + c.slope26*3 
      * Intercept at Session=6 Time=5
      lincom _cons*1 + c.slope12*1 + c.slope26*4 
      * Difference between slope12 and slope26
      lincom c.slope12*-1 + c.slope26*1 
      estimates store FitFix12Fix26, 
 
display as result "Ch 6: 4b: Random Slope12, Fixed Slope26 Model" 
mixed rt c.slope12 c.slope26, ///
           || personid: slope12, variance reml covariance(unstructured) ///
           residuals(independent,t(session)), 
      estat ic, n(101), 
      estat recovariance, relevel(personid), 
      estat recovariance, relevel(personid) correlation, 
      estat wcorrelation, covariance, 
      estat wcorrelation, 
      * Intercept at Session=1 Time=0
      lincom _cons*1 + c.slope12*0 + c.slope26*0 
      * Intercept at Session=2 Time=1
      lincom _cons*1 + c.slope12*1 + c.slope26*0 
      * Intercept at Session=3 Time=2
      lincom _cons*1 + c.slope12*1 + c.slope26*1 
      * Intercept at Session=4 Time=3
      lincom _cons*1 + c.slope12*1 + c.slope26*2 
      * Intercept at Session=5 Time=4
      lincom _cons*1 + c.slope12*1 + c.slope26*3 
      * Intercept at Session=6 Time=5
      lincom _cons*1 + c.slope12*1 + c.slope26*4 
      * Difference between slope12 and slope26
      lincom c.slope12*-1 + c.slope26*1 
      estimates store FitRand12Fix26, 
      lrtest FitRand12Fix26 FitFix12Fix26,  
 
display as result "Ch 6: 4c: Random Slope12, Random Slope26 Model" 
mixed rt c.slope12 c.slope26, ///
           || personid: slope12 slope26, variance reml covariance(unstructured) ///
           residuals(independent,t(session)), 
      estat ic, n(101), 
      estat recovariance, relevel(personid), 
      estat recovariance, relevel(personid) correlation, 
      estat wcorrelation, covariance, 
      estat wcorrelation, 
      * Intercept at Session=1 Time=0
      lincom _cons*1 + c.slope12*0 + c.slope26*0 
      * Intercept at Session=2 Time=1
      lincom _cons*1 + c.slope12*1 + c.slope26*0 
      * Intercept at Session=3 Time=2
      lincom _cons*1 + c.slope12*1 + c.slope26*1 
      * Intercept at Session=4 Time=3
      lincom _cons*1 + c.slope12*1 + c.slope26*2 
      * Intercept at Session=5 Time=4
      lincom _cons*1 + c.slope12*1 + c.slope26*3 
      * Intercept at Session=6 Time=5
      lincom _cons*1 + c.slope12*1 + c.slope26*4 
      * Difference between slope12 and slope26
      lincom c.slope12*-1 + c.slope26*1 
      estimates store FitRand12Rand26, 
      lrtest FitRand12Rand26 FitRand12Fix26,  
 
display as result "Ch 6: Random Slope12, Random Slope26 Model + Fixed Quadratic Slope26" 
mixed rt c.slope12 c.slope26 c.slope26#c.slope26, ///
           || personid: slope12 slope26, variance reml covariance(unstructured) ///
           residuals(independent,t(session)), 
      estat ic, n(101), 
      estat recovariance, relevel(personid), 
      estat recovariance, relevel(personid) correlation, 
      estat wcorrelation, covariance, 
      estat wcorrelation, 
 
display as result "Ch 6: 5a: Fixed Time, Fixed Slope26, Random Intercept Model" 
mixed rt c.time1 c.slope26, ///
           || personid: , variance reml covariance(unstructured) ///
           residuals(independent,t(session)), 
      estat ic, n(101), 
      estat recovariance, relevel(personid), 
      estat recovariance, relevel(personid) correlation, 
      estat wcorrelation, covariance, 
      estat wcorrelation, 
      * Intercept at Session=1 Time=0
      lincom _cons*1 + c.time1*0 + c.slope26*0 
      * Intercept at Session=2 Time=1
      lincom _cons*1 + c.time1*1 + c.slope26*0 
      * Intercept at Session=3 Time=2
      lincom _cons*1 + c.time1*2 + c.slope26*1 
      * Intercept at Session=4 Time=3
      lincom _cons*1 + c.time1*3 + c.slope26*2 
      * Intercept at Session=5 Time=4
      lincom _cons*1 + c.time1*4 + c.slope26*3 
      * Intercept at Session=6 Time=5
      lincom _cons*1 + c.time1*5 + c.slope26*4 
      * Rate of change from session 2 to 6
      lincom c.time1*1 + c.slope26*1 
      estimates store FitFix16Fix26, 
 
display as result "Ch 6: 5b: Random Time, Fixed Slope26 Model" 
mixed rt c.time1 c.slope26, ///
           || personid: time1, variance reml covariance(unstructured) ///
           residuals(independent,t(session)), 
      estat ic, n(101), 
      estat recovariance, relevel(personid), 
      estat recovariance, relevel(personid) correlation, 
      estat wcorrelation, covariance, 
      estat wcorrelation, 
      * Intercept at Session=1 Time=0
      lincom _cons*1 + c.time1*0 + c.slope26*0 
      * Intercept at Session=2 Time=1
      lincom _cons*1 + c.time1*1 + c.slope26*0 
      * Intercept at Session=3 Time=2
      lincom _cons*1 + c.time1*2 + c.slope26*1 
      * Intercept at Session=4 Time=3
      lincom _cons*1 + c.time1*3 + c.slope26*2 
      * Intercept at Session=5 Time=4
      lincom _cons*1 + c.time1*4 + c.slope26*3 
      * Intercept at Session=6 Time=5
      lincom _cons*1 + c.time1*5 + c.slope26*4 
      * Rate of change from session 2 to 6
      lincom c.time1*1 + c.slope26*1 
      estimates store FitRand16Fix26, 
      lrtest FitRand16Fix26 FitFix16Fix26,  
 
display as result "Ch 6: 5c: Random Time, Random Slope26 Model" 
mixed rt c.time1 c.slope26, ///
           || personid: time1 slope26, variance reml covariance(unstructured) ///
           residuals(independent,t(session)), 
      estat ic, n(101), 
      estat recovariance, relevel(personid), 
      estat recovariance, relevel(personid) correlation, 
      estat wcorrelation, covariance, 
      estat wcorrelation, 
      * Intercept at Session=1 Time=0
      lincom _cons*1 + c.time1*0 + c.slope26*0 
      * Intercept at Session=2 Time=1
      lincom _cons*1 + c.time1*1 + c.slope26*0 
      * Intercept at Session=3 Time=2
      lincom _cons*1 + c.time1*2 + c.slope26*1 
      * Intercept at Session=4 Time=3
      lincom _cons*1 + c.time1*3 + c.slope26*2 
      * Intercept at Session=5 Time=4
      lincom _cons*1 + c.time1*4 + c.slope26*3 
      * Intercept at Session=6 Time=5
      lincom _cons*1 + c.time1*5 + c.slope26*4 
      * Rate of change from session 2 to 6
      lincom c.time1*1 + c.slope26*1 
      estimates store FitRand16Rand26, 
      lrtest FitRand16Rand26 FitRand16Fix26,  
 
display as result "Ch 6: 0: Saturated Means, Unstructured Variance Model" 
display as result "Using ML Instead of REML" 
mixed rt i.session, ///
           || personid: , noconstant variance mle covariance(unstructured) ///
           residuals(unstructured,t(session)), 
      estat ic, n(101), 
      estimates store FitMLSatUN, 
 
display as result "Ch 6: 1b: Empty Means, Random Intercept Model" 
display as result "Using ML Instead of REML" 
mixed rt , ///
           || personid: , variance mle covariance(unstructured) ///
           residuals(independent,t(session)), 
      estat ic, n(101), 
      estat icc, 
      estimates store FitMLEmptyRI, 
      lrtest FitMLSatUN FitMLEmptyRI,  
 
display as result "Ch 6: 3b: Random Quadratic Time Model" 
display as result "Using ML instead of REML" 
mixed rt c.time1 c.time1#c.time1, ///
           || personid: time1 time1sq, variance mle covariance(unstructured) ///
           residuals(independent,t(session)), 
      estat ic, n(101), 
      estimates store FitMLRandQuad, 
      lrtest FitMLRandQuad FitMLEmptyRI,  
      lrtest FitMLSatUN FitMLRandQuad,  
 
display as result "Ch 6: 4c: Random Slope12, Random Slope26 Model" 
display as result "Using ML instead of REML" 
mixed rt c.slope12 c.slope26, ///
           || personid: slope12 slope26, variance mle covariance(unstructured) ///
           residuals(independent,t(session)), 
      estat ic, n(101), 
      estimates store FitMLRand12Rand26, 
      lrtest FitMLRand12Rand26 FitMLEmptyRI,  
      lrtest FitMLSatUN FitMLRand12Rand26,  
  
****** END CHAPTER 6 MODELS ******
 
* Close log
log close STATA_Chapter6 
* Convert log to html using custom downloaded package
log2html $filesave\STATA_Chapter6_Output, replace