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 ;
... View more
