// Fit a couple of models
// Model 1
eststo clear
eststo: xtmixed wellbeing || cluster: , ml var
qui estadd scalar dev = -2*e(ll) // Deviance
qui matrix foo = e(N_g) // Number of level 2 units
qui estadd scalar nc = foo[1,1] // Number of level 2 units
qui estadd scalar v1 = exp(2*[lns1_1_1]_b[_cons]) // Intercept variance
qui estadd scalar v_e = exp(2*[lnsig_e]_b[_cons]) // Residual variance
// Model 2
eststo: xtmixed wellbeing hourscgm || cluster: hourscgm, ml var cov(uns)
qui estadd scalar dev = -2*e(ll) // Deviance
qui matrix foo = e(N_g) // Number of level 2 units
qui estadd scalar nc = foo[1,1] // Number of level 2 units
qui estadd scalar v2 = exp(2*[lns1_1_1]_b[_cons]) // Slope variance
qui estadd scalar v1 = exp(2*[lns1_1_2]_b[_cons]) // Intercept variance
qui estadd scalar cov = tanh([atr1_1_1_2]_b[_cons]) * /// Slope-intercept covariance
exp([lns1_1_1]_b[_cons]) * ///
exp([lns1_1_2]_b[_cons])
qui estadd scalar v_e = exp(2*[lnsig_e]_b[_cons]) // Residual variance
// Model 3
eststo: xtmixed wellbeing hourscwc || cluster: hourscwc, ml var cov(uns)
qui estadd scalar dev = -2*e(ll) // Deviance
qui matrix foo = e(N_g) // Number of level 2 units
qui estadd scalar nc = foo[1,1] // Number of level 2 units
qui estadd scalar v2 = exp(2*[lns1_1_1]_b[_cons]) // Slope variance
qui estadd scalar v1 = exp(2*[lns1_1_2]_b[_cons]) // Intercept variance
qui estadd scalar cov = tanh([atr1_1_1_2]_b[_cons]) * /// Slope-intercept covariance
exp([lns1_1_1]_b[_cons]) * ///
exp([lns1_1_2]_b[_cons])
qui estadd scalar v_e = exp(2*[lnsig_e]_b[_cons]) // Residual variance
// Model 4
eststo: xtmixed wellbeing sizecgm hourscgm || cluster: hourscgm, ml var cov(uns)
qui estadd scalar dev = -2*e(ll) // Deviance
qui matrix foo = e(N_g) // Number of level 2 units
qui estadd scalar nc = foo[1,1] // Number of level 2 units
qui estadd scalar v2 = exp(2*[lns1_1_1]_b[_cons]) // Slope variance
qui estadd scalar v1 = exp(2*[lns1_1_2]_b[_cons]) // Intercept variance
qui estadd scalar cov = tanh([atr1_1_1_2]_b[_cons]) * /// Slope-intercept covariance
exp([lns1_1_1]_b[_cons]) * ///
exp([lns1_1_2]_b[_cons])
qui estadd scalar v_e = exp(2*[lnsig_e]_b[_cons]) // Residual variance
// Model 5
eststo: xtmixed wellbeing sizecgm hourscwc || cluster: hourscwc, ml var cov(uns)
qui estadd scalar dev = -2*e(ll) // Deviance
qui matrix foo = e(N_g) // Number of level 2 units
qui estadd scalar nc = foo[1,1] // Number of level 2 units
qui estadd scalar v2 = exp(2*[lns1_1_1]_b[_cons]) // Slope variance
qui estadd scalar v1 = exp(2*[lns1_1_2]_b[_cons]) // Intercept variance
qui estadd scalar cov = tanh([atr1_1_1_2]_b[_cons]) * /// Slope-intercept covariance
exp([lns1_1_1]_b[_cons]) * ///
exp([lns1_1_2]_b[_cons])
qui estadd scalar v_e = exp(2*[lnsig_e]_b[_cons]) // Residual variance
// Create table via -esttab-

esttab est1 est2 est3 est4 est5 ///
, se ///
stats(v1 v2 v_e cov dev nc N, /// Add variance components to table
labels("Var(Intercept)" ///
"Var(Slope)" ///
"Var(Residual)" ///
"Cov(Int., Slope)" ///
"Deviance" ///
"No. clusters" ///
"No. individuals")) ///
label /// Use variable labels
keep(wellbeing:) // Drop variance components in weird shapes

// Give it another go
esttab est1 est2 est4 est3 est5 ///
, se ///
stats(v1 v2 v_e cov dev nc N, ///
labels("Var(Intercept)" ///
"Var(Slope)" ///
"Var(Residual)" ///
"Cov(Int., Slope)" ///
"Deviance" ///
"No. clusters" ///
"No. individuals")) ///
rename(hourscwc hours hourscgm hours) /// Match up both hour variables
varlabels(hours "Work hours (CGM/CWC)" /// Relabel variables
sizecgm "Workgroup size (CGM)" ///
_cons "Constant") ///
wrap /// Wrap variable labels
mgroups("Null model" ///
"Work hours CGM" /// Group models
"Work hours CWC", pattern(1 1 0 1 0) span) ///
eqlabels("") /// Suppress equation label
nomtitle /// Suppress model titles
keep(wellbeing:) // Drop weird variance components