/******************************************************************************* The covid dataset *******************************************************************************/ use covid, clear describe *****CWM K=2 cwmglm y x1 x2 x3 n1 female, xnormal(x1 x2 x3) xpoisson(n1) xbin(female) k(2) /// posterior(z2z) family(poisson) matlist e(mu), title("mean of multivariate Gaussian covariates ") matlist e(sigma),title("variance-covariance matrix of multivariate Gaussian covariates") matlist e(lambda), title("mean of Poisson covariates") matlist e(p_binomial), title("mean of binomial covariates") matrix results= ( e(ic)' \ e(R2)[.,"Overall"] \ e(globaldeviance)["Normalized Deviance", 2..3]') *****CWM K=3 cwmglm y x1 x2 x3 n1 female, xnormal(x1 x2 x3) xpoisson(n1) xbin(female) k(3) /// posterior(z3z) family(poisson) matlist e(mu), title("mean of multivariate Gaussian covariates ") matlist e(sigma),title("variance-covariance matrix of multivariate Gaussian covariates") matlist e(lambda), title("mean of Poisson covariates") matlist e(p_binomial), title("mean of binomial covariates") matrix results= results, ( e(ic)' \ e(R2)[.,"Overall"] \ e(globaldeviance)["Normalized Deviance", 2..3]') *****CWM K=4 cwmglm y x1 x2 x3 n1 female, xnormal(x1 x2 x3) xpoisson(n1) xbin(female) k(4) /// posterior(z4z) family(poisson) matlist e(mu), title("mean of multivariate Gaussian covariates ") matlist e(sigma),title("variance-covariance matrix of multivariate Gaussian covariates") matlist e(lambda), title("mean of Poisson covariates") matlist e(p_binomial), title("mean of binomial covariates") matrix results= results, ( e(ic)' \ e(R2)[.,"Overall"] \ e(globaldeviance)["Normalized Deviance", 2..3]') ***summary of model fit matrix colnames results= "k=2" "k=3" "k=4" matlist results, title("summary of model fit") /******************************************************************************* The students dataset *******************************************************************************/ use students, clear cwmglm weight height heightf, k(2) posterior(z) xnormal(height heightf) eee matlist e(mu) matlist e(sigma) cap drop group predict group tab group gender tw (scatter height weight if group==1 & gender=="M", mcolor(red) msymbol(T)) /// (scatter height weight if group==2 & gender=="M", mcolor(blue) msymbol(T)) /// (scatter height weight if group==1 & gender=="F", mcolor(green) msymbol(O)) /// (scatter height weight if group==2 & gender=="F", mcolor(red) msymbol(O)), /// legend(off) saving(g1,replace) tw (scatter height heightf if group==1 & gender=="M", mcolor(red) msymbol(T)) /// (scatter height heightf if group==2 & gender=="M", mcolor(blue) msymbol(T)) /// (scatter height heightf if group==1 & gender=="F", mcolor(green) msymbol(O)) /// (scatter height heightf if group==2 & gender=="F", mcolor(red) msymbol(O)), /// legend(off) saving(g2,replace) tw (scatter heightf weight if group==1 & gender=="M", mcolor(red) msymbol(T)) /// (scatter heightf weight if group==2 & gender=="M", mcolor(blue) msymbol(T)) /// (scatter heightf weight if group==1 & gender=="F", mcolor(green) msymbol(O)) /// (scatter heightf weight if group==2 & gender=="F", mcolor(red) msymbol(O)), /// legend(off) saving(g3,replace) graph combine g1.gph g2.gph g3.gph, rows(1) cap drop z1 z2 cwmglm w height heightf, k(2) posterior(z) xnormal(height heightf) eee cwmbootstrap, nreps(100) **see the returned results matlist r(b) matlist r(mu) /******************************************************************************* The multinorm dataset *******************************************************************************/ use multinorm, clear tw (scatter x1 x2 if group==1) /// (scatter x1 x2 if group==2) /// (scatter x1 x2 if group==3), /// legend(order (1 "comp. 1" 2 "comp. 2" 3 "comp. 3")) /// legend( rows(1)) title(Artificial Data) local models vev evv vvv eei vei evi vvi eii vii eee vee eve vve eev local bestbic=10e20 local bestaic=10e20 cap matrix drop res foreach model of local models { forval i=2/5 { cap drop _tau* qui cwmglm, xnorm(x1 x2) k(`i') posterior(_tau) `model' if (e(converged)==1) { matrix ic=(e(ic),`i', e(ll)) matrix rownames ic= "`model'" matrix res = nullmat(res) \ ic local current_BIC=e(ic)[1,2] if (`current_BIC'<`bestbic') { local bestbic=`current_BIC' local bestk_BIC=`i' local bestmodel_BIC `model' } local current_AIC=e(ic)[1,1] if (`current_AIC'<`bestaic') { local bestaic=`current_AIC' local bestk_AIC=`i' local bestmodel_AIC `model' } } else di in red /// "model `model' with `i' mixture component did not converge" } } di as result "best model according to BIC: k=`bestk_BIC' type `bestmodel_BIC'" di as result "best model according to AIC: k=`bestk_AIC' type `bestmodel_AIC'" ***res is the matrix that collects all the cwm that have converged matrix colnames res=AIC BIC k ll matlist res ***** re-estimating selected model cap drop _tau* cwmglm, xnorm(x1 x2) k(`bestk_BIC') posterior(_tau) `bestmodel_BIC' matlist e(mu) matlist e(sigma) predict map tw (hist x1) (kdensity x1 [aw=_tau1], lcolor(navy)) (kdensity x1 [aw=_tau2], /// lcolor(maroon)) (kdensity x1 [aw=_tau3], lcolor(forest_green)), /// legend(rows(1)) /// legend(order(1 "Observed PDF" 2 "comp.1" 3 "comp.2" 4 "comp.3")) /// saving(gg1,replace) title(x1) tw (hist x2) (kdensity x2 [aw=_tau1], lcolor(navy)) (kdensity x2 [aw=_tau2], /// lcolor(maroon)) (kdensity x2 [aw=_tau3], lcolor(forest_green)), /// legend(rows(1)) /// legend(order(1 "Observed PDF" 2 "comp.1" 3 "comp.2" 4 "comp.3")) /// saving(gg2,replace) title(x2) tab map group tw (scatter x1 x2 if map==1) (scatter x1 x2 if map==2) /// (scatter x1 x2 if map==3), /// legend(order (1 "CWM Comp. 1" 2 "CWM Comp. 2" 3 "CWM Comp. 3")) /// legend( rows(1)) title(Artificial Data) subtitle(Estimated components) graph combine g1.gph g2.gph g3.gph, rows(1) ycommon /******************************************************************************* The gsem_mixture dataset *******************************************************************************/ webuse gsem_mixture, clear set seed 435162 ***supplying the same starting values to cwmgl and fmm gen start=rbinomial(1,.5) replace start=1+start quietly tab start, gen(class) cwmglm drvisits private medicaid c.age##c.age actlim chronic, /// family(poisson) k(2) posterior(tau) start(custom) initial(class1 class2) predict map_cwm fmm 2, emopts(iterate(500)) startvalues(classid start): /// poisson drvisits private medicaid c.age##c.age actlim chronic predict z1, classposteriorpr class(1) predict z2, classposteriorpr class(2) estat lcprob estat ic cap drop map_fmm corr z1 tau1 tw (scatter z1 tau1) (function y=x), /// ytitle(CWM posterior probability) xtitle(FMM posterior probability) /// title(class 1) legend(off) saving(g1,replace) subtitle({&rho}=0.992) tw (scatter z2 tau2) (function y=x), /// ytitle(CWM posterior probability) xtitle("FMM posterior probability") /// title(class 2) legend(off) saving(g2,replace) subtitle({&rho}=0.992) gen map_fmm=cond(z1>z2,1,2) tab map_fmm map_cwm graph combine g1.gph g2.gph