proc mcmc: An exception has been encountered - Read Access Violation - Write Access Violation

[aaraujo]

I am trying to use "proc mcmc" to specify a correlated random intercept and slope model, similar to the one exemplified on this page:

http://support.sas.com/documentation/cdl/en/statug/68162/HTML/default/viewer.htm#statug_mcmc_example...

The purpose is to specify more complex covariance structures of the residual error. I have written a program. The problem does not run. Immediately after submitting I get Read Access Violation and Write Access Violation errors. So, what is wrong?

Here is the program with the dataset:

data WORK.ndata;
	infile datalines;
	input Patient Treatment $ Cycle $ Period $ Y;
	datalines;
		1 A 1 1 2394
		1 B 1 2 2686
		1 A 2 3 2515
		1 B 2 4 2675
		1 A 3 5 2583
		1 B 3 6 2802
		2 A 1 1 2746
		2 B 1 2 2726
		2 A 2 3 2592
		2 B 2 4 2867
		2 A 3 5 2743
		2 B 3 6 2742
		3 A 1 1 2668
		3 B 1 2 2560
		3 A 2 3 2542
		3 B 2 4 2584
		3 A 3 5 2491
		3 B 3 6 2737
		4 A 1 1 2397
		4 B 1 2 2696
		4 A 2 3 2411
		4 B 2 4 2895
		4 A 3 5 2499
		4 B 3 6 2760
		5 A 1 1 3179
		5 B 1 2 3221
		5 A 2 3 2952
		5 B 2 4 3096
		5 A 3 5 2600
		5 B 3 6 3192
		6 A 1 1 2643
		6 B 1 2 2496
		6 A 2 3 2759
		6 B 2 4 2847
		6 A 3 5 2651
		6 B 3 6 2860
		7 A 1 1 2678
		7 B 1 2 2843
		7 A 2 3 2492
		7 B 2 4 2763
		7 A 3 5 2801
		7 B 3 6 2890
		8 A 1 1 2887
		8 B 1 2 2862
		8 A 2 3 2875
		8 B 2 4 3083
		8 A 3 5 2689
		8 B 3 6 2967
		9 A 1 1 2490
		9 B 1 2 2841
		9 A 2 3 2648
		9 B 2 4 3044
		9 A 3 5 2688
		9 B 3 6 2914
		10 A 1 1 2268
		10 B 1 2 2576
		10 A 2 3 2413
		10 B 2 4 2493
		10 A 3 5 2344
		10 B 3 6 2699
		11 A 1 1 2617
		11 B 1 2 2923
		11 A 2 3 2629
		11 B 2 4 2832
		11 A 3 5 2732
		11 B 3 6 2866
		12 A 1 1 2627
		12 B 1 2 2759
		12 A 2 3 2712
		12 B 2 4 2698
		12 A 3 5 2572
		12 B 3 6 2826
run;

proc sort data=WORK.ndata;
	by Patient Period;
run;

* Determine number of patients and number of observations ;
data _NULL_;
	set WORK.ndata nobs=_n;
	by Patient;
	retain num 0;
	if first.Patient then num+1;
	call symput('PatientNumber', num);
	call symputx('nobs', _n);
run;

* Create random effects design matrix ;
data WORK.design_matrix;
	set WORK.ndata;
	array b_[&PatientNumber];
	array d_[&PatientNumber];
	do i = 1 to &PatientNumber by 1;
		if Patient = i then do;
			b_[i] = 1;
			if Treatment = "B" then d_[i] = 1;
			else d_[i] = 0;
		end;
		else do;
			b_[i] = 0;
			d_[i] = 0;
		end;
	end;
	keep b_: d_:;
run;

* Create dataset with first and last index ;
data WORK.index;
	set WORK.ndata;
	by Patient;
	retain num;
	if first.Patient then num = _n_;
	if last.Patient then do;
		first = num;
		last = _n_;
		output;
	end;
	keep first last;
run;

* Create dummy dataset ;
data dummy;
run;

proc mcmc data=dummy outpost=outAR1_jointllike seed=3141593 nmc=10000 thin=1 nbi=0 jointllike;
	array ALPHA[%eval(2*&PatientNumber)] b_1-b_%eval(&PatientNumber) d_1-d_%eval(&PatientNumber);
	array MU[&nobs];
	array COV[&nobs,&nobs];
	
	array M[%eval(2*&PatientNumber)];
	array G[%eval(2*&PatientNumber),%eval(2*&PatientNumber)];
	
	begincnst;
		call zeromatrix(G);
		do i = 1 to &PatientNumber by 1;
			M[i] = gamma_00;
			M[i+&PatientNumber] = gamma_10;
			G[i,i] = sigma_u0;
			G[i+&PatientNumber,i+&PatientNumber] = sigma_u1;
			G[i,i+&PatientNumber] = sigma_u0u1;
			G[i+&PatientNumber,i] = sigma_u1u0;
		end;
	endcnst;
	
	array TAU[2] gamma_00 gamma_10;
	array SIGMA[2,2] sigma_u0 sigma_u0u1 sigma_u1u0 sigma_u1;
	
	array THETA[2];
	array DELTA[2,2];
	array PSI[2,2];
	
	array data[&nobs] / nosymbols;
	array F[&PatientNumber] / nosymbols;
	array L[&PatientNumber] / nosymbols;
	array XZ[&nobs,%eval(2*&PatientNumber)] / nosymbols;
	
	begincnst;
		call zeromatrix(THETA);
		call identity(DELTA);
		call mult(DELTA, 1e7, DELTA);
		call identity(PSI);
		call mult(PSI, 0.01, PSI);
		rc = read_array("WORK.ndata", data, "Y");
		rc = read_array("WORK.index", F, "first");
		rc = read_array("WORK.index", L, "last");
		rc = read_array("WORK.design_matrix", XZ);
	endcnst;
	
	parms TAU;
	parms SIGMA;
	parms ALPHA;
	parms sigma_e;
	parms rho 0.5;
	
	hyperprior TAU ~ mvn(THETA, DELTA);
	hyperprior SIGMA ~ iwish(2, PSI);
	
	prior ALPHA ~ mvn(M, G);
	prior sigma_e ~ igamma(0.01, scale=0.01);
	prior rho ~ uniform(left=-1, right=1);
	
	call zeromatrix(COV);
	do k = 1 to &PatientNumber by 1;
		do i = F[k] to L[k] by 1;
			COV[i,i] = sigma_e;
			do j = i+1 to L[k] by 1;
				COV[i,j] = sigma_e*rho**abs(i-j);
				COV[j,i] = COV[i,j];
			end;
		end;
	end;
	
	call mult(XZ, ALPHA, MU);
	
	llike = lpdfmvn(data, MU, COV);
	
	model general(llike);
run;

%symdel PatientNumber nobs;

proc sql;
	drop table WORK.design_matrix, WORK.index, WORK.dummy;
quit;

[JacobSimonsen]

The error comes because of your way of reading in data. You dont give the procedure your data in the usual way (data=..). Instead, you attemp to read in data with the read_array function. According to the documentation, the read_array function can only be used inside a function or call routine, which in means that you only can use read_array inside proc fcmp.

This is from the documentation:

"The READ_ARRAY function attempts to dynamically resize the array to match the dimensions of the input data set. This means that the array must be dynamic. That is, the array must be declared either in a function or CALL routine or declared with the /NOSYMBOLS option. "

[aaraujo]

Thanks a lot for helping out.

I am not sure that is the problem. The arrays are declared with the /nosymbols option. The problem must lie in the way the random effects are declared. I have written an alternative program where a distinct model is specified. Here the random effects are defined distinctively. The program uses the "read_array" function inside "proc mcmc". The difference is that this program runs perfectly fine. Check it out.

data WORK.ndata;
	infile datalines;
	input Patient Treatment $ Cycle $ Period $ Y;
	datalines;
		1 A 1 1 2394
		1 B 1 2 2686
		1 A 2 3 2515
		1 B 2 4 2675
		1 A 3 5 2583
		1 B 3 6 2802
		2 A 1 1 2746
		2 B 1 2 2726
		2 A 2 3 2592
		2 B 2 4 2867
		2 A 3 5 2743
		2 B 3 6 2742
		3 A 1 1 2668
		3 B 1 2 2560
		3 A 2 3 2542
		3 B 2 4 2584
		3 A 3 5 2491
		3 B 3 6 2737
		4 A 1 1 2397
		4 B 1 2 2696
		4 A 2 3 2411
		4 B 2 4 2895
		4 A 3 5 2499
		4 B 3 6 2760
		5 A 1 1 3179
		5 B 1 2 3221
		5 A 2 3 2952
		5 B 2 4 3096
		5 A 3 5 2600
		5 B 3 6 3192
		6 A 1 1 2643
		6 B 1 2 2496
		6 A 2 3 2759
		6 B 2 4 2847
		6 A 3 5 2651
		6 B 3 6 2860
		7 A 1 1 2678
		7 B 1 2 2843
		7 A 2 3 2492
		7 B 2 4 2763
		7 A 3 5 2801
		7 B 3 6 2890
		8 A 1 1 2887
		8 B 1 2 2862
		8 A 2 3 2875
		8 B 2 4 3083
		8 A 3 5 2689
		8 B 3 6 2967
		9 A 1 1 2490
		9 B 1 2 2841
		9 A 2 3 2648
		9 B 2 4 3044
		9 A 3 5 2688
		9 B 3 6 2914
		10 A 1 1 2268
		10 B 1 2 2576
		10 A 2 3 2413
		10 B 2 4 2493
		10 A 3 5 2344
		10 B 3 6 2699
		11 A 1 1 2617
		11 B 1 2 2923
		11 A 2 3 2629
		11 B 2 4 2832
		11 A 3 5 2732
		11 B 3 6 2866
		12 A 1 1 2627
		12 B 1 2 2759
		12 A 2 3 2712
		12 B 2 4 2698
		12 A 3 5 2572
		12 B 3 6 2826
run;

proc sort data=WORK.ndata;
	by Patient Period;
run;

* Determine number of patients and number of observations ;
data _NULL_;
	set WORK.ndata nobs=_n;
	by Patient;
	retain num 0;
	if first.Patient then num+1;
	call symput('PatientNumber', num);
	call symputx('nobs', _n);
run;

* Create fixed effects random effects design matrix ;
data WORK.design_matrix;
	set WORK.ndata;
	B0 = 1;
	if Treatment = "B" then B1 = 1;
	else B1 = 0;
	array b_[&PatientNumber];
	array d_A_[&PatientNumber];
	array d_B_[&PatientNumber];
	do i = 1 to &PatientNumber by 1;
		if Patient = i then do;
			b_[i] = 1;
			if Treatment = "B" then do;
				d_B_[i] = 1;
				d_A_[i] = 0;
			end;
			else do;
				d_B_[i] = 0;
				d_A_[i] = 1;
			end;
		end;
		else do;
			b_[i] = 0;
			d_B_[i] = 0;
			d_A_[i] = 0;
		end;
	end;
	keep B0 B1 b_: d_:;
run;

* Create dataset with first and last index ;
data WORK.index;
	set WORK.ndata;
	by Patient;
	retain num;
	if first.Patient then num = _n_;
	if last.Patient then do;
		first = num;
		last = _n_;
		output;
	end;
	keep first last;
run;

* Create dummy dataset ;
data dummy;
run;

proc mcmc data=dummy outpost=outAR1_jointllike seed=3141593 nmc=10000 thin=1 nbi=0
	jointllike monitor=(B0 B1 sb sd sigma rho ite);
	array ALPHA[%eval(3*&PatientNumber+2)] B0 B1 b_1-b_%eval(&PatientNumber)
		d_A_1-d_A_%eval(&PatientNumber) d_B_1-d_B_%eval(&PatientNumber);
	array MU[&nobs];
	array COV[&nobs,&nobs];
	
	array ite[&PatientNumber];
	
	array data[&nobs] / nosymbols;
	array F[&PatientNumber] / nosymbols;
	array L[&PatientNumber] / nosymbols;
	array XZ[&nobs,%eval(3*&PatientNumber+2)] / nosymbols;
	
	begincnst;
		rc = read_array("WORK.ndata", data, "Y");
		rc = read_array("WORK.index", F, "first");
		rc = read_array("WORK.index", L, "last");
		rc = read_array("WORK.design_matrix", XZ);
	endcnst;
	
	parms B0 B1;
	parms sb 1000;
	parms sd 1000;
	parms sigma 1000 / slice;
	parms rho 0.5 / slice;
	parms b_1-b_%eval(&PatientNumber);
	parms d_A_1-d_A_%eval(&PatientNumber) d_B_1-d_B_%eval(&PatientNumber);
	
	hyperprior sb ~ general(0, lower=0);
	hyperprior sd ~ general(0, lower=0);
	
	prior B0 B1 ~ normal(0, var=1e7);
	prior b_: ~ normal(0, sd=sb);
	prior d_: ~ normal(0, sd=sd);
	prior sigma ~ general(0, lower=0);
	prior rho ~ uniform(left=-1, right=1);
	
	call zeromatrix(COV);
	call mult(XZ, ALPHA, MU);
	
	do k = 1 to &PatientNumber by 1;
		ite[k] = ALPHA[2] + ALPHA[2+2*&PatientNumber+k] - ALPHA[2+&PatientNumber+k];
		do i = F[k] to L[k] by 1;
			COV[i,i] = sigma**2;
			do j = i+1 to L[k] by 1;
				COV[i,j] = sigma**2*rho**abs(i-j);
				COV[j,i] = COV[i,j];
			end;
		end;
	end;
	
	llike = lpdfmvn(data, MU, COV);
	
	model general(llike);
run;

%symdel PatientNumber nobs;

proc sql;
	drop table WORK.design_matrix, WORK.index, WORK.dummy;
quit;

[aaraujo]

This is another attempt which does not work. This is based on a Stan program. This results in read access and write access violations.

* Options for macro debugging ;
options symbolgen mprint mlogic;

data WORK.ndata;
	infile datalines;
	input Patient Treatment $ Cycle Period Y;
	datalines;
		1 A 1 1 2394
		1 B 1 2 2686
		1 A 2 3 2515
		1 B 2 4 2675
		1 A 3 5 2583
		1 B 3 6 2802
		2 A 1 1 2746
		2 B 1 2 2726
		2 A 2 3 2592
		2 B 2 4 2867
		2 A 3 5 2743
		2 B 3 6 2742
		3 A 1 1 2668
		3 B 1 2 2560
		3 A 2 3 2542
		3 B 2 4 2584
		3 A 3 5 2491
		3 B 3 6 2737
		4 A 1 1 2397
		4 B 1 2 2696
		4 A 2 3 2411
		4 B 2 4 2895
		4 A 3 5 2499
		4 B 3 6 2760
		5 A 1 1 3179
		5 B 1 2 3221
		5 A 2 3 2952
		5 B 2 4 3096
		5 A 3 5 2600
		5 B 3 6 3192
		6 A 1 1 2643
		6 B 1 2 2496
		6 A 2 3 2759
		6 B 2 4 2847
		6 A 3 5 2651
		6 B 3 6 2860
		7 A 1 1 2678
		7 B 1 2 2843
		7 A 2 3 2492
		7 B 2 4 2763
		7 A 3 5 2801
		7 B 3 6 2890
		8 A 1 1 2887
		8 B 1 2 2862
		8 A 2 3 2875
		8 B 2 4 3083
		8 A 3 5 2689
		8 B 3 6 2967
		9 A 1 1 2490
		9 B 1 2 2841
		9 A 2 3 2648
		9 B 2 4 3044
		9 A 3 5 2688
		9 B 3 6 2914
		10 A 1 1 2268
		10 B 1 2 2576
		10 A 2 3 2413
		10 B 2 4 2493
		10 A 3 5 2344
		10 B 3 6 2699
		11 A 1 1 2617
		11 B 1 2 2923
		11 A 2 3 2629
		11 B 2 4 2832
		11 A 3 5 2732
		11 B 3 6 2866
		12 A 1 1 2627
		12 B 1 2 2759
		12 A 2 3 2712
		12 B 2 4 2698
		12 A 3 5 2572
		12 B 3 6 2826
run;

data WORK.outdata;
	set WORK.ndata;
	if Treatment = "B" then TreatmentB = 1;
	else TreatmentB = 0;
	if Cycle = 2 then Cycle2 = 1;
	else Cycle2 = 0;
	if Cycle = 3 then Cycle3 = 1;
	else Cycle3 = 0;
	drop Treatment Cycle;
run;

proc sort data=WORK.outdata;
	by Patient Period;
run;

* Determine number of patients and number of observations ;
data _NULL_;
	set WORK.outdata nobs=_n;
	by Patient;
	retain num 0;
	if first.Patient then num+1;
	call symput('PatientNumber', num);
	call symputx('nobs', _n);
run;

* Create dataset with first and last index ;
data WORK.index;
	set WORK.outdata;
	by Patient;
	retain num;
	if first.Patient then num = _n_;
	if last.Patient then do;
		first = num;
		last = _n_;
		output;
	end;
	keep first last;
run;

* Create dummy dataset ;
data dummy;
run;

%macro mcmc_array;
	%do i = 1 %to &PatientNumber %by 1;
		array ALPHA&i._[2];
	%end;
%mend;

%macro mcmc_parms;
	parms
	%do i = 1 %to &PatientNumber %by 1;
		ALPHA&i._
	%end;
	;
%mend;

%macro mcmc_prior;
	prior
	%do i = 1 %to &PatientNumber %by 1;
		ALPHA&i._
	%end;
	~ mvn(MU, PSI);
%mend;

%macro mcmc_average;
	%do k = 1 %to &PatientNumber %by 1;
		do i = F[&k] to L[&k] by 1;
			average[i] = ALPHA&k._[1]+ALPHA&k._[2]*treat[i];
			COV[i,i] = sigma;
			do j = i+1 to L[&k] by 1;
				COV[i,j] = sigma*rho**abs(i-j);
				COV[j,i] = COV[i,j];
			end;
		end;
	%end;
%mend;

proc mcmc data=dummy outpost=outAR1_jointllike seed=3141593 nmc=10000 thin=1 nbi=0 jointllike;
	%mcmc_array
	
	array average[&nobs];
	array COV[&nobs,&nobs];
	
	array MU[2];
	array PSI[2,2];
	
	array THETA[2];
	array ETA[2,2];
	array DELTA[2,2];
	
	array treat[&nobs] / nosymbols;
	array outcome[&nobs] / nosymbols;
	array F[&PatientNumber] / nosymbols;
	array L[&PatientNumber] / nosymbols;
	
	begincnst;
		call zeromatrix(THETA);
		call identity(ETA);
		call mult(ETA, 1e7, ETA);
		call identity(DELTA);
		call mult(DELTA, 0.01, DELTA);
		rc = read_array("WORK.outdata", treat, "TreatmentB");
		rc = read_array("WORK.outdata", outcome, "Y");
		rc = read_array("WORK.index", F, "first");
		rc = read_array("WORK.index", L, "last");
	endcnst;
	
	parms MU;
	%mcmc_parms
	parms PSI;
	parms sigma;
	parms rho 0.5;
	
	hyperprior MU ~ mvn(THETA, ETA);
	hyperprior PSI ~ iwish(2, DELTA);
	
	%mcmc_prior
	prior sigma ~ igamma(0.01, scale=0.01);
	prior rho ~ uniform(left=-1, right=1);
	
	call zeromatrix(COV);
	%mcmc_average
	
	llike = lpdfmvn(outcome, average, COV);
	
	model general(llike);
run;

%symdel PatientNumber nobs;

* Delete macros from WORK.Sasmacr ;
proc catalog cat=WORK.Sasmacr;
	delete mcmc_array mcmc_parms mcmc_prior mcmc_average / et=macro;
quit;

proc sql;
	drop table
		WORK.ndata,
		WORK.outdata,
		WORK.index,
		WORK.dummy,
		WORK.outAR1_jointllike;
quit;

%put _GLOBAL_;
%put _LOCAL_;

[aaraujo · Accepted solution]

I eventually was able to write a program to fit a correlated random intercept slope with AR1 within subject residual errors, using proc mcmc. I had to resort to the macro facility to declared arrays for each subject. The following program works.

title "Correlated random intercept and slope";
title2 "AR1 stochastic process";
title3 "Bayesian inference";

* define dataset names here ;
%let datain=WORK.gaba_placebo; * input dataset ;
%let datatemp=WORK.gaba_placebo_temp; * temporary dataset ;
%let datadummy=WORK.dummy; * dummy dataset ;
%let dataout=WORK.outCRISAR1; * output dataset ;

* define variable names here ;
%let Subject=id; * subject variable ;
%let Treatment=Treatment; * treatment variable ;
%let Period=Period; * period variable ;
%let Outcome=Pain; * outcome variable ;
%let ref=Placebo; * reference treatment ;
%let seed=3141593; * random seed for simulation ;

* Determine and count unique treatments ;
proc sql noprint;
	select distinct &Treatment
			into :Treatment1 - :Treatment999999999
	from &datain
	where &Outcome is not null;
	%let TreatmentNumber = &SqlObs; * number of treatments ;
quit;

%macro findRef;
	%if &Treatment1 ne &ref %then %do;
		%local t f;
		%let t = 1;
		%do %until (&t eq &TreatmentNumber or &&Treatment&t eq &ref);
			%let t=%eval(&t+1);
		%end;
		%let f=&Treatment1;
		%let Treatment1=&ref;
		%let Treatment&t=&f;
	%end;
%mend findRef;

* Set first treatment to reference treatment ;
%findRef

%macro dataTemp;
	%local t;
	proc sql;
		create table &datatemp as
		select 	&Subject,
				case
				%do t = 1 %to &TreatmentNumber %by 1;
					when &Treatment eq "&&Treatment&t" then &t
				%end;
				end as &Treatment,
				&Period,
				&Outcome
		from &datain
		where &Outcome is not null
		order by &Subject, &Period;
	quit;
%mend dataTemp;

* Convert treatment variable to numeric,
and remove missing observations;
%dataTemp

* Determine and count unique subjects,
and number of observations ;
proc sql noprint;
	select	distinct &Subject
			into :Subject1 - :Subject999999999
	from &datatemp;
	%let SubjectNumber = &SqlObs; * number of subjects ;
	select count(&Outcome) into :size1 - :size999999999
	from &datatemp
	group by &Subject; * number of observations per subject ;
	select count(*) into :nobs /* total number of observations */
	from &datatemp;
quit;

* Create dummy dataset ;
data &datadummy;
run;

%macro mcmcArray;
	%local s t;
	%do s = 1 %to &SubjectNumber %by 1;
		array Y_&&Subject&s..[&&size&s];
		array MU_&&Subject&s..[&&size&s];
		array COV_&&Subject&s..[&&size&s,&&size&s];
		array R_&&Subject&s..[&TreatmentNumber] bi_&&Subject&s
		%do t = 2 %to &TreatmentNumber %by 1;
			ite_&&Subject&s.._&&Treatment&t
		%end;
		;
	%end;
	array B[&TreatmentNumber] intercept
	%do t = 2 %to &TreatmentNumber %by 1;
		pte_&&Treatment&t
	%end;
	;
	array G[&TreatmentNumber,&TreatmentNumber];
	array M[&TreatmentNumber];
	array S[&TreatmentNumber,&TreatmentNumber];
	array U[&TreatmentNumber,&TreatmentNumber];
	array treat[&nobs] / nosymbols;
	array data[&nobs] / nosymbols;
	begincnst;
		rc = read_array("&datatemp", treat, "&Treatment");
		rc = read_array("&datatemp", data, "&Outcome");
	endcnst;
%mend mcmcArray;

%macro mcmcParms;
	%local s;
	parms B;
	parms G;
	parms var_residual rho;
	parms
	%do s = 1 %to &SubjectNumber %by 1;
		R_&&Subject&s
	%end;
	;
%mend mcmcParms;

%macro mcmcPrior;
	%local s;
	beginnodata;
		hyperprior B ~ mvn(M, S);
		hyperprior G ~ iwish(&TreatmentNumber, U);
		prior
		%do s = 1 %to &SubjectNumber %by 1;
			R_&&Subject&s
		%end;
		~ mvn(B, G);
		prior var_residual ~ igamma(0.01, scale=0.01);
		prior rho ~ uniform(left=-1, right=1);
	endnodata;
%mend mcmcPrior;

%macro mcmcLogLike;
	%local s t;
	index = 1;
	llike = 0;
	%do s = 1 %to &SubjectNumber %by 1;
		do i = 1 to &&size&s by 1;
			Y_&&Subject&s..[i] = data[index];
			MU_&&Subject&s..[i] = bi_&&Subject&s
			%do t = 2 %to &TreatmentNumber %by 1;
				+ ite_&&Subject&s.._&&Treatment&t*(treat[index] eq &t)
			%end;
			;
			COV_&&Subject&s..[i,i] = var_residual;
			do j = i+1 to &&size&s by 1;
				COV_&&Subject&s..[i,j] = var_residual*rho**abs(i-j);
				COV_&&Subject&s..[j,i] = COV_&&Subject&s..[i,j];
			end;
			index = index + 1;
		end;
		llike = llike + lpdfmvn(Y_&&Subject&s, MU_&&Subject&s, COV_&&Subject&s);
	%end;
	model general(llike);
%mend mcmcLogLike;

* Fit linear mixed-effects model ;
proc mcmc
		data=&datadummy
		outpost=&dataout
		nbi=1000 /* number of burn-in iterations */
		nmc=10000 /* number of mcmc iterations */
		ntu=1000 /* number of turning iterations */
		seed=&seed /* random seed for simulation */
		thin=1 /* thinning rate */
		jointmodel; /* specify joint log-likelihood */
	%mcmcArray
	
	begincnst;
		call zeromatrix(M);
		call identity(S);
		call mult(S, 1e7, S);
		call identity(U);
		call mult(U, 0.01, U);
	endcnst;
	
	%mcmcParms
	
	%mcmcPrior
	
	%mcmcLogLike
run;

* Delete datasets from system ;
proc sql;
	drop table &datatemp, &datadummy, &dataout;
quit;

%symdel datain datatemp datadummy dataout;
%symdel Subject Treatment Period Outcome ref seed nobs;

%macro deleteVariables;
	%local i;
	%do i = 1 %to &SubjectNumber %by 1;
		%symdel Subject&i size&i;
	%end;
	%do i = 1 %to &TreatmentNumber %by 1;
		%symdel Treatment&i;
	%end;
	%symdel SubjectNumber TreatmentNumber;
%mend deleteVariables;

* Delete macro variables from system ;
%deleteVariables

* Delete macros from WORK.Sasmacr ;
proc catalog cat=WORK.Sasmacr;
	delete findRef dataTemp mcmcArray
		mcmcParms mcmcPrior mcmcLogLike
		deleteVariables / et=macro;
quit;

title; * clear all titles ;