Multi-year data analysis (with crop rotation) in proc mixed

[Olanike]

I have a 4-year experiment, with four crops in 4 rotations. There were a total of 48 plots, 12 plots per block. Plots with G_Manure are terminated at flowering which means no yield while other crops were harvested and have yield data. I tried to do a multi year analysis, but found oot I was getting yield data for tyhe Green manure (GM). The code I used is below:

Data _2015 ; set Work.FieldData2015;
proc print ;
run;

ods graphics on;
Data AllYield2015;
Set _2015;
proc sort data = AllYield2015;
by Year TillageTrt;
/* Running Factorial Design without Rep effect but with year effect, The Covtest provides the covariance parameter estimates*/

proc mixed data= AllYield2015 covtest;
   class Year Rep T_System;
   model After_Cleaning=T_System;
   random  year Rep(year) year*T_System/solution;
   Lsmeans T_System/diff=all;
run;

and my LS means output is:

  Type 3 Tests of Fixed Effects
 
                                                         Num     Den
                                           Effect         DF      DF    F Value    Pr > F

                                           T_System       13      29       7.30    <.0001


                                                        Least Squares Means
 
                                                                 Standard
                            Effect      T_System     Estimate       Error      DF    t Value    Pr > |t|

                            T_System    HIGH DC_W     2317.75      352.41      29       6.58      <.0001
                            T_System    HIGH DG_M     1024.92      595.48      29       1.72      0.0959
                            T_System    HIGH D_FL     1408.64      415.28      29       3.39      0.0020
                            T_System    HIGH D_PE      965.27      466.90      29       2.07      0.0477
                            T_System    HIGH Dlen      721.62      415.28      29       1.74      0.0929
                            T_System    HIGH Dmus     1155.97      415.28      29       2.78      0.0094
                            T_System    HIGH SC_W     2756.74      352.41      29       7.82      <.0001
                            T_System    LOW DC_W      1761.43      352.41      29       5.00      <.0001
                            T_System    LOW DG_M      1245.06      595.48      29       2.09      0.0454
                            T_System    LOW D_FL      1089.20      415.28      29       2.62      0.0138
                            T_System    LOW D_PE       841.79      468.57      29       1.80      0.0828
                            T_System    LOW Dlen      1088.95      415.28      29       2.62      0.0138
                            T_System    LOW Dmus       849.83      415.28      29       2.05      0.0499
                            T_System    LOW SC_W      2277.74      352.41      29       6.46      <.0001

I belive I should not have yield results for G-M which is my green manure.

ALSO MY DATA IS ATTACHED

Does anyone knows the code for multi year analysis of a rotation experiment?

[SteveDenham]

Try this:

 proc mixed data= AllYield2015 ;
 class Year Rep T_System;
 model After_Cleaning=T_System year year*T_System;
 random rep(year);
 repeated year /subject=rep type=ar(1) solution;
 Lsmeans T_System year year*T_System/diff=all;
 run;

and see if all of the main effect LSmeans are estimable. If so, then forcing the structural zero response for GM is going to be harder.  If the main effects LSMeans are NOT estimable, then this can be approached by using what Milliken and Johnson call a "means mode", with the use of ESTIMATE and LSMESTIMATE statements to get tests and values of interest.

I removed the covtest option as the implementation in PROC MIXED isn't optimal, especially the testing aspect.  I selected ar(1) as the covariance error structure, which also assumes that the rep is the subject measured from year to year.   

Missing from this appears to be the blocking factor.

Anyway, this should get you started on the repeated measures analysis.

Steve Denham

[Olanike]

Hi Steve, Thanks for the response. I tried the code but I got the error below:

repeated year/subject=rep type=ar(1) solution;
________
22
202
ERROR 22-322: Syntax error, expecting one of the following: ;, GROUP, HLM, HLPS, LDATA, LOCAL, LOCALVAR, LOCALW, NONLOCALW, R, RC,
RCI, RCORR, RI, SSCP, SUBJECT, TYPE.
I did not know what was really missing from the code.

Thanks.

[SteveDenham]

This is my error.  Remove the 'solution' option.  It is not available on the REPEATED statement.  Replace it with the R or RCORR option.

Steve Denham

[Olanike]

Thanks Steve, it worked, the log was clean (no error). However, no values
were displayed for the LSmeans, I will like to make contrasts statements
based on the LSmeans.

##- Please type your reply above this line. Simple formatting, no
attachments. -##

[SteveDenham]

You may want to attach your output listing.  It makes a difference as to whether the LSmeans are non-estimable, or just not presented because PROC MIXED failed to converge.

Steve Denham

[Olanike]

Hi Steve,

My log info is below:

There were 192 observations read from the data set WORK.ALLYIELD2015.
NOTE: The data set WORK.ALLYIELD2015 has 192 observations and 14 variables.
NOTE: PROCEDURE SORT used (Total process time):
      real time           0.01 seconds
      cpu time            0.01 seconds
      

24         proc mixed data= AllYield2015;
25            class Year Rep T_System;
26            model After_Cleaning = T_System year year*T_System;
27            random Rep(year);
28            repeated year/subject=rep type=ar(1) RCORR;
29            Lsmeans T_System year year*T_System/diff=all;
30            run;

NOTE: 1 observation is not included because of missing values.
NOTE: An infinite likelihood is assumed in iteration 0 because of a nonpositive definite estimated R matrix for Rep 1.
NOTE: The PROCEDURE MIXED printed page 2.
NOTE: PROCEDURE MIXED used (Total process time):
      real time           0.03 seconds
      cpu time            0.03 seconds
      

      endsas;

AND IN THE OUTPUT (LST) IT GAVE ME MY DATA (because of proc print) WHICH IS COMPLETELY FINE AND THEN THE STATEMENTS BELOW:

                                       Data Set                     WORK.ALLYIELD2015        
                                       Dependent Variable           After_Cleaning           
                                       Covariance Structures        Variance Components,     
                                                                    Autoregressive           
                                       Subject Effect               Rep                      
                                       Estimation Method            REML                     
                                       Residual Variance Method     Profile                  
                                       Fixed Effects SE Method      Model-Based              
                                       Degrees of Freedom Method    Containment              


                                                      Class Level Information
 
                                        Class       Levels    Values

                                        Year             6    2010 2011 2012 2013 2014 2015
                                        Rep              4    1 2 3 4                       
                                        T_System        12    HIGH DC_W HIGH D_FL HIGH D_PE
                                                              HIGH Dlen HIGH Dmus HIGH SC_W
                                                              LOW DC_W LOW D_FL LOW D_PE LOW
                                                              Dlen LOW Dmus LOW SC_W        


                                                            Dimensions

                                                Covariance Parameters             3
                                                Columns in X                     67
                                                Columns in Z                     24
                                                Subjects                          1
                                                Max Obs Per Subject             192


                                                      Number of Observations

                                            Number of Observations Read             192
                                            Number of Observations Used             191
                                            Number of Observations Not Used           1

[SteveDenham]

I will bet that the subject needs to be more correctly specified--this error is associated with multiple observations per time point for the subject which here is rep.  I will wager it needs to be rep*T_system.  Once you get that straightened out, everything should work.

Steve Denham

[Olanike]

Thanks Steve. I think I am still having difficulty with the code. I specified the subject to be rep*T_System, but my log gave the error below:

24         proc mixed data= AllYield2015;
NOTE: PROCEDURE MIXED used (Total process time):
      real time           0.01 seconds
      cpu time            0.01 seconds
      
NOTE: The SAS System stopped processing this step because of errors.
25            class Year Rep T_System;
26            model After_Cleaning = T_System year year*T_System;
27            random Rep(year);
28            repeated year/rep*T_System=rep type=ar(1) RCORR;
                            ___
                            1
                               _
                               22
                               76
WARNING 1-322: Assuming the symbol R was misspelled as rep.
ERROR 22-322: Syntax error, expecting one of the following: ;, =, GROUP, HLM, HLPS, LDATA, LOCAL, LOCALVAR, LOCALW, NONLOCALW, R,
              RC, RCI, RCORR, RI, SSCP, SUBJECT, TYPE.  
ERROR 76-322: Syntax error, statement will be ignored.
29            Lsmeans T_System year year*T_System/diff=all;
30            run;
31            endsas;

I WANTED TO USE REPEATED MEASURE BECAUSE THE EXPERIMENT WAS REPEATED 6 TIMES (I MEAN 6 YEARS) , ALTHOUGH ALL THE PARAMETERS (SEEDLING DENSITY, YIELD ETC) WERE MEASURED ONCE PER GROWING PERIOD. IS THERE ANY OTHER WAY AROUND THIS OR REPEATED MEASURE IS THE BETTER.

ALSO, IS THERE A WAY FOR ME TO FIND THE YEARLY DIFFERENCE

THANKS

[SteveDenham · Accepted solution]

Line 28 from the log says that the code you are running needs to be fixed up, and I apologize if my code wasn't clear.

proc mixed data= AllYield2015 ;
 class Year Rep T_System;
 model After_Cleaning=T_System year year*T_System;
 random rep(year);
 repeated year /subject=rep*T_system type=ar(1) solution;
/* I suspect that you have identically numbered reps in the various T-Systems */
 Lsmeans T_System year year*T_System/diff=all;
 run;

Try the above.

Steve Denham 

[Olanike]

 Hi Steve,

It finally worked. Thanks for the clarification. However, some parameters were non-estimable, I was thinking maybe there was not enough degree of freedom - PLEASE FIND ATTACHED.

 IS THERE A WAY FOR ME TO FIND THE YEARLY DIFFERENCE?.  ALSO IS THERE A WAY TO SPECIFY MY ALPHA TO BE 10% INSTEAD OF THE CONSTANT 5%

THANKS

[Olanike]

 IS THERE A WAY FOR ME TO FIND THE YEARLY DIFFERENCE?.  ALSO IS THERE A WAY TO SPECIFY MY ALPHA TO BE 10% INSTEAD OF THE CONSTANT 5%

THANKS

[Olanike]

Hi Steve,

When I used repeated measures statement, my log is clean free of error but some of the LSmeans for T_System and years are non-estimable.

If I do "Lsmeans T_System" WITHOUT YEAR AND YEAR*T_SYSTEM, all parameters are estimable but I don't get the yearly difference.

Is there a statement to find the yearly difference (e.g High DC_W 2010 vs High DC_W 2011) without getting a long list of Lsmeans of all the years.

Also some of the plots have crops that are different every other year (e.g Flax in 2010, wheat in 2011 and Flax in 2012) , is there a way to get the yearly difference or effect?

Thanks.

[Olanike]

Hi Steve,
The code above worked perfectly well but I have a very long list of LSmeans for years and year*T_System. When I use the code below:
proc mixed data= AllYield2015T;
class Year Rep T_System;
model &var = T_System;
random Rep(year)year year*T_System/solution;
repeated year/subject=rep*T_System type=ar(1) RCORR;
Lsmeans T_System/diff=all;
IT WORKS FINE AS WELL AS WELL BUT I DO NOT HAVE THE YEARLY DIFFERENCE. PLEASE HOW DO I COMPARE TWO YEARS FOR INSTANCE HIGH DC_W 2010 VS DC_W 2011, 2012 AND 2013.

THANKS.