mixreg.gss

gauss 离散型计量估计源代码,直接下载下来就可以使用

	else; @ no observations, no updating of prior @
		fnk   = f0;
		gnki  = g0i;
	endif;
	gnk   = invpd(gnki);
	{lambdaik, lambdak} 	= wishart(rankx,fnk,gnk);

	dlambik  = det(lambdaik);          @ determinant of lambdaik @
	@ store lambda's  @
	lambdai[iptgp[k,1]:iptgp[k,2],.] 	= lambdaik;
	lambda[iptgp[k,1]:iptgp[k,2],.]  	= lambdak;
	bvars[iptgp[k,1]:iptgp[k,2]]		= diag(lambdak); @ Save diagonals for plotting @

	@ Compute the probability of class k for all subjects @
	if mmod > 1;
		@ Start computing P(z[i] = k) @
		zprob[.,k] = 0.5*ln(dlambik)*ones(nsub,1);
		resid	= beta - thetak;
		resid2	= lambdaik*resid;
		for i0 (1,nsub,1); i = i0;
			zprob[i,k] = zprob[i,k] - 0.5*resid[.,i]'resid2[.,i];
		endfor;	
	endif;
endfor;
 
@ Normalize zprob to avoid overflow @
zmaxp = maxc(zprob');
zprob = exp(zprob - zmaxp + 3).*(psi');
zprob  = zprob./sumc(zprob');
/*
************************************************************
*  Generate the z's
*  z[i] is MN(1,p_i)
*  p_ik propto det(lambda_k)^{-0.5}
*     *exp(-0.5(beta_i-theta_k)'lambda_k^{-1}(beta_i-theta_k))
*     *psi_k
*  Only accept a new  psi,  if they generate
*  a ordering of psi_1 <= psi_2 <= ....
************************************************************
*/
if mmod > 1;
	z      = rndzmn(zprob);
	classn = sumc( z .== mint' );		@ Number of subjects assigned to the classes @
	/*
	************************************************************
	*  DIRORD is ordered dirichlet. See /gauss/src/plrand.src
	*  {psi, xgam} = dirord(alpha, xgam);
	*   	alpha	= k x 1 vector of parameters.
	*	xgam	= n x k matrix of ordered gamma random deviates.
	*   	psi		= n x k matrix of ordered dirichlet probs.
	************************************************************
	*/
	{psi, xgam } = dirord(w0+classn,xgam);
	psi = psi';

	@ Enforce order restrictions on classn.  @
	@ Allow run of pcount violations before reorder @
	rclass = rankindx(classn,1);
	if not rclass == mint;			@ counts violate order restrictions @
		pflag = 1;					@ set flag for next iteration            @
		ipcount = ipcount + 1;
		if ipcount >= pcount;		@ let labels switch  @
			pflag  = 0;
			ipcount = 0;
			if nmcmc <= nblow;
				@ permute the assignments @
				classn = sortc(classn,1);
		       		z      = recode(z, z.== mint', rclass);
		       		psi[rclass] = psi;  @ reorder psi  @
					xgam[1, rclass] = xgam; 
		       		theta[.,rclass] = theta;  
		       		j = 1;
		       		lamb0  = lambda;
		       		lamb2 = lambdai;
		      		do while j <= mmod;
		    			lambda[iptgp[rclass[j],1]:iptgp[rclass[j],2],.] =
		              		lamb0[iptgp[j,1]:iptgp[j,2],.];  
		       			lambdai[iptgp[rclass[j],1]:iptgp[rclass[j],2],.] =
		              		lamb2[iptgp[j,1]:iptgp[j,2],.] ;  
		       			j = j + 1;
					endo;
			endif;			@ if igibbs <= nblow  @
		endif;				@ if ipcount >= pcount @
	else;				@ rclass .== mint -> reset counter & flag @
		ipcount = 0;
		pflag   = 0;
	endif;				@ if not rclass == mint @

endif;   					@ if mmod > 1 @

endp;

/*
*****************************************************************************************************
*  OUTPUTANAL
*	Does analysis of output and save some results 
****************************************************************************************************
*/
PROC (0) = outputanal;
local  bout, sout, ebeta, sbeta, cb, rmse, fmtn1, fmtn2, fmts1, fmts2, 
 knames,knames2,mmod2,a, a1,b1,iptgpt,
lambdatk,lambdamk,lambdask,fk,k,   clrate, hbclass, hbclassk,
betat, sigmat, classt, psit, thetat, lambdat, mmodt;


if flagtrue == 1;		@ Know true paramaters from simulation @
	load betat		= betat;
	load sigmat 	= sigmat;
	load classt		= classt;
	load psit		= psit;
	load thetat		= thetat;
	load lambdat	= lambdat;
	mmodt			= maxc(classt);		@ True number of mixture components @
endif;



knames = 0 $+ "Group " $+ ftocv(mint,1,0);
if flagtrue == 1;
	mmodt = maxc(classt);
	knames2	= 0 $+ "Group " $+ ftocv(seqa(1,1,mmodt),1,0);
endif;
@ Define formats for fancy printing @
@ Used to print a matrix of alpha & numeric variables @
let fmtn1[1,3]	= "*.*lf" 10 5;			@ Format for printing numeric variable 		@
let fmtn2[1,3]	= "*.*lf" 10  0;		@ Format for numeric variable, no decimal	@

let fmts1[1,3]	= "-*.*s" 10 9;			@ Format for alpha, left justify		 	@
let fmts2[1,3]  = "*.*s" 10 9;			@ Format for alpha, right justify		  	@
format 10,5;							@ Default pring format						@

output file = ^outfile reset;			@ outfile is the file handle for the output file 		@
								@ Route printed output to the defined by outfile 	@
print "Results from MIXREG.GSS";
print "Hierarchical Bayesian linear regression model using MCMC.";
print "Use mixture distribution for heterogeneity";
print "Different design matrices for each subject";
print "Y_i = X_i*beta_i + epsilon_i";
print "beta_i = sum_k psi_k N(beta_i | theta_k, lambda_k)";
print "epsilon_i is N(0,sigma2*I)";
print;
print "Number of components is fixed at: " mmod;
print;
print "Ouput file: " getpath(outfile);		@ File assigned to file handle outfile @
datestr(date);								@ Print the current data				@
print;
print;
print "-----------------------------------------------------------------------------------";
print;
print "MCMC Analysis";
print;
print "Total number of MCMC iterations           	   = " nmcmc;
print "Number of iterations used in the analysis 	   = " smcmc;
print "Number in transition period               	   = " nblow;
print "Number of iterations between saved iterations  = " skip-1;
print;
print "Number of subjects	                = " nsub;
print "Mean # of observations per subject  = " meanc(yrows);
print "STD  # of observations per subject  = " stdc(yrows);
print "MIN  # of observations per subject  = " minc(yrows);
print "MAX  # of observations per subject  = " maxc(yrows);
print "Total number of observations        = " ntot;
print "Number of independent variables X   = " xdim " (excluding intercept)";
print;
print "Dependent variable is " $ynames;
print;
print "Independent variables in first level equation: Y_i = X_i*beta_i + epsilon_i";
call sumstats(xnames,xdata,fmts1,fmts2,fmtn1);			@ Print summary statisticcs @
print;
print "Loglikelihood form MCMC:";
print "Number of mixture components = " mmod;
print "Posterior mean               = " loglikem;
print "Posterior STD                = " loglikes;
print; 
print "-----------------------------------------------------------------------------------";
print;
print "Statistics of Fit Measures for each Subject";
print "Average Predictive Correlation (Muptiple R) = " meanc(br);
print "STD of Predictive Correlations              = " stdc(br);
print "Average R-Squared                           = " meanc(br^2);
print "STD of R-Squared                            = " stdc(br^2);
print "Average Error Standard Deviation            = " meanc(estd);
print "STD of Error Standard Deviation             = " stdc(estd);
print;
print "-----------------------------------------------------------------------------------";
print "Estimation of the error STD sigma";
if flagtrue == 1;
	print "True Sigma     = " sigmat;
endif; 
print "MLE            = " sighat;
print "Posterior Mean = " sigmam;
print "Posterior STD  = " sigmas;
print;
print "-----------------------------------------------------------------------------------";
print "Statistics for Individual-Level Regression Coefficients";

if flagtrue == 1;
	ebeta	= meanc(betat');
	sbeta	= stdc(betat');
	print "True Beta";
	a	= {"Variable" "Mean" "STD" };
	call outitle(a,fmts1,fmts2);

	bout 	= xnames~ebeta~sbeta;
	call outmat(bout,fmts1,fmtn1);

endif;

print "MLE Coefficients ";
a	= {"Variable" "MeanMLE" "STDMLE"};
call outitle(a,fmts1,fmts2);
bout = xnames~meanc(bhat')~stdc(bhat');
call outmat(bout,fmts1,fmtn1);

print "HB Coefficients ";
a 	= {"Variable"  "PostMean"  "PostSTD"};
call outitle(a,fmts1,fmts2);
ebeta	= meanc(betam');
sbeta	= sqrt( meanc( (betas^2)') + stdc( betam')^2);

bout 	= xnames~ebeta~sbeta;
call outmat(bout,fmts1,fmtn1);

if flagtrue == 1;
	print "Comparison of True Beta to Individual Level Estimates";
	for i0 (1,rankx,1); i = i0;
		print "Variable is " $ xnames[i];
		cb 		= corrx( betat[i,.]'~betam[i,.]' );
		rmse 	= betat[i,.] - betam[i,.];
		rmse	= rmse*rmse';
		rmse	= sqrt(rmse/nsub);
		print "Correlation between true and HB  = " cb[1,2];
		print "RMSE between true and HB         = " rmse;
		print;
		cb 		= corrx( betat[i,.]'~bhat[i,.]' );
		rmse 	= betat[i,.] - bhat[i,.];
		rmse	= rmse*rmse';
		rmse	= sqrt(rmse/nsub);
		print "Correlation between true and MLE = " cb[1,2];
		print "RMSE between true and MLE        = " rmse;
		print;
	endfor;
endif;
print "-----------------------------------------------------------------------------------";
if mmod > 1;
	print "Estimated Group Probabilities psi";

	if flagtrue == 1;

		a		= "  "~(knames2');
		call outitle(a,fmts1,fmts2);
		bout	= "True"~(psit');
		call outmat(bout,fmts1,fmtn1);
	endif;
	a		= "   "~(knames');
	call outitle(a,fmts1,fmts2);

	bout	= "HB Mean"~(psim');
	bout	= bout|("HB STD"~(psis'));
	call outmat(bout,fmts1,fmtn1);
endif;

if flagtrue == 1;			@ Do Misclassification @
print "-----------------------------------------------------------------------------------";

	clrate = zeros(mmod,mmodt);
	hbclass	= maxindc(zprobm');			@ Get index that correponds to the maximum @
	for fk (1,mmodt,1); k =fk;
		hbclassk = selif(hbclass, classt .== k);
		if rows(hbclassk) > 0;
			clrate[.,k]	= sumc( hbclassk .== mint');
		endif;
	endfor;

	print "Classification Rates:  True versus Maximum HB Posterior Probability";
	a = 0 $+ "HB Groups ";
	a = a|(0 $+ "True " $+ ftocv(seqa(1,1,mmodt),1,0));
	a = a|"Total";
	a = a';
	call outitle(a,fmts1,fmts2);

	a		= knames|"Total";
	bout	= clrate~(sumc(clrate'));
	bout	= bout|(sumc(bout)');
	bout	= a~bout;
	call outmat(bout,fmts1,fmtn2);
	print;
endif;
print "-----------------------------------------------------------------------------------";
print "HB Estimates of Theta";

if flagtrue == 1;
	print "True Theta";
	a	= "Variable"~(knames2');
	call outitle(a,fmts1,fmts2);
	bout = xnames~thetat;
	call outmat(bout,fmts1,fmtn1);
endif;
print "Posterior Mean of Theta";
a 		= "Variable"~(knames');
call outitle(a,fmts1,fmts2);
bout = xnames~thetam;
call outmat(bout,fmts1,fmtn1);

print "Posterior STD of Theta";
call outitle(a,fmts1,fmts2);
bout = xnames~thetas;
call outmat(bout,fmts1,fmtn1);
print "-----------------------------------------------------------------------------------";
a		= "Variable"~(xnames');

print "HB Estimate of Lambda";

if flagtrue == 1;
	b1	= seqa(rankx,rankx,mmodt);
	a1	= 1|(1+b1[1:mmodt-1]);
	iptgpt	= a1~b1;
	for fk (1,mmodt,1); k = fk;
		lambdatk	= lambdat[iptgpt[k,1]:iptgpt[k,2],.];
		print "True Lambda for group " k;
		call outitle(a,fmts1,fmts2);
		bout	= xnames~lambdatk;
		call outmat(bout,fmts1,fmtn1);
	endfor;
endif;
print "-----------------------------------------------------------------------------------";
for fk (1,mmod,1); k = fk;
print "Posterior Mean of Lambda for group " k;
	lambdamk	= lambdam[iptgp[k,1]:iptgp[k,2],.];
	call outitle(a,fmts1,fmts2);
	bout	= xnames~lambdamk;
	call outmat(bout,fmts1,fmtn1);
print "Posterior STD of Lambda for group " k;
	lambdask	= lambdas[iptgp[k,1]:iptgp[k,2],.];
	call outitle(a,fmts1,fmts2);
	bout	= xnames~lambdask;
	call outmat(bout,fmts1,fmtn1);
print "-----------------------------------------------------------------------------------";
endfor;
print "===================================================================================";

output off;
closeall;
endp;

/*
*****************************************************************************************
* OUTITLE
*	Prints header for columns of numbers.
*	INPUT
*		a 	= character row vector of column names
*		fmts1	= format for first column
*		fmts2	= format for second column
*	OUTPUT
*		None
******************************************************************************************
*/
PROC (0) = outitle(a,fmt1,fmt2);
local mask, fmt, flag, ncols;
ncols	= cols(a);
mask	= zeros(1,ncols);
fmt		= fmt1|(ones(ncols-1,1).*.fmt2);
flag	= printfm(a,mask,fmt);
print;
endp;
/*
***************************************************************************************
* OUTMAT
*	Outputs a matrix:
*	(Character Vector)~(Numeric matrix);
*	The entries in the numeric matrix have the same format
*	INPUT
*		bout		= matrix to be printed
*		fmts		= format for string
*		fmtn		= format for numeric matrix
*	OUTPUT
*		None
******************************************************************************************
*/
PROC (0) = outmat(bout,fmts,fmtn);
local fmt,mask,flag,ncols, nrows;
ncols		= cols(bout);
nrows		= rows(bout);
fmt			= fmts|(ones(ncols-1,1).*.fmtn);
mask		= zeros(nrows,1)~ones(nrows,ncols-1);
flag		= printfm(bout,mask,fmt);
print;
endp;

/*
*****************************************************************************************
* SUMSTATS
*	Prints summary statistics for a data matrix
*	INPUT
*		names		= charater vector of names
*		data		= data matrix
*		fmts1		= format for string
*		fmts2		= format for string
*		fmtn		= format for numbers
*	OUTPUT
*		None
********************************************************************************************
*/
PROC (0) = sumstats(names,data,fmts1,fmts2,fmtn);
local a, bout;
a	= {"Variable" "Mean" "STD" "MIN" "MAX"};
call outitle(a,fmts1,fmts2);
bout	= names~meanc(data)~stdc(data)~minc(data)~maxc(data);
call outmat(bout,fmts1,fmtn);
endp;