Turn on suggestions

Auto-suggest helps you quickly narrow down your search results by suggesting possible matches as you type.

Showing results for

- Home
- /
- Programming
- /
- Graphics
- /
- Re: ANOVA using SAS 9.3 Proc Mixed

Options

- RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Mute
- Printer Friendly Page

- Mark as New
- Bookmark
- Subscribe
- Mute
- RSS Feed
- Permalink
- Report Inappropriate Content

Posted 02-11-2014 04:17 PM
(2706 views)

I am doing ANOVA using SAS 9.3 Proc Mixed on a field data (RCBD with four treatments (One control vs. three different land uses). The response variable is water table depth (WTD) which is collected every 15-minute. I would like to know which of the two sets of codes is more appropriate to use. I would also welcome any suggestion to improve the set of codes. Seas = Seasons (traditional: winter, spring, summer, and fall), Day = Julian Day, WTD = water table depth averaged over a day. Here are the codes:

**A)**

ODS GRAPHICS ON;

PROC MIXED DATA=WORK.WTDWOR_2014 plots(only)=(ResidualPanel(marginal))

PLOTS(MAXPOINTS= 2000000);

CLASS Block Treat Year Seas Day;

MODEL WTD = Treat|Year|Seas/ddfm=satterthwaite residual;

RANDOM Block Block*Treat Block*Treat*Year Block*Treat*Year*Seas;

LSMEANS Treat Treat*Year Treat*Year*Seas/adjust = tukey;

RUN;

ODS GRAPHICS OFF;

**B)**

ODS GRAPHICS ON;

PROC MIXED DATA=WORK.WTDWOR_2014 plots(only)=(ResidualPanel(marginal))

PLOTS(MAXPOINTS= 2000000);

CLASS Block Treat Year Seas Day;

MODEL WTD = Treat|Year|Seas Seas(Year) Treat*Seas(Year)/ddfm=satterthwaite residual;

RANDOM Block Block*Treat Block*Treat*Year Block*Treat*Year*Seas;

LSMEANS Treat Treat*Year Treat*Seas(Year)/adjust = tukey;

RUN;

ODS GRAPHICS OFF;

P.S. On set of codes B, I was trying to let SAS know that season(Seas) is embedded within a year. Is this correct?

Thanks.

9 REPLIES 9

- Mark as New
- Bookmark
- Subscribe
- Mute
- RSS Feed
- Permalink
- Report Inappropriate Content

- Mark as New
- Bookmark
- Subscribe
- Mute
- RSS Feed
- Permalink
- Report Inappropriate Content

Thanks Rick...

- Mark as New
- Bookmark
- Subscribe
- Mute
- RSS Feed
- Permalink
- Report Inappropriate Content

Some interesting points to consider here:

1. You include DAY as a class variable, but do not use it in the model. As the analysis now stands, you are averaging over all days, even though they are in different years and seasons. This may lead to problems.

2. Of the two models presented, I would think they are identical. Your code for the LSMEANS is certainly appropriate in either case, given how SAS models nested factors (see the documentation on this) SEAS(YEAR) and SEAS*YEAR set up exactly the same matrix.

3. I guess my main concern is that this is a repeated measures design of some complexity, and the correlations within a plot over time are not being modeled. It may be that this just can't be done easily. You may wish to think of season-day, if the measurements are spaced out at the same times within each season for each year. However, I have a strong hunch that this level of complexity may not be amenable to modeling, due to memory and other constraints.

Let us know what you think.

Steve Denham

- Mark as New
- Bookmark
- Subscribe
- Mute
- RSS Feed
- Permalink
- Report Inappropriate Content

Thanks Steve. As you may have recognized, I am collecting the data every 15 minutes. I am both interested in main effects and simple effects (of treatments on WTD on annual basis, and on seasonal basis within a given year, if any). I average 240 samples of 15-minute WTD data to come up with an average for a day of a given season per treatment. The raw data are already on in terms of average per day and not the 15-minute interval. Analyzing by season (traditional 4 seasons) means averaging 120 days of WTD into one value per treatment. My point is, I could add a repeated statement in the model; however, the fact that I am averaging the data from 15-minute interval (240) to daily (120) to season makes the autocorrelation issue less of an issue.

You are right with point number 2. I looked at the covariance parameters estimates and fit statistics (e.g. AIC values), the two models were very similar. I just would like to tell the model that DAY is nested within SEAS, and SEAS within Year and I am not sure if I am doing it right. Based on your comment about similarity in matrices, I think I am close.

You mentioned about DAY. Any suggestion as to how to incorporate that into the model (DAY as nested into SEAS) as a code?

Again, thanks for your time.

- Mark as New
- Bookmark
- Subscribe
- Mute
- RSS Feed
- Permalink
- Report Inappropriate Content

TD,

The nested terms in your model are already included when you 'pipe' the fixed effects, so I am now in favor of Model A (* foolish* consistency is the hobgoblin of small minds).

If you have the computing power to handle DAY as a repeated measure, I would try it. How many days of data do you have?

Things to consider:

I would shift to HPMIXED, and then port to GLIMMIX for comparisons of lsmeans (see EX. 46.3 in the HPMIXED documentation)

I would index the days within season, rather than using the Julian date. If they don't align, then forget using a repeated measures approach. If they do align, or would with only minor shifts, then code that looks like:

CLASS Treat Year Seas Day Block Plotid;

MODEL WTD = Treat|Year|Seas|Day ;

RANDOM Block Block*Treat Block*Treat*Year Block*Treat*Year*Seas;

REPEATED Day/subject=plotid(block treat year seas) type=AR(1);

Hope this helps.

Steve Denham

Get everything in ODS datasets and then port it all into GLIMMIX.

Steve Denham

- Mark as New
- Bookmark
- Subscribe
- Mute
- RSS Feed
- Permalink
- Report Inappropriate Content

I have 92 days in 2009; 365 x 3 = 1095 days in 2010, 2011, and 2013; and 366 in 2012. So in each treatment I have 1553 days of data. By the way, I didn't include Julian date in my dataset but rather Julian Day. In short, I have a column for year, separate column for season, and day. For example in 2010, I have column/variable for year which is 2010, another column for season (winter, spring, summer, and fall) and another season for Day (1 to 365).

Thanks again.

- Mark as New
- Bookmark
- Subscribe
- Mute
- RSS Feed
- Permalink
- Report Inappropriate Content

Given that data set up, and the already "binned' nature of the data, would yet another binning by week within season make sense? Then you would have 13 repeated measures on each combination of year, season, block and treatment. If you believe that this is something that reflects an interesting approach, I think it would be the best way to incorporate changes over time. The problem I have with both season and day in the model is that the number of observations per season (and per year) isn't consistent because 4 doesn't divide 365 or 366 equally.

An alternative would be to fit day rather than season, and then address season effects through the use of LSMESTIMATE or CONTRAST statements, viewing the season effect as a composite of the day effect. It will require HPMIXED, followed by GLIMMIX to get this to work though.

Steve Denham

- Mark as New
- Bookmark
- Subscribe
- Mute
- RSS Feed
- Permalink
- Report Inappropriate Content

Hi Steve:

I followed your latest advice (add week, analyze with HPMIXED then port to GLIMMIX). I seemed to have successfully run HPMIXED, but not with GLIMMIX. The error was about PARMSDATA which I highlighted below. Any suggestion how to fix this error? Pardon the long post. I just wanted to give you as much background information as I can.

Thanks.

TD

55 PROC HPMIXED DATA = WORK.WTDWOR_Final2014; *plots(only)=(ResidualPanel(marginal)

NOTE: Writing HTML Body file: sashtml1.htm

56 StudentPanel(conditional)) PLOTS(MAXPOINTS= 2000000);

57 CLASS Block Treat Year Seas Week Day;

58 MODEL WTD = Treat|Year|Seas|Week;

59 RANDOM INT Treat/SUBJECT=Block;

60 REPEATED DAY/SUBJECT=Block*Treat*Year*Seas*Week type=ar(1) R;

61 TEST Treat|Year|Seas|Week;

62 ODS OUTPUT COVPARMS=HPMEstimate;

63 RUN;

NOTE: 129 observations are excluded because of: missing response values (n=129).

NOTE: Convergence criterion (GCONV=1E-8) satisfied.

NOTE: The data set WORK.HPMESTIMATE has 4 observations and 3 variables.

NOTE: PROCEDURE HPMIXED used (Total process time):

real time 2:32:48.21

cpu time 1:24:56.25

64

65 ODS GRAPHICS ON;

66 PROC GLIMMIX DATA = WORK.WTDWOR_Final2014;

67 CLASS Block Treat Year Seas Week Day;

68 MODEL WTD = Treat|Year|Seas|Week;

69 RANDOM INT Treat/SUBJECT=Block;

70 PARMS/PDATA=HPMEstimate hold=1,2 noiter;

71 LSMEANS Treat|Year|Seas|Week/ADJUST=Tukey;

72 RUN;

NOTE: Some observations are not used in the analysis because of: missing response values (n=129).**ERROR: The PARMSDATA= data set must have either 2 more or 1 less observations.**

NOTE: The SAS System stopped processing this step because of errors.

NOTE: PROCEDURE GLIMMIX used (Total process time):

real time 40.35 seconds

cpu time 37.02 seconds

- Mark as New
- Bookmark
- Subscribe
- Mute
- RSS Feed
- Permalink
- Report Inappropriate Content

It looks to me like you aren't fitting the same model in GLIMMIX. Try adding a second random statement

RANDOM DAY/residual SUBJECT=Block*Treat*Year*Seas*Week type=ar(1) R;

and changing the PARMS statement to:

PARMS/PDATA=HPMEstimate hold=1,2,3,4 noiter;

Steve Denham

**Available on demand!**

Missed SAS Innovate Las Vegas? Watch all the action for free! View the keynotes, general sessions and 22 breakouts on demand.

How to Concatenate Values

Learn how use the CAT functions in SAS to join values from multiple variables into a single value.

Find more tutorials on the SAS Users YouTube channel.