topic Re: Bayesian implementation of UN@AR(1) model in Statistical Procedures

Bayesian implementation of UN@AR(1) model

JOADUTWUM — Thu, 10 Nov 2022 13:07:38 GMT

Dear all:

I hope you all are doing well and keeping safe.

I would like to implement this mixed model in Bayesian framework. The codes are below:

***********************************************************************************************

/*Multivariate RM Models with Kronecker product covariance*/
proc mixed data=sg2 covtest cl;
class tree aspect height time ;
model y = height|time / noint ;

random tree ;
repeated height time / type = un@ar(1) subject=tree(aspect) r rcorr;
run;

********************************************************************************************

Height factor has 6 levels; time factor has 10 levels and aspect has 2 levels.

Thank you.

Re: Bayesian implementation of UN@AR(1) model

SteveDenham — Mon, 21 Nov 2022 18:51:37 GMT

The only thing I can think of is to look at PROC BGLIMM, and model the unstructured part of the Kronecker expression with a RANDOM statement and the AR(1) part with a REPEATED statement.

proc bglimm data=sg2;
class tree aspect height time ;
model y = height|time / noint ;

random intercept height/type=un subject=tree ;
repeated  time / type =ar(1) subject=tree(aspect*height) r rcorr ;
run;

Absolutely no guarantees on this, but if you do get it to give reasonable results, I would really, REALLY like to know of any tweaks you might have used.

SteveDenham

Re: Bayesian implementation of UN@AR(1) model

JOADUTWUM — Tue, 22 Nov 2022 22:14:54 GMT

Dear @SteveDenham

Thank you for the helpful reply.

I adjusted the statements as:

/*PROC BGLIMM default*/
proc bglimm data=rrdsg seed=475193 nmc=500000 nbi=5000 thin=100 nthreads = 4 dic diag = (mcse ess raftery geweke)
outpost = pr_out plots = (trace autocorr density) STATISTICS(PERCENT=(2.5 50 97.5));
class tree aspect height time;
model pbsg = height|rrd / noint cprior=normal(var=1e4);
random intercept / subject = tree nuts;
random intercept height / type=un subject=tree(aspect) gcorr nuts;
repeated time / type=ar(1) subject=tree(aspect*height) rcorr;
run;

One problem is the diagnostics for the fixed height effects. Dependence is too high.

Re: Bayesian implementation of UN@AR(1) model

SteveDenham — Mon, 28 Nov 2022 17:59:27 GMT

The only thing in the code I see that might be causing that are these two statements:

random intercept / subject = tree nuts;
random intercept height / type=un subject=tree(aspect) gcorr nuts;

What happens if you shift those around a bit to:

random intercept / subject = tree nuts;
random height / type=un subject=tree(aspect) gcorr nuts;

This ought to treat height more like a G side repeated measure, and remove the two level intercepts for tree and aspect within tree, as I suspect you may be short of enough data to fit both of those and height.

SteveDenham