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