I am trying to plot the cohort group effect by country after estimating Hierarchical Age Period Cohort Model (HAPC):
This is the model - binary outcome variable with four levels (country, cohort group, period and individuals), cross-classified:
Code: Select all
sort cntry_n cohortgrp essround ind
runmlwin vote cons age_cen age_censq, ///
level4(cntry_n: cons) ///
level3(cohortgrp: cons) ///
level2(essround: cons) ///
level1(ind:) ///
discrete(distribution(binomial) link(logit) denom(cons) pql2)
display "`c(current_time)' `c(current_date)'"
runmlwin vote cons age_cen age_censq, ///
level4(cntry_n: cons, residuals(k)) ///
level3(cohortgrp: cons, residuals(v)) ///
level2(essround: cons, residuals(u)) ///
level1(ind:) ///
discrete(distribution(binomial) link(logit) denom(cons)) ///
mcmc(cc burnin(5000) chain(500000) hcen(3)) initsprevious nopause nogroup
display "`c(current_time)' `c(current_date)'"
save hapc, replace
estimates save hapc, replace
Code: Select all
egen pickone = tag(cohortgrp)
mat A = e(b)
scalar B = A[1,1]
gen newv0 = v0 + B
replace newv0 = invlogit(newv0)
gen newv0seHI = v0 + B + (1.96*v0se)
gen HIv0= invlogit(newv0seHI)
gen newv0seLO = v0 + B - (1.96*v0se)
gen LOv0= invlogit(newv0seLO)
sort cohort
twoway (line newv0 cohortgrp if pickone==1) ///
(line LOv0 cohortgrp if pickone==1, lpattern(dash)) ///
(line HIv0 cohortgrp if pickone==1, lpattern(dash)), ///
legend(off) scheme(s2mono) graphregion(color(white)) ///
ytitle("") xtitle("Birth Year") ///
xscale(range(1900 2000)) yscale(range(0.7 0.8)) ///
xlabel(1909(20)1990) ylabel(0.5(0.1)0.8)
graph save hacp_cohortgrp, replace
Afterwards, I try to plot by country simply adding the "by country, compact" to the twoway code:
Code: Select all
sort cohort
twoway (line newv0 cohortgrp if pickone==1, by(cntry, compact)) ///
(line LOv0 cohortgrp if pickone==1, lpattern(dash)) ///
(line HIv0 cohortgrp if pickone==1, lpattern(dash)), ///
legend(off) scheme(s2mono) graphregion(color(white)) ///
ytitle("") xtitle("Birth Year") ///
xscale(range(1900 2000)) yscale(range(0.7 0.8)) ///
xlabel(1909(20)1990) ylabel(0.5(0.1)0.8)
However, it shows only one country. Please share your suggestions.
Looking forward.
Regards,
Rza