10-10-2016 06:43 PM
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).
Proposed plan for the Growth data (blood serum similar):
/*AA: FULL MODEL*/
OPTIONS PAGESIZE=60 LINESIZE=90;
DATA LAMBGROW; SET GROW;
PROC SORT; BY DAY ID JUN UREA;
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);
/*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;
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);
ODS OUTPUT lsmeans=lsmeans;
PROC PRINT; RUN; quit;
PROC PRINT; RUN;
10-17-2016 07:50 PM
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;
10-20-2016 09:13 PM
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);
11-02-2016 12:53 PM - edited 11-02-2016 01:00 PM
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.
11-30-2016 03:39 PM
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.