topic PROX MIXED & GLIMMIX Repeated measures, interactions, and linear/quad. contrasts in Statistical Procedures

PROX MIXED & GLIMMIX Repeated measures, interactions, and linear/quad. contrasts

AgReseach7 — Mon, 10 Oct 2016 22:43:29 GMT

Appologies in advance for the long discussion, but here it goes (some of the data can be uploaded if needed):

Study design (lamb feeding trial):

Effects of using 2 different feed ingredients (juniper and urea) in supplements fed to ewe lambs on the following dependent variables: intake, growth, and blood serum (e.g., glucose).

Animal is exp. unit (lambs fed in individual pens). Each lamb has own unique ID
Randomized design
All lambs fed hay. Each lamb also fed respective treatment: 1 of 8 different feeds in a 4×2 factorial: 4 juniper levels (15, 30, 45, or 60%) and 2 urea levels (1 or 3%).
- TRT: J15U1, J15U3, J30U1, J30U3, J45U1, J45U3, J60U1, J60U3
Feed Intake and Growth evaluated on d 0, 5, 12, 19, 26, 33, and 40.
Blood serum evaluated on d 6, 7, 8 (PERIOD 1) and 20, 21, and 22 (PERIOD 2). At each day, blood was collected & analyzed at hour 0 and 6. Thus, 12 total blood draws. Mostly interested in Hour within PERIOD effect for blood serum, not so much in differences of day within period. Thus, will leave DAY out of this blood serum model, but include PERIOD and HOUR.

Objectives:

Does increasing level of juniper in the supplement result in linear or quadratic trends in intake, growth, or blood serum constituents (e.g., glucose)?
Does the response (increasing juniper) change over days on trial?
Does the response (increasing juniper) change due to level of urea?
Additional objective for blood serum: does sampling hour make a difference in the serum?

Proposed plan for the Growth data (blood serum similar):

Run the full model: DAY|JUN|UREA and evaluate linear/quad. contrasts as far down the line as needed:
1. if DAYxJUNxUREA and JUNxUREA are significant: run linear/quad. contrasts for the 4 JUN levels at each DAY & for each level of UREA.
2. if DAYxJUNxUREA significant but JUNxUREA not sign.: run contrasts for the 4 juniper levels at each day (averaged across both levels of urea) and also a contrast for UREA1 vs. UREA3
3. if DAYxJUNxUREA not significant and JUNxUREA significant: run juniper contrasts (averaged across days) at each level of UREA.
4. if DAYxJUNxUREA and JUNxUREA not significant: run juniper contrasts (averaged across urea and days), and also a contrast for UREA1 vs. UREA3
For Growth Data, I’m running PROC MIXED (SAS Editor below). Data look good (UNIVARIATE) except for gain to feed ratio; using ASIN Square root.

Questions:

Do I stick with basic linear & Quadratic coefficient terms and just look at JUN trends according to what the model significance is (As discussed above) OR do I run the Full Model and code contrast statements according to above discussion. My current plan is to test whether or not DAYxJUNxUREA is significant and if not, drop it from the model.
Can I get help with contrast statements if it involves more than just JUN (I know how to run PROC IML)? I’ve never made contrast statements for interactions (much less 3 interactions). Studied it the past couple of days & only confused myself.
What about LSMESTIMATE? I’ve never used it before, but…
To reduce complexity, I may go back to original plan & analyze all possible comparisons; just drop the codes JUN and UREA and code the treatments as J15U1, J15U3, J30U1… etc. Use “Simulate” to protect the “familywise error rate.”
Code to back-transform ARSIN SQRT data in Model AA below?

/*AA: FULL MODEL*/

OPTIONS PAGESIZE=60 LINESIZE=90;

DATA LAMBGROW; SET GROW;

asinGF= ARSIN(SQRT((GFkgkg/1000)+.5));

PROC SORT; BY DAY ID JUN UREA;

RUN;QUIT;

PROC MIXED;

CLASS DAY ID JUN UREA;

MODEL Intake = DAY JUN UREA JUN*UREA JUN*DAY UREA*DAY JUN*UREA*DAY/DDFM=KR;

REPEATED DAY/SUBJECT=ID TYPE=un; /*“picked the structure with the best AIC & BIC”*/

CONTRAST 'LINEAR' JUN 0 -0.668116 -0.224944 0.218229 0.6748311;

CONTRAST 'QUADRADIC' JUN 0 0.5022684 -0.491797 -0.508131 0.4976596;

CONTRAST 'UREA1 VS. UREA1' UREA 1 -1 0 0 0; /*Need help coding these correctly*/

LSMEANS DAY JUN UREA JUN*UREA JUN*DAY UREA*DAY JUN*UREA*DAY/DIFF

ADJUST=SIMULATE (REPORT SEED=121211) cl adjdfe=row SLICE=(JUN*UREA JUN*DAY UREA*DAY);

RUN;QUIT;

/*BB. FULL MODEL for Serum; really only intersted in JUN and UREA effected by sampling hour (hour 0 or 6), thus left DAY and PERIOD out of model. Repeated the blood draw on multiple days and periods just to "strengthen the arugument"*/

DATA LAMBserum; SET SERUM;

PROC SORT; BY DAY HOUR ID JUN UREA;

RUN;QUIT;

PROC GLIMMIX;

CLASS DAY PERIOD HOUR ID JUN UREA;

MODEL Glucose = HOUR JUN UREA JUN*UREA JUN*HOUR UREA*HOUR JUN*UREA*HOUR/dist=lognormal ddfm=kr solution; /*QUESTION: do I need DAY and/or PERIOD and it’s interactions in here or is it OK to just average across both periods?*/

Random day/residual subject = ID(PERIOD) type =CSH; /*don’t think I’m specifying this correctly. Do I need "subject=ID(Period)?*/

/*Still need to figure out the contrast statements here since HOUR is now in model"

LSMEANS PERIOD JUN UREA JUN*UREA JUN*PERIOD UREA*PERIOD JUN*UREA*PERIOD/DIFF

ADJUST=SIMULATE (REPORT SEED=121211) cl adjdfe=row SLICE=(PERIOD JUN UREA);

/*"backtransform"*/

ODS OUTPUT lsmeans=lsmeans;

PROC PRINT; RUN; quit;

data btlsmeans;

set lsmeans;

omega=exp(stderr*stderr);

btlsmean=exp(estimate)*sqrt(omega);

btvar=exp(2*estimate)*omega*(omega-1);

btsem=sqrt(btvar);

PROC PRINT; RUN;

Re: PROX MIXED & GLIMMIX Repeated measures, interactions, and linear/quad. contrasts

AgReseach7 — Mon, 17 Oct 2016 23:50:25 GMT

Initial post may have be confusing, so I'll try & narrow it down. What I need help with is setting up a linear & quadratic contrast statement for 4 levels of treatment A (Juniper) x 2 levels of treatment B (urea). When I have this, I'll be able to evalute linear trends in Juniper concentration, at each level of urea (if P-value says there is an interaction). If there isn't an interaction, then I'll just run a contrast for juniper or urea separately.

Class level info.

Juniper 15 30 45 60

Urea 1 3

CONTRAST 'LINEAR' JUN*UREA XXXXXXXXXXXXXX;

CONTRAST 'QUADRADIC' JUN*UREA XXXXXXXXXXXXXX;

If no JUN*UREA interaction, then I'll run the following

CONTRAST 'LINEAR' JUN XXXXXXXXXXXXXX;

CONTRAST 'QUADRADIC' JUN XXXXXXXXXXXXXX;

CONTRAST 'UREA1 VS. UREA3' UREA XXXXXXXXXXXXXX;

Contrasts with an interaction

AgReseach7 — Fri, 21 Oct 2016 01:13:36 GMT

Update. I'm getting closer, but still need help. For example, I'll focus on body weight of lambs. The full model was ran and here are the P-values: day(<0.001); juniper(0.53); Urea(0.06); Jun x urea(0.93); juniper x day(<0.001); urea x day(0.90) JxUxday(0.61)

Thus, for this variable, I need to run linear/quad/cubic contrast statements for JUN x Day. Any idea what I'm doing wrong below?

I discovered the "e" option at the end of the contrast statement, to get contrast coefficients (see below).

However, contrasts are not running correctly (log: .. is not estimatable) & I can't figure it out.

Juniper (15, 30, 45, 60)

Urea (1, 3)

day (0, 5, 12, 19, 26, 33, 40)

DATA LAMB; SET BW;
PROC SORT; BY DAY ID JUN UREA;
PROC PRINT; RUN;QUIT;

PROC MIXED; BY DAY;
CLASS JUN DAY UREA ID;
MODEL BWkg = JUN|DAY|UREA/DDFM=KR;
REPEATED DAY/SUBJECT=ID TYPE = CS;
CONTRAST 'LINEAR JUNIPER'            JUN       -3 -1   1 3/e;
CONTRAST 'QUADRADIC JUNIPER'         JUN        1 -1 -1 1/e;
CONTRAST 'CUBIC JUNIPER'             JUN       -1   3 -3 1/e;
CONTRAST 'LINEAR J*DAY'        JUN*DAY -.429 -.429 -.429 -.429 -.429 -.429 -.429 -.143 -.143 -.143 -.143 -.143 -.143 -.143   .1429 .1429 .1429 .1429 .1429 .1429 .1429 .4286 .4286 .4286 .4286 .4286 .4286 .4286;
CONTRAST 'QUADRADIC J*DAY'     JUN*DAY   .1429 .1429 .1429 .1429 .1429 .1429 .1429 -.143 -.143 -.143 -.143 -.143 -.143 -.143 -.143 -.143 -.143 -.143 -.143 -.143 -.143 .1429 .1429 .1429 .1429 .1429 .1429 .1429;
CONTRAST 'CUBIC J*DAY'         JUN*DAY -.143 -.143 -.143 -.143 -.143 -.143 -.143   .4286 .4286 .4286 .4286 .4286 .4286 .4286 -.429 -.429 -.429 -.429 -.429 -.429 -.429 .1429 .1429 .1429 .1429 .1429 .1429 .1429;
CONTRAST 'UREA'                    UREA        1 -1;
LSMEANS JUN UREA DAY JUN*UREA JUN*DAY UREA*DAY/DIFF ADJUST=SIMULATE (REPORT SEED=121211) cl adjdfe=row SLICE=(DAY JUN UREA);
*/RUN;QUIT;

Re: Contrasts with an interaction

SteveDenham — Wed, 02 Nov 2016 17:00:06 GMT

This is where the LSMESTIMATE statement might be more convenient to use than a CONTRAST statement. I think roundoff errors are killing this. If you stay with the CONTRAST statement, try to use integer values, with a DIVISOR= option to avoid the roundoff/nonestimable thing.

Also someplace you refer to gain/feed, and using an arcsin square root transformation. Since the numerator and denominator are both random normals, the ratio has an ugly distribution. Since I know the value is constained between 0 and 1, I would consider using PROC GLIMMIX with a BETA distribution. I would not construct the rather ungainly (ratio/1000 + 0.5). Transformations are normally meant to stabilize variances as normal. When you have tools that don't require normal distributions, you can let that go, and just analyze GFkgkg.

Steve Denham

Re: Contrasts with an interaction

AgReseach7 — Fri, 18 Nov 2016 16:57:07 GMT

What if the G:F has negative values?

Can you still use beta or do you have to add a constant to all data?

Re: Contrasts with an interaction

SteveDenham — Wed, 30 Nov 2016 20:39:33 GMT

In a word, crap. I forgot sheep lose weight when you feed them shrubs (or ammonia processed rice straw)...

So gain/feed isn't bounded below or above, and is the ratio of two normals. That means that a Cauchy distribution might work, but, of course, GLIMMIX doesn't fit a Cauchy. A T distribution is approximatelya Cauchy (at least it has heavy tails), so you might consider that.

Really, fitting it as a variable with a normally distributed error makes as much sense as anything at this point.

Steve Denham