MIXED EFFECTS MODELING

[JOADUTWUM]

I want to clarify the effect of growth rate on wood specific gravity. 

Three trees were partitioned into two halves each with the central core as the reference (this is affected by the growth rate and consequently named it plot). Wood samples are collected at different height levels at 1m interval starting from 0.3m from ground level. Wood specific gravity is measured at relative radial interval (RRD) (0.1 to 1) standardizing the different radii and aligning them to the same age when wood was formed. 

I assume PLOT to be fixed factor, RDD as a covariate and Repeated variable, Height as Random factor and repeated subject Height(Plot).

I wrote this procedure to execute the plan;

proc mixed data=teaksggrowthrate;
class plot rrd height;
model SG = plot rrd plot*rrd;
repeated rrd / subject= height (plot) r rcorr;
run; 

Am I doing the right thing? 

Thank you in advance.

[SteveDenham · Accepted solution]

If you wish to treat height as a random effect (and estimate the variance component) you may wish to fit a random intercept model for height.  I would adapt your code to try the following:

proc mixed data=teaksggrowthrate;
class plot rrd height;
model SG = plot rrd plot*rrd;
random intercept/subject=height;
repeated rrd / subject= height (plot) r rcorr;
lsmeans plot rrd plot*rrd;
run; 

I also added an LSMEANS statement to get the marginal means (averaged over all instances of height).

I am a bit curious as to treating height as a random effect.  You may want to consider some sort of 3 dimensional spatial correlation for rrd and height, jointly.

SteveDenham

[JOADUTWUM]

Dear Steve,
I also did random intercept model for height.
With height growth, wood is formed by the same cambium but in successive years.
PS: how do I implement the 3 dimensional spatial correlation for rrd and height?
Thank you.

[SteveDenham]

Well, a 3 dimensional spatial analysis isn't really needed, as you have radial symmetry within each plot.  Using the information found here , I would try something like this:

proc mixed data=teaksggrowthrate;
class plot ;
model SG = plot / solution;
repeated  / subject= intercept type=sp(powa)(height rrd) /r rcorr;
lsmeans plot ;
run; 

If you want LSMEANS at various points on the rrd axis, you have to employ a bit of trickery.  First, create a variable in your dataset that is exactly equal to rrd.  For this example, I will call it rrd_continuous.  Then fit the following:

proc mixed data=teaksggrowthrate;
class plot  rrd;
model SG = plot rrd plot*rrd/ solution;
repeated rrd / subject= intercept type=sp(powa)(height rrd_continuous) /r rcorr;
lsmeans plot ;
run; 

The covariance type in both models is an anisotropic (not the same correlation in both dimensions) power relationship.  Other covariance types that might be appropriate would be anisotropic exponential [ SP(POWA)(height rrd_continuous) ] and 2D exponential, geometrically anisotropic [ SP(EXPGA)(height rrd_continuous).  You could select amongst these by selecting the model with the smallest AIC or corrected AIC. Short descriptions of these structures are found in the documentation under the REPEATED statement at the TYPE= section.

Also, be aware that the arguments to the spatial covariance structures MUST be continuous, and so should not appear in the CLASS statement, which is why you have to have rrd_continuous in the dataset for the second MIXED example.

SteveDenham

[JOADUTWUM]

Dear Steve, 

Thank you for the support. 

With the first proc,  I got an error: 

ERROR 22-322: Syntax error, expecting one of the following: ;, GROUP, HLM, HLPS, LDATA, LOCAL, LOCALVAR, LOCALW, NONLOCALW, R, RC,

RCI, RCORR, RI, SSCP, SUBJECT, TYPE.

ERROR 200-322: The symbol is not recognized and will be ignored.

[SteveDenham]

That error is my fault.  There are two slashes in the REPEATED statements.  Remove the second one and that should take care of the issue.

SteveDenham

[JOADUTWUM]

PROC MIXED DATA=sgonage INFO IC ASYCORR ASYCOV ;
CLASS TREE RADII AGE ;
MODEL PBSG = RADII AGE RADII*AGE / S CHISQ ;
REPEATED / TYPE = un SUBJECT = TREE(RADII) R RCORR ;
RUN;

The statements above could not make Hessian positive definite.

Any help?

Thank you

[SteveDenham]

The most likely reason is that you do not have enough data to fit a fully unstructured covariance model.  For instance, 6 levels of AGE requires fitting 21 parameters, and at a rule of thumb of at least 10 obs per parameter, you would need around 210 observations.  You'll likely have to find a more parsimonious covariance structure.

If you change the REPEATED statement to something like:

REPEATED AGE/ TYPE = csh SUBJECT = TREE(RADII) R RCORR ;

to see if that avoided the problem, it would be a start.  Time series related structures- AR(1), ARH(1), ANTE(1)-are good choices for observations equally spaced in time.  For unequally spaced in time observations, a spatial power structure has worked for us in the past.  You can select among the various choices by selecting the one with the smallest AIC or AICC among those that converge and don't have the Hessian issue.

SteveDenham

[JOADUTWUM]

Good day Steve Denham ,

I trust this message finds you well. 

I have a question regarding doubly repeated measures.

I used UN@AR(1) and I got a negative covariance in the UN structure. 

Is it okay ?

How would I explain that? 

proc mixed data=shproportion method=reml covtest cl info ic plots=vcirypanel ;
class tree aspect height ;
model HW_ = height aspect height*aspect / s chisq cl influence residual;
random tree / g gcorr v vcorr ;
repeated aspect height / type = UN@AR(1) subject=tree r rcorr ;
lsmeans aspect height / pdiff=all ;
contrast 'aspect 1 vs 2 intercept' aspect 1 -1 ;
contrast 'aspect 1 vs 2 slope' aspect 1 -1 ;
run;

2,1 is negative. 

Thank you.

[jiltao]

A negative value of UN(2,1) indicates that the covariance is negative. A negative covariance is okay.