BookmarkSubscribeRSS Feed
☑ This topic is solved. Need further help from the community? Please sign in and ask a new question.
Laser_Taco_
Fluorite | Level 6

Hello, 

Exp details: 24 pens, 3 feed types (800, 400, 100), 8 replicates, RCBD

 

The purpose is to evaluate a preference for feed type (800, 400, 100).

To make the comparisons, three "treatments" were created (A = 800 vs 100, B= 800 vs 400, C =400 vs 100), the A,B, C treatments were used to allocate 2 of the 3 feeds to 1 pen in a RCBD.

 

For 5 days, a known amount of feed was added to each feeder/pen, and after 24 h the feed was weighed back to determine the feed disappearance.  Preference was defined as  ((disappearance of ___ feed (800, 400, 100; respectively))/ Total feed disappearance (combined total of the 2 feeders/pen)*100). 

Giving the proportion of of each feed /pen/day.

 

In the proposal (was proved to me...), the Statistical Analysis is suppose to be "...preference test as the fixed effect and replicate as the random effect. LSmean and means are to be separated using pdiff with Tukey's adjustment with pen as the EU." IMHO this somewhat lacking and not very descriptive. 

 

CODE: 

input

PENDIETTreatmentRepFeedpref1pref2pref3pref4pref5
113180068.750057.647167.857157.954578.9157
133110031.250042.352932.142942.045521.0843
211180079.710144.736844.444441.304348.3333
221140020.289955.263255.555658.695751.6667
3221400100.000067.567685.185280.582556.0847
33211000.000032.432414.814819.417543.9153
411280092.592657.746565.517273.913059.0278
4212

400

7.407442.253534.482826.087040.9722

.... 

;
Run;

Proc Print;
Run;

%macro model_loop;
%let yvar1 = pref1;
%let yvar2 = pref2;
%let yvar3 = pref3;
%let yvar4 = pref4;
%let yvar5 = pref5;

%do i=1 %to 5 ;

 

 

Proc Mixed data=Exp200 NOBOUND ASYCOV method=REML;
Title "&&yvar&i";
class Pen Trt Rep Feed;
model &&yvar&i = Feed/ddfm=Kenwardroger;
random Rep / subject=Pen;
lsmeans Feed/pdiff adjust=tukey;
store out=work.glm&&yvar&i;
run;

proc plm restore=work.glm&&yvar&i;
Title "&&yvar&i PLM";
lsmeans Feed/ pdiff lines adjust=tukey;
run;

%end;
%mend model_loop;
%model_loop;
quit;

 

 

The issue I'm having is a failure to converge . I have a feeling that it has to do with the specification of the subject. 

 

I guess I don't have a super specific question. I just don't feel like this is the correct way to analyze the data and looking for advice. 

Regards,

 

 

 

 

1 ACCEPTED SOLUTION

Accepted Solutions
SteveDenham
Jade | Level 19

Rather than a ratio of each component to the total you might consider a ratio of consumption of the feed with the larger feed level to the smaller feed level for each pen. You get a single value that is indexed by the variable 'Treatment'. Or consider the difference in consumption calculated in much the same way. That might even be better as it is much more likely to collectively have residuals that are distributed normally. The lsmeans by treatment are then a measure of how much more was consumed of the higher valued feed - which is a clear measure of preference.

 

SteveDenham 

View solution in original post

9 REPLIES 9
jiltao
SAS Super FREQ

It seems to me that for every Pen, you have two observations. The values for these two observations for each of Pref1-Pref5 always add up to 100. That might be a problem. You might want to reconsider how you want to arrange your data and set up your model.

Thanks,

Jill

Laser_Taco_
Fluorite | Level 6
Hi Jill,

Thank you for the response. I also think that is the issue.

Since these values are percentages that add to 100% (when combined), they aren't independent.

Any idea or suggestion on how I can re-arrange?

Regards,
Kyle
jiltao
SAS Super FREQ

One approach might be to just use one observation for each Pen.

 

Laser_Taco_
Fluorite | Level 6
Thank you again for the response. How might only using 1 or the observations impact the interpretation or the results?
Regards,

SteveDenham
Jade | Level 19

I hope I understand this correctly. Each pen is randomly assigned to one of 3 treatments (A, B, C). Thus, a pen will have one value. I don't see the advantage of calculating a value for each as a proportion. Have you considered fitting the raw consumption values (which should be independent), and then calculating the ratios from the least squares means?

 

The other point to consider is that ratios probably do not yield normally distributed residuals, as the cut-offs at 0 and 1 may be a problem. If you are going to analyze the proportions as the dependent variable, a beta distribution for the values may be more appropriate. That would mean shifting to GLIMMIX or GEE.

 

And last, failure to converge. That very likely can be fixed by adding MAXITER= and MAXFUNC= options to the PROC MIXED statement, provided there are no other issues. Check both the .log and .lst file for notes, warnings and notifications that indicate that the algorithm has difficulties other than slow convergence.

 

SteveDenham

Laser_Taco_
Fluorite | Level 6
Hello @stevendenham,
Thank you for your response!

Yes, you are correct. There are 3 pens/block (i used the word "Rep" in my initial post) with 8 blocks in total.

The reason I have two values for each pen is due to the different consumptions from the two feeders in each pen. If I were to use only 1 value from the pen, I am unsure how I would differentiate the consumption between the two feeders within the pen (800 v 100 etc.).

To your point on proportions, I also considered this and agree especially with a handful of cases were the consumption for the feeders within a pen are 0 and 1, respectively. Admittedly, I am a bit unfamiliar with GLIMMIX and GEE, but should the code look something like this:

Proc GLIMMIX data=Exp200 ;
class Pen Trt Block Feed;
model y= Feed / ddfm=KR;
random block/ subject=Pen; **Unsure if " /subject=pen" is necessary**
lsmeans Feed/pdiff lines adjust=tukey;
run;

As for the failure to converge, I was able to get it to converge (with no errors) by removing the NOBOUND option, but the Cov Parm Estimate for Block/subject=pen was zero which is odd.

Regards,
Kyle



SteveDenham
Jade | Level 19

hi Kyle @Laser_Taco_ ,

 

Just a few things that occurred to me over the weekend. In GLIMMIX, for ratios you should probably specify a distribution, and change the method to LAPLACE with tech=nrridg. So perhaps something like this:

 

Proc GLIMMIX data=Long_version method=laplace;
nloptions maxiter=1000 tech=nrridg; class Pen Trt diet Block Feed instance; model y=Trt Feed instance Feed*instance/dist=beta ddfm=KR; random instance/subject=pen*diet type=ar(1); random feed/subject=pen type=chol; random block; lsmeans Feed/pdiff lines adjust=tukey; run;

The first thing to do is to re-arrange the data into a long format, where there is one dependent value in each observation (here it is y). The other variables in each line would be pen, trt, diet (although this looks completely confounded with feed), block, feed and instance, which is indexed from 1 to 5. So instead of:

PEN DIET Treatment Rep Feed pref1 pref2 pref3 pref4 pref5
1 1 3 1 800 68.7500 57.6471 67.8571 57.9545 78.9157
1 3 3 1 100 31.2500 42.3529 32.1429 42.0455 21.0843

 

You would have (for the first pen, and all other pens would follow this pattern:

PEN DIET Treatment Rep Feed instance y
1 1 3 1 800 1 0.6875
1 1 3 1 800 2 0.576471
1 1 3 1 800 3 0.678571
1 1 3 1 800 4 0.579545
1 1 3 1 800 5 0.789157
1 3 3 1 100 1 0.3125
1 3 3 1 100 2 0.423529
1 3 3 1 100 3 0.321429
1 3 3 1 100 4 0.420455
1 3 3 1 100 5 0.210843

 And so on. Note that all of the y values have been rescaled to the (0, 1) interval. So now there are two levels of feed for each pen and five levels of instance for each pen by diet combination.

 

If I were better informed about how to use it, I might recommend PROC COPULA for the analysis, but that is only based on what the examples look like.

I thought I had a way to address this, but so long as each pref within a treatment is such that the values for the two diets within a treatment sum to 1, all of the above will have some major issues due to lack of independence. If the preference response was converted to a ratio or some other type of value, so that a single value is obtained per pen at each preference measurement (time), you should be able to improve the analysis.

 

SteveDenham

 

 

Laser_Taco_
Fluorite | Level 6

Hi @SteveDenham ,

 

I'm not sure if I understand what you mean by "converted to a ratio or some other type of value, so that a single value is obtained per pen at each preference measurement (time), you should be able to improve the analysis.".  

I understand that having multiple responses from the same pen is my issue here. For converting the values to a ratio, each value per diet/pen is already a ratio == Diet__ :Total Consumed. 

 

Regards,

SteveDenham
Jade | Level 19

Rather than a ratio of each component to the total you might consider a ratio of consumption of the feed with the larger feed level to the smaller feed level for each pen. You get a single value that is indexed by the variable 'Treatment'. Or consider the difference in consumption calculated in much the same way. That might even be better as it is much more likely to collectively have residuals that are distributed normally. The lsmeans by treatment are then a measure of how much more was consumed of the higher valued feed - which is a clear measure of preference.

 

SteveDenham 

SAS Innovate 2025: Save the Date

 SAS Innovate 2025 is scheduled for May 6-9 in Orlando, FL. Sign up to be first to learn about the agenda and registration!

Save the date!

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.

Discussion stats
  • 9 replies
  • 2427 views
  • 3 likes
  • 3 in conversation