Turn on suggestions

Auto-suggest helps you quickly narrow down your search results by suggesting possible matches as you type.

Showing results for

- Home
- /
- Analytics
- /
- Stat Procs
- /
- PROX MIXED & GLIMMIX Repeated measures, interactions, and linear/quad....

Options

- RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Mute
- Printer Friendly Page

🔒 This topic is **solved** and **locked**.
Need further help from the community? Please
sign in and ask a **new** question.

- Mark as New
- Bookmark
- Subscribe
- Mute
- RSS Feed
- Permalink
- Report Inappropriate Content

Posted 10-10-2016 06:43 PM
(1991 views)

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:
- if DAYxJUNxUREA and JUNxUREA are
**significant**: run linear/quad. contrasts for the 4 JUN levels at each DAY & for each level of UREA. - 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 - if DAYxJUNxUREA
**not significant**and JUNxUREA**significant**: run juniper contrasts (averaged across days) at each level of UREA. - if DAYxJUNxUREA and JUNxUREA
**not significant**: run juniper contrasts (averaged across urea and days), and also a contrast for UREA1 vs. UREA3

- if DAYxJUNxUREA and JUNxUREA are
- 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**;

1 ACCEPTED SOLUTION

Accepted Solutions

- Mark as New
- Bookmark
- Subscribe
- Mute
- RSS Feed
- Permalink
- Report Inappropriate Content

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

5 REPLIES 5

- Mark as New
- Bookmark
- Subscribe
- Mute
- RSS Feed
- Permalink
- Report Inappropriate Content

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;

- Mark as New
- Bookmark
- Subscribe
- Mute
- RSS Feed
- Permalink
- Report Inappropriate Content

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;

- Mark as New
- Bookmark
- Subscribe
- Mute
- RSS Feed
- Permalink
- Report Inappropriate Content

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

- Mark as New
- Bookmark
- Subscribe
- Mute
- RSS Feed
- Permalink
- Report Inappropriate Content

What if the G:F has negative values?

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

- Mark as New
- Bookmark
- Subscribe
- Mute
- RSS Feed
- Permalink
- Report Inappropriate Content

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

**Available on demand!**

Missed SAS Innovate Las Vegas? Watch all the action for free! View the keynotes, general sessions and 22 breakouts on demand.

What is ANOVA?

ANOVA, or Analysis Of Variance, is used to compare the averages or means of two or more populations to better understand how they differ. Watch this tutorial for more.

Find more tutorials on the SAS Users YouTube channel.