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).
Objectives:
Proposed plan for the Growth data (blood serum similar):
Questions:
/*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;
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
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;
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;
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
What if the G:F has negative values?
Can you still use beta or do you have to add a constant to all data?
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
Are you ready for the spotlight? We're accepting content ideas for SAS Innovate 2025 to be held May 6-9 in Orlando, FL. The call is open until September 25. Read more here about why you should contribute and what is in it for you!
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.