use cntry happy mnactic agea gndr using ESS1e06_4.dta, clear
// Prepare variables
recode agea (999 = .), gen(age)
recode gndr (1 = 0 "Male") (2 = 1 "Female") (9 = .), gen(female)
recode mnactic (77 88 99 = .), gen(active)
label var female "Sex (ref. Male)"
label var age "Age in decades"
label var happy "Happiness"
label var active "Labor market status"
label val active mnactic
center happy age // Center outcome to make intercept smaller for plot
replace c_age = c_age / 10
// Listwise deletion
drop if missing(female, c_age, c_happy, active, cntry)
drop happy gndr agea mnactic // Unecessary variables can go
// Fit model
mixed c_happy c_age i.female i.active || cntry: c_age, ml cov(uns)
// Save random part
local var_age = round(exp(_b[lns1_1_1:_cons])^2, .0001)
local var_int = round(exp(_b[lns1_1_2:_cons])^2, .0001)
local covaria = round(tanh(_b[atr1_1_1_2:_cons]) * ///
exp(_b[lns1_1_1:_cons]) * ///
exp(_b[lns1_1_2:_cons]), .0001)
// Plot fixed part
coefplot, xline(0) ///
xtitle(" " "Estimates and 95% CI's") ///
scheme(s1mono) ///
ciopts(recast(rcap)) ///
coeflabels(_cons = "{bf:Intercept}" ///
c_age = "{bf:Age} in decades" ///
8.active = "Housework") /// // Shorten label
headings(c_age = " " ///
1.female = "{bf:Sex} ({it:ref.} Male)" ///
2.active = `""{bf:Labor market status}" "({it:ref.} Paid work)""' ///
_cons = " ") ///
mlabel format(%9.2f) mlabposition(12) mlabgap(*1.5) ///
xscale(alt) msymbol(x) ///
ysize(8) xsize(4) ///
title("Fixed part", span) ///
name(fixed_part, replace) nodraw
/* The default sizes of the available area are -ysize(4)- and -xsize(5.5)-,
by the way. Letter size is -ysize(11)- and -xsize(8.5)-*/
// Plot random part
// Estimate residuals
capture drop u1* u0*
predict u1 u0, reffects
predict u1se u0se, reses
// Plot intercept variation
preserve
egen pickone = tag(cntry)
keep if pickone
egen order_ = rank(-u0), unique
labmask order_, value(cntry)
gen high = u0 + (1.96 * u0se)
gen low = u0 - (1.96 * u0se)
twoway (rcap u0 u0 order_, dsymbol(x)) ///
(rspike high low order_) , ///
yline(0) xlabel(1/22, val ang(v)) ///
title("Intercept variance = `var_int'") ///
xtitle("") ///
ytitle("Random intercept residuals") ///
legend(off) ///
name(rand_int, replace) nodraw
restore
// Plot slope variation
preserve
egen pickone = tag(cntry)
keep if pickone
egen order_ = rank(-u1), unique
labmask order_, value(cntry)
gen high = u1 + (1.96 * u1se)
gen low = u1 - (1.96 * u1se)
twoway (rcap u1 u1 order_, dsymbol(x)) ///
(rspike high low order_) , ///
yline(0) xlabel(1/22, val ang(v)) ///
xtitle("") ///
title("Slope variance = `var_age'") ///
ytitle("Random slope residuals") ///
ylabel(-1 (.5) 1) /// // Make y-axis identical to other plot,
legend(off) /// // otherwise x-axis looks different, too
name(rand_slope, replace) nodraw
restore
// Plot intercept-slope covariance
preserve
gen predRandomSlope= (_b[_cons] + u0) + ((_b[c_age] + u1) * c_age)
qui sum c_age
local hi1 = 1*r(sd)
local hi2 = 2*r(sd)
local hi3 = 3*r(sd)
local lo1 = -1*r(sd)
local lo2 = -2*r(sd)
sort cntry c_age
twoway (line predRandomSlope c_age, connect(ascending)), ///
ytitle("Predicted Happiness") ///
xtitle("") ///
title("Intercept{c 150}slope covariance = `covaria'") ///
xlabel(`lo1' "-1 SD" 0 "Average age" `hi1' "+1 SD" `hi2' "+2 SD" `hi3' "+3 SD") ///
ylabel(,format(%6.1f)) ///
name(covariance, replace) nodraw
restore
// Combine Figures
graph combine rand_int rand_slope covariance, col(1) ysize(8) title("Random part") name(random_part, replace) nodraw
graph combine fixed_part random_part, col(2) ysize(8) xsize(8) altshrink title("Random coefficient model")
