* Prevent "more" messages from appearing
set more off 
* Control line length
set linesize 150 
 
********************************************************************************
******* BEGIN DATA MANIPULATION FOR CHAPTER 10b THREE-LEVEL TIME EXAMPLE *******
*******               CHANGE "filesave" to your directory                *******
********************************************************************************
 
* Defining global variable for file location to be replaced in code below
global filesave "C:\Dropbox\PilesOfVariance\Chapter10b\STATA" 
 
* Import chapter 10b stacked data and create centered predictors for analysis
use "$filesave\STATA_Chapter10b.dta", clear 
* Piecewise slopes for session
gen slope12 = session  
gen slope26 = session  
recode slope12 (1=-1) if session==1  
recode slope12 (2=0)  if session==2  
recode slope12 (3=0)  if session==3  
recode slope12 (4=0)  if session==4  
recode slope12 (5=0)  if session==5  
recode slope12 (6=0)  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: Initial Slope (Session 1-2)"  
label variable slope26 "slope26: Later Slope (Session 2-6)" 
* Linear slope for burst.
gen burst1 = burst-1 
gen burst1sq = burst1*burst1 
* Difference between early and later bursts
gen b1or2=0 
recode b1or2 (0=1) if burst==1 
recode b1or2 (0=1) if burst==2 
label variable slope12 "slope12: Initial Slope (Session 1-2)" 
label variable slope26 "slope26: Later Slope (Session 2-6)" 
label variable burst1 "burst1: Measurement Burst (0=1)" 
label variable burst1sq "burst1sq: Quadratic Measurement Burst (0=1)" 
label variable b1or2 "b1or2: Is Burst 1 or 2 (0=no, 1=yes)" 
 
* Subset sample to complete cases for all predictors
egen nummiss = rowmiss(session burst) 
drop if nummiss>0 
     
********************************************************************************
*******          BEGIN CHAPTER 10b THREE-LEVEL TIME MODELS               *******
******* NOTE: No multivariate longitudinal models are possible in MIXED  *******
********************************************************************************
 
* Save results to separate file
log using $filesave\STATA_Chapter10b_Output, replace name(STATA_Chapter10b) 
 
display as result "Chapter 10b: Descriptive Statistics for Time-Varying Variables" 
summarize symptoms posaff 
 
display as result "Ch 10b: Empty Means, Single-Level Model for the Variance for Symptoms" 
display as result "Independent Observations" 
mixed symptoms , ///
                 || personid: , noconstant variance mle covariance(unstructured) ///
                 || burst: , noconstant covariance(unstructured), 
      estat ic, n(108), 
      estimates store FitEmpty1S, 
 
display as result "Ch 10b: Empty Means, Two-Level Model for the Variance for Symptoms" 
display as result "Sessions Within Burst*Persons" 
mixed symptoms , ///
                 || personid: , noconstant variance mle covariance(unstructured) ///
                 || burst: , covariance(unstructured), 
      estat ic, n(108), 
      estat wcorrelation, covariance, 
      estat wcorrelation, 
      estimates store FitEmpty2S, 
 
display as result "Eq 10b.5: Empty Means, Three-Level Model for the Variance for Symptoms" 
display as result "Level-1 Sessions Within Level-2 Bursts Within Level-3 Persons" 
mixed symptoms , ///
                 || personid: , variance mle covariance(unstructured) ///
                 || burst: , covariance(unstructured), 
      estat ic, n(108), 
      estat icc, 
      estat wcorrelation, covariance, 
      estat wcorrelation, 
      estimates store FitEmpty3S, 
      lrtest FitEmpty3S FitEmpty2S,  
 
display as result "Ch 10b: Empty Means, Single-Level Model for the Variance for Positive Affect" 
display as result "Independent Observations" 
mixed posaff , ///
               || personid: , noconstant variance mle covariance(unstructured) ///
               || burst: , noconstant covariance(unstructured), 
      estat ic, n(108), 
      estimates store FitEmpty1P, 
 
display as result "Ch 10b: Empty Means, Two-Level Model for the Variance for Positive Affect" 
display as result "Sessions Within Burst*Persons" 
mixed posaff , ///
               || personid: , noconstant variance mle covariance(unstructured) ///
               || burst: , covariance(unstructured), 
      estat ic, n(108), 
      estat wcorrelation, covariance, 
      estat wcorrelation, 
      estimates store FitEmpty2P, 
 
display as result "Eq 10b.5: Empty Means, Three-Level Model for the Variance for Positive Affect" 
display as result "Level-1 Sessions Within Level-2 Bursts Within Level-3 Persons" 
mixed posaff , ///
               || personid: , variance mle covariance(unstructured) ///
               || burst: , covariance(unstructured), 
      estat ic, n(108), 
      estat icc, 
      estat wcorrelation, covariance, 
      estat wcorrelation, 
      estimates store FitEmpty3P, 
      lrtest FitEmpty3P FitEmpty2P,  
 
display as result "Eq 10b.6: Saturated Means for Burst by Session" 
display as result "Three-Level Model for the Variance for Symptoms" 
mixed symptoms i.session i.burst i.session#i.burst, ///
                 || personid: , variance mle covariance(unstructured) ///
                 || burst: , covariance(unstructured), 
      estat ic, n(108), 
      contrast i.session,  
      margins  i.session,  
      contrast i.burst,  
      margins  i.burst,  
      contrast i.session#i.burst,  
      margins  i.session#i.burst,  
      margins  i.session@i.burst,  
      margins  i.burst@i.session,  
      estimates store FitSatAllS, 
 
display as result "Eq 10b.7: Piecewise Session Slopes by Observed Burst" 
display as result "Three-Level Model for the Variance for Symptoms" 
mixed symptoms c.slope12 c.slope26 i.burst ///
               c.slope12#i.burst c.slope26#i.burst, ///
                 || personid: , variance mle covariance(unstructured) ///
                 || burst: , covariance(unstructured), 
      estat ic, n(108), 
      * Session 2 at Burst 1
      lincom _cons*1 + i1.burst 
      * Session 2 at Burst 2
      lincom _cons*1 + i2.burst 
      * Session 2 at Burst 3
      lincom _cons*1 + i3.burst 
      * Session 2 at Burst 4
      lincom _cons*1 + i4.burst 
      * Session 2 at Burst 5
      lincom _cons*1 + i5.burst 
      * Slope12 at Burst 1
      lincom c.slope12*1 + c.slope12#i1.burst 
      * Slope12 at Burst 2
      lincom c.slope12*1 + c.slope12#i2.burst 
      * Slope12 at Burst 3
      lincom c.slope12*1 + c.slope12#i3.burst 
      * Slope12 at Burst 4
      lincom c.slope12*1 + c.slope12#i4.burst 
      * Slope12 at Burst 5
      lincom c.slope12*1 + c.slope12#i5.burst 
      * Slope26 at Burst 1
      lincom c.slope26*1 + c.slope26#i1.burst 
      * Slope26 at Burst 2
      lincom c.slope26*1 + c.slope26#i2.burst 
      * Slope26 at Burst 3
      lincom c.slope26*1 + c.slope26#i3.burst 
      * Slope26 at Burst 4
      lincom c.slope26*1 + c.slope26#i4.burst 
      * Slope26 at Burst 5
      lincom c.slope26*1 + c.slope26#i5.burst 
      estimates store FitPiecebyBurstMeansS, 
 
display as result "Eq 10b.8: Piecewise Session Slopes by Quadratic Burst" 
display as result "Three-Level Model for the Variance for Symptoms" 
mixed symptoms c.burst1 c.burst1#c.burst1 c.slope12 c.slope26 ///
               c.slope12#c.burst1 c.slope26#c.burst1 ///
               c.slope12#c.burst1#c.burst1 c.slope26#c.burst1#c.burst1, /// 
                 || personid: , variance mle covariance(unstructured) ///
                 || burst: , covariance(unstructured), 
      estat ic, n(108), 
      * Slope12 at Burst 1
      lincom c.slope12*1 + c.slope12#c.burst1*0 + c.slope12#c.burst1#c.burst1*0 
      * Slope12 at Burst 2
      lincom c.slope12*1 + c.slope12#c.burst1*1 + c.slope12#c.burst1#c.burst1*1 
      * Slope12 at Burst 3
      lincom c.slope12*1 + c.slope12#c.burst1*2 + c.slope12#c.burst1#c.burst1*4 
      * Slope12 at Burst 4
      lincom c.slope12*1 + c.slope12#c.burst1*3 + c.slope12#c.burst1#c.burst1*9 
      * Slope12 at Burst 5
      lincom c.slope12*1 + c.slope12#c.burst1*4 + c.slope12#c.burst1#c.burst1*16 
      * Slope26 at Burst 1
      lincom c.slope26*1 + c.slope26#c.burst1*0 + c.slope26#c.burst1#c.burst1*0 
      * Slope26 at Burst 2
      lincom c.slope26*1 + c.slope26#c.burst1*1 + c.slope26#c.burst1#c.burst1*1 
      * Slope26 at Burst 3
      lincom c.slope26*1 + c.slope26#c.burst1*2 + c.slope26#c.burst1#c.burst1*4 
      * Slope26 at Burst 4
      lincom c.slope26*1 + c.slope26#c.burst1*3 + c.slope26#c.burst1#c.burst1*9 
      * Slope26 at Burst 5
      lincom c.slope26*1 + c.slope26#c.burst1*4 + c.slope26#c.burst1#c.burst1*16 
      estimates store FitPiecebyQuadBurstS, 
      lrtest FitPiecebyBurstMeansS FitPiecebyQuadBurstS,  
 
display as result "Eq 10b.9: Piecewise Session Slopes (Slope26 by Quadratic Burst Only)" 
display as result "Three-Level Model for the Variance for Symptoms" 
mixed symptoms c.burst1 c.burst1#c.burst1 c.slope12 c.slope26 ///
               c.slope26#c.burst1 c.slope26#c.burst1#c.burst1, ///
                 || personid: , variance mle covariance(unstructured) ///
                 || burst: , covariance(unstructured), 
      estat ic, n(108), 
      * Slope26 at Burst 1
      lincom c.slope26*1 + c.slope26#c.burst1*0 + c.slope26#c.burst1#c.burst1*0 
      * Slope26 at Burst 2
      lincom c.slope26*1 + c.slope26#c.burst1*1 + c.slope26#c.burst1#c.burst1*1 
      * Slope26 at Burst 3
      lincom c.slope26*1 + c.slope26#c.burst1*2 + c.slope26#c.burst1#c.burst1*4 
      * Slope26 at Burst 4
      lincom c.slope26*1 + c.slope26#c.burst1*3 + c.slope26#c.burst1#c.burst1*9 
      * Slope26 at Burst 5
      lincom c.slope26*1 + c.slope26#c.burst1*4 + c.slope26#c.burst1#c.burst1*16 
      estimates store FitPieceQuadBurst26S, 
      lrtest FitPiecebyBurstMeansS FitPieceQuadBurst26S,  
      lrtest FitSatAllS FitPieceQuadBurst26S,  
 
display as result "Ch 10b: Piecewise Session Slopes (Slope26 by Quadratic Burst Only)" 
display as result "Add Random Linear Burst Slope across Persons for Symptoms" 
mixed symptoms c.burst1 c.burst1#c.burst1 c.slope12 c.slope26 ///
               c.slope26#c.burst1 c.slope26#c.burst1#c.burst1, ///
                 || personid: burst1, variance mle covariance(unstructured) ///
                 || burst: , covariance(unstructured), 
      estat ic, n(108), 
      estat recovariance, relevel(personid), 
      estat recovariance, relevel(personid) correlation, 
      estimates store FitRandBurstLin3S, 
      lrtest FitRandBurstLin3S FitPieceQuadBurst26S,  
 
display as result "Ch 10b: Piecewise Session Slopes (Slope26 by Quadratic Burst Only)" 
display as result "Add Random Quadratic Burst Slope across Persons for Symptoms" 
mixed symptoms c.burst1 c.burst1#c.burst1 c.slope12 c.slope26 ///
               c.slope26#c.burst1 c.slope26#c.burst1#c.burst1, ///
                 || personid: burst1 burst1sq, variance mle covariance(unstructured) ///
                 || burst: , covariance(unstructured), 
      estat ic, n(108), 
      estat recovariance, relevel(personid), 
      estat recovariance, relevel(personid) correlation, 
      estimates store FitRandBurstQuad3S, 
      lrtest FitRandBurstQuad3S FitRandBurstLin3S,  
 
display as result "Ch 10b: Piecewise Session Slopes (Slope26 by Quadratic Burst Only)" 
display as result "Add Random Linear Slope12 Across Bursts for Symptoms" 
mixed symptoms c.burst1 c.burst1#c.burst1 c.slope12 c.slope26 ///
               c.slope26#c.burst1 c.slope26#c.burst1#c.burst1, ///
                 || personid: burst1, variance mle covariance(unstructured) ///
                 || burst: slope12, covariance(unstructured), 
      estat ic, n(108), 
      estat recovariance, relevel(personid), 
      estat recovariance, relevel(personid) correlation, 
      estat recovariance, relevel(burst), 
      estat recovariance, relevel(burst) correlation, 
      estimates store FitRandSlope12at2S, 
      lrtest FitRandSlope12at2S FitRandBurstLin3S,  
 
display as result "Ch 10b: Piecewise Session Slopes (Slope26 by Quadratic Burst Only)" 
display as result "Add Random Linear Slope12 Across Persons for Symptoms" 
mixed symptoms c.burst1 c.burst1#c.burst1 c.slope12 c.slope26 ///
               c.slope26#c.burst1 c.slope26#c.burst1#c.burst1, ///
                 || personid: burst1 slope12, variance mle covariance(unstructured) ///
                 || burst: slope12, covariance(unstructured), 
      estat ic, n(108), 
      estat recovariance, relevel(personid), 
      estat recovariance, relevel(personid) correlation, 
      estat recovariance, relevel(burst), 
      estat recovariance, relevel(burst) correlation, 
      estimates store FitRandSlope12at23S, 
      lrtest FitRandSlope12at23S FitRandSlope12at2S, force 
 
display as result "Eq 10b.10: Piecewise Session Slopes (Slope26 by Quadratic Burst Only)" 
display as result "Add Random Linear Slope26 Across Bursts for Symptoms" 
mixed symptoms c.burst1 c.burst1#c.burst1 c.slope12 c.slope26 ///
               c.slope26#c.burst1 c.slope26#c.burst1#c.burst1, ///
                 || personid: burst1, variance mle covariance(unstructured) ///
                 || burst: slope12 slope26, covariance(unstructured), 
      estat ic, n(108), 
      estat recovariance, relevel(personid), 
      estat recovariance, relevel(personid) correlation, 
      estat recovariance, relevel(burst), 
      estat recovariance, relevel(burst) correlation, 
      * Slope26 at Burst 1
      lincom c.slope26*1 + c.slope26#c.burst1*0 + c.slope26#c.burst1#c.burst1*0 
      * Slope26 at Burst 2
      lincom c.slope26*1 + c.slope26#c.burst1*1 + c.slope26#c.burst1#c.burst1*1 
      * Slope26 at Burst 3
      lincom c.slope26*1 + c.slope26#c.burst1*2 + c.slope26#c.burst1#c.burst1*4 
      * Slope26 at Burst 4
      lincom c.slope26*1 + c.slope26#c.burst1*3 + c.slope26#c.burst1#c.burst1*9 
      * Slope26 at Burst 5
      lincom c.slope26*1 + c.slope26#c.burst1*4 + c.slope26#c.burst1#c.burst1*16 
      estimates store FitRandSlope26at2S, 
      lrtest FitRandSlope26at2S FitRandSlope12at2S,  
      predict PredUncS, xb, 
      margins, at (c.slope12=-1 c.slope26=0 c.burst1=(0(1)4)) vsquish, 
      margins, at (c.slope12=0 c.slope26=(0(1)4) c.burst1=(0(1)4)) vsquish, 
corr symptoms PredUncS 
 
display as result "Ch 10b: Piecewise Session Slopes (Slope26 by Quadratic Burst Only)" 
display as result "Add Random Linear Slope26 Across Persons for Symptoms" 
mixed symptoms c.burst1 c.burst1#c.burst1 c.slope12 c.slope26 ///
               c.slope26#c.burst1 c.slope26#c.burst1#c.burst1, ///
                 || personid: burst1 slope26, variance mle covariance(unstructured) ///
                 || burst: slope12 slope26, covariance(unstructured), 
      estat ic, n(108), 
      estat recovariance, relevel(personid), 
      estat recovariance, relevel(personid) correlation, 
      estat recovariance, relevel(burst), 
      estat recovariance, relevel(burst) correlation, 
      estimates store FitRandSlope26at23S, 
      lrtest FitRandSlope26at23S FitRandSlope26at2S,  
 
display as result "Ch 10b: Final Unconditional Model for Symptoms" 
display as result "Remove Level-2 Random Effects Variances and Covariances" 
mixed symptoms c.burst1 c.burst1#c.burst1 c.slope12 c.slope26 ///
               c.slope26#c.burst1 c.slope26#c.burst1#c.burst1, ///
                 || personid: burst1, variance mle covariance(unstructured), 
      estat ic, n(108), 
      estat recovariance, relevel(personid), 
      estat recovariance, relevel(personid) correlation, 
      estimates store FitNo2S, 
      lrtest FitRandSlope26at2S FitNo2S,  
 
display as result "Eq 10b.6: Saturated Means for Burst by Session" 
display as result "Three-Level Model for the Variance for Positive Affect" 
mixed posaff i.session i.burst i.session#i.burst, ///
               || personid: , variance mle covariance(unstructured) ///
               || burst: , covariance(unstructured), 
      estat ic, n(108), 
      contrast i.session,  
      margins  i.session,  
      contrast i.burst,  
      margins  i.burst,  
      contrast i.session#i.burst,  
      margins  i.session#i.burst,  
      margins  i.session@i.burst,  
      margins  i.burst@i.session,  
      estimates store FitSatAllP, 
 
display as result "Eq 10b.7: Piecewise Session Slopes by Observed Burst" 
display as result "Three-Level Model for the Variance for Positive Affect" 
mixed posaff c.slope12 c.slope26 i.burst ///
             c.slope12#i.burst c.slope26#i.burst, ///
               || personid: , variance mle covariance(unstructured) ///
               || burst: , covariance(unstructured), 
      estat ic, n(108), 
      * Session 2 at Burst 1
      lincom _cons*1 + i1.burst 
      * Session 2 at Burst 2
      lincom _cons*1 + i2.burst 
      * Session 2 at Burst 3
      lincom _cons*1 + i3.burst 
      * Session 2 at Burst 4
      lincom _cons*1 + i4.burst 
      * Session 2 at Burst 5
      lincom _cons*1 + i5.burst 
      * Slope12 at Burst 1
      lincom c.slope12*1 + c.slope12#i1.burst 
      * Slope12 at Burst 2
      lincom c.slope12*1 + c.slope12#i2.burst 
      * Slope12 at Burst 3
      lincom c.slope12*1 + c.slope12#i3.burst 
      * Slope12 at Burst 4
      lincom c.slope12*1 + c.slope12#i4.burst 
      * Slope12 at Burst 5
      lincom c.slope12*1 + c.slope12#i5.burst 
      * Slope26 at Burst 1
      lincom c.slope26*1 + c.slope26#i1.burst 
      * Slope26 at Burst 2
      lincom c.slope26*1 + c.slope26#i2.burst 
      * Slope26 at Burst 3
      lincom c.slope26*1 + c.slope26#i3.burst 
      * Slope26 at Burst 4
      lincom c.slope26*1 + c.slope26#i4.burst 
      * Slope26 at Burst 5
      lincom c.slope26*1 + c.slope26#i5.burst 
      estimates store FitPiecebyBurstMeansP, 
      lrtest FitSatAllP FitPiecebyBurstMeansP,  
 
display as result "Eq 10b.11and12: Piecewise Session Slopes, Linear Burst, Slope12 by Burst1or2" 
display as result "Three-Level Model for the Variance for Positive Affect" 
mixed posaff c.burst1 c.slope12 c.slope12#c.b1or2 c.slope26, ///
               || personid: , variance mle covariance(unstructured) ///
               || burst: , covariance(unstructured), 
      estat ic, n(108), 
      lincom c.slope12*1 + c.slope12#c.b1or2*1 // Slope12 at Burst 1 or 2
      estimates store FitPieceBurst1or2P, 
      lrtest FitPiecebyBurstMeansP FitPieceBurst1or2P,  
      lrtest FitSatAllP FitPieceBurst1or2P,  
 
display as result "Ch 10b: Piecewise Session Slopes, Linear Burst, Slope12 by Burst1or2" 
display as result "Add Random Linear Burst across Persons for Positive Affect" 
mixed posaff c.burst1 c.slope12 c.slope12#c.b1or2 c.slope26, ///
               || personid: burst1, variance mle covariance(unstructured) ///
               || burst: , covariance(unstructured), 
      estat ic, n(108), 
      estat recovariance, relevel(personid), 
      estat recovariance, relevel(personid) correlation, 
      estimates store FitRandBurstLin3P, 
      lrtest FitRandBurstLin3P FitPieceBurst1or2P,  
 
display as result "Ch 10b: Piecewise Session Slopes, Linear Burst, Slope12 by Burst1or2" 
display as result "Add Fixed and Random Quadratic Burst across Persons for Positive Affect" 
mixed posaff c.burst1 c.slope12 c.slope12#c.b1or2 c.slope26 ///
             c.burst1#c.burst1, ///
               || personid: burst1 burst1sq, variance mle covariance(unstructured) ///
               || burst: , covariance(unstructured), 
      estat ic, n(108), 
      estat recovariance, relevel(personid), 
      estat recovariance, relevel(personid) correlation, 
      estimates store FitRandBurstQuad3P, 
      lrtest FitRandBurstQuad3P FitRandBurstLin3P,  
 
display as result "Ch 10b: Piecewise Session Slopes, Linear Burst, Slope12 by Burst1or2" 
display as result "Add Random Linear Slope12 across Bursts for Positive Affect" 
mixed posaff c.burst1 c.slope12 c.slope12#c.b1or2 c.slope26, ///
               || personid: burst1, variance mle covariance(unstructured) ///
               || burst: slope12, covariance(unstructured), 
      estat ic, n(108), 
      estat recovariance, relevel(personid), 
      estat recovariance, relevel(personid) correlation, 
      estat recovariance, relevel(burst), 
      estat recovariance, relevel(burst) correlation, 
      estimates store FitRandSlope12at2P, 
      lrtest FitRandSlope12at2P FitRandBurstLin3P,  
 
display as result "Ch 10b: Piecewise Session Slopes, Linear Burst, Slope12 by Burst1or2" 
display as result "Add Random Linear Slope12 across Persons for Positive Affect" 
mixed posaff c.burst1 c.slope12 c.slope12#c.b1or2 c.slope26, ///
               || personid: burst1 slope12, variance mle covariance(unstructured) ///
               || burst: slope12, covariance(unstructured), 
      estat ic, n(108), 
      estat recovariance, relevel(personid), 
      estat recovariance, relevel(personid) correlation, 
      estat recovariance, relevel(burst), 
      estat recovariance, relevel(burst) correlation, 
      estimates store FitRandSlope12at23P, 
      lrtest FitRandSlope12at23P FitRandSlope12at2P, force 
 
display as result "Ch 10b: Piecewise Session Slopes, Linear Burst, Slope12 by Burst1or2" 
display as result "Add Random Linear Slope26 across Bursts for Positive Affect" 
mixed posaff c.burst1 c.slope12 c.slope12#c.b1or2 c.slope26, ///
               || personid: burst1 slope12, variance mle covariance(unstructured) ///
               || burst: slope12 slope26, covariance(unstructured), 
      estat ic, n(108), 
      estat recovariance, relevel(personid), 
      estat recovariance, relevel(personid) correlation, 
      estat recovariance, relevel(burst), 
      estat recovariance, relevel(burst) correlation, 
      estimates store FitRandSlope26at2P, 
      lrtest FitRandSlope26at2P FitRandSlope12at23P,  
 
display as result "Eq 10b.13: Piecewise Session Slopes, Linear Burst, Slope12 by Burst1or2" 
display as result "Add Random Linear Slope26 across Persons for Positive Affect" 
mixed posaff c.burst1 c.slope12 c.slope12#c.b1or2 c.slope26, ///
               || personid: burst1 slope12 slope26, variance mle covariance(unstructured) ///
               || burst: slope12 slope26, covariance(unstructured), 
      estat ic, n(108), 
      estat recovariance, relevel(personid), 
      estat recovariance, relevel(personid) correlation, 
      estat recovariance, relevel(burst), 
      estat recovariance, relevel(burst) correlation, 
      lincom c.slope12*1 + c.slope12#c.b1or2*1 // Slope12 at Burst 1 or 2
      estimates store FitRandSlope26at23P, 
      lrtest FitRandSlope26at23P FitRandSlope26at2P,  
      predict PredUncP, xb, 
      margins, at (c.slope12=-1 c.slope26=0 c.b1or2=1 c.burst1=(0(1)1)) vsquish, 
      margins, at (c.slope12=-1 c.slope26=0 c.b1or2=0 c.burst1=(2(1)4)) vsquish, 
      margins, at (c.slope12=0 c.slope26=(0(1)4) c.b1or2=1 c.burst1=(0(1)1)) vsquish, 
      margins, at (c.slope12=0 c.slope26=(0(1)4) c.b1or2=0 c.burst1=(2(1)4)) vsquish, 
corr posaff PredUncP 
 
display as result "Ch 10b: Final Unconditional Model for Positive Affect" 
display as result "Removing Level-2 Random Effects Variances and Covariances" 
mixed posaff c.burst1 c.slope12 c.slope12#c.b1or2 c.slope26, ///
               || personid: burst1 slope12 slope26, variance mle covariance(unstructured), 
      estat ic, n(108), 
      estat recovariance, relevel(personid), 
      estat recovariance, relevel(personid) correlation, 
      estimates store FitNo2P, 
      lrtest FitRandSlope26at23P FitNo2P,  
  
****** END CHAPTER 10b MODELS ******
 
* Close log
log close STATA_Chapter10b 
* Convert log to html using custom downloaded package
log2html $filesave\STATA_Chapter10b_Output, replace