use deleteme.dta
// Cross-level interaction model
qui mixed lsat c.wfc##c.gdp || country: wfc , cov(uns)
estimates store cross_lvl
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
esttab cross_lvl ///
, 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")) ///
coeflabel(wfc "Work-family conflict" ///
gdp "GDP" ///
c.wfc#c.gdp "Work-family conflict X GDP" ///
_cons "Intercept") ///
mtitles("Cross-level interaction") ///
nonumbers varwidth(28) modelwidth(24) ///
keep(lsat:) // Drop variance components in weird shapes
// Interaction plot (A)
preserve
estimates restore cross_lvl
qui sum wfc // get standard deviation and mean
local wfcminsd = r(mean) - r(sd)
local wfcmean = r(mean)
local wfcplusd = r(mean) + r(sd)
qui sum gdp // get standard deviation
local gdpmin2sd = r(mean) - 2*r(sd)
local gdpminsd = r(mean) - r(sd)
local gdpmean = r(mean)
local gdpplusd = r(mean) + r(sd)
local gdpplu2sd = r(mean) + 2*r(sd)
margins, at(gdp = (`gdpmin2sd' `gdpminsd' `gdpmean' `gdpplusd' `gdpplu2sd') ///
wfc = (`wfcminsd' `wfcmean' `wfcplusd')) vsquish
marginsplot, x(gdp) noci ///
xlabel(`gdpmin2sd' "-2 SD" ///
`gdpminsd' "-1 SD" ///
`gdpmean' `""Average" "GDP""' ///
`gdpplusd' "+1 SD" ///
`gdpplu2sd' "+2 SD") ///
xtitle("") ///
ylabel(, format(%6.1f)) ///
ytitle("Predicted happiness") ///
plotopts(msymbol(none)) /// // Turn off markers
plot1opts(lpattern(longdash)) /// // Define line types here
plot2opts(lpattern(solid)) ///
plot3opts(lpattern(shortdash)) ///
legend(subtitle("Work{c 150}family conflict" , size(small)) ///
order(1 "- 1 SD" 2 "Mean" 3 "+ 1 SD") size(small) ///
pos(11) ring(0)) ///
title("(A) Interaction plot") ///
name(xlvl_plot1, replace)
restore
// Interaction plot (B)
preserve
estimates restore cross_lvl
predict u1 u0, reffect
capture drop pickone
egen pickone = tag(country)
keep if pickone
gen wfceffect = u1 + _b[wfc] + _b[c.wfc#c.gdp] * gdp
qui sum gdp // get standard deviation
local gdpmin2sd = r(mean) - 2*r(sd)
local gdpminsd = r(mean) - r(sd)
local gdpmean = r(mean)
local gdpplusd = r(mean) + r(sd)
local gdpplu2sd = r(mean) + 2*r(sd)
twoway (scatter wfceffect gdp, mlabel(country) mlabpos(0) msymbol(none)) ///
(function y = _b[wfc] + _b[c.wfc#c.gdp] * x, range(gdp)) ///
, ytitle("Work{c 150}family conflict coefficient") ///
xlabel(`gdpmin2sd' "-2 SD" `gdpminsd' "-1 SD" ///
`gdpmean' `""Average" "GDP""' ///
`gdpplusd' "+1 SD" `gdpplu2sd' "+2 SD") ///
title("(B) Interaction plot") ///
legend(order(2 "Work{c 150}family conflict coefficient by GDP") ///
pos(7) ring(0) size(small)) ///
name(xlvl_plot2, replace)
restore
// Combine graphs
graph combine xlvl_plot1 xlvl_plot2, row(1) ysize(3) xsize(5.5) altshrink
