topic Re: Split-Plot design replicated in time in Statistical Procedures

Split-Plot design replicated in time

Ana79 — Mon, 29 Jul 2013 16:52:06 GMT

I want to analyze a data set of soybean yield from a split plot RCBD that was conducted over 2 years.

The fixed factors are irrigation (whole plot, 2 levels) and fertilizer (subplot 3 levels)

Random factors are blocks (3) and years (2).

I want to test whether the years are different (so I need to address each year separately) or I can pool the years to get a more robust understanding of both water and fertilizer treatments.

So, I'm trying to use proc mixed, but I'm not sure I am testing year right (should it be nested or not?):

proc mixed data=soy;

class year blk irr fert;

model Y = irr fert irr*fert;

random intercept blk blk*irr / subject=year;

Thanks

Re: Split-Plot design replicated in time

SteveDenham — Mon, 29 Jul 2013 18:13:56 GMT

This fits a split plot, with three random effects - year, year by blk, and year by blk by irr. Everything should work out, except for one thing. You say that you want to test whether the years are different. The only way to do this would be to shift year from a random effect to a repeated effect, and include it, and its interactions, in the model statement. The following should work:

proc mixed data=soy;

class year blk irr fert;

model y = irr fert irr*fert year year*irr year*fert year*irr*fert;

random intercept irr/subject=blk;

repeated year/type=un subject=blk*irr*fert;

lsmeans irr fert irr*fert year year*irr year*fert year*irr*fert;

run;

This assumes that the same blk*irr*fert combination is measured in successive years.

A good reference for this sort of model is Littell et al. SAS System for Mixed Models, 2nd ed.

Steve Denham

Re: Split-Plot design replicated in time

Ana79 — Mon, 29 Jul 2013 21:31:41 GMT

Thanks Steve.

The blk*irr*fert combination is not the same in the sense that the experiment was moved to a new location in the 2nd year. So, I can not test years.

An option may be to quantify the amount of the total random variation in the data associated with the random year effect from the covariance parameter estimates,

So, do I need to modify the random statement to include all possible combinations? Do I need to nest blocks within years?

proc mixed data=soy covtest;

class year blk irr fert;

model y = irr fert irr*fert;

random year blk year*blk year*irr blk*irr year*irr(block);

lsmeans irr fert irr*fert;

run;

Re: Split-Plot design replicated in time

SteveDenham — Tue, 30 Jul 2013 12:09:50 GMT

Nesting within years actually results in the same matrix as crossing with year (see Nested Effects in the Details section of the PROC MIXED documentation). Given that, your code should be good. I would be very hesitant, however, to use the covtest results to check for significance of any of the variance components. Covtest uses a Wald chi-square which really needs a large sample size to be accurate for tests on variances. Instead, fit "reduced" models, and consider likelihood ratio tests or information criteria (such as AICc) as a way of selecting various random parts of the model If you do go the way of information criteria, be sure that the fixed effects remain constant. Also, deleting random effects in a split-plot may result in changes in the denominator degrees of freedom that upset the testing of the fixed effects.

Steve Denham

Re: Split-Plot design replicated in time

Ana79 — Tue, 30 Jul 2013 18:43:22 GMT

Thanks Steve.

I used:

proc mixed data=soy2;

class year bl water trt;

model RTO = water trt water*trt;

random year bl year*bl year*water bl*water year*water*bl;

lsmeans water trt water*trt;

run;

I run it but I get zero variance estimation.

year	1.17E-12
Bl	0
year*Bl	0
year*water	0
Bl*water	0
yearBlwater	26056
Residual	84540

I read the multi location experiment in Ch10 (Littel et al) that uses BLUPs and the split plot RCBD in Ch 2.

I'm wondering what mess I'm doing with the error terms

From the chapters:

For the split plot the variance is composed of BL, Whole-Plot error (Bl*water), and Split-Plot error (residual)

For a multi location RCBD (no split plot) the variance would be composed of: location (year in my case), rep(location) (bl(year) in my case), and location*trt (year*water*fert in my case, right?)

So, adding years to the split plot would be like adding another blocking term (to the WP and SP), but I'm not so sure the way I did it is the right one.

I think the variance would be made up by :

year

bl(year)

WP error (bl*water(year))

SP error (residual)

proc mixed data=soy2;

class year bl water trt;

model RTO = water trt water*trt;

random year bl(year) bl*water(year);

lsmeans water trt water*trt;

run;

But using this I get the following estimates:

year	0
Bl(year)	0
Bl*water(year)	26057
Residual	84540

So, not sure if I am missing any random error term or there is no variance among years or no variance for bl(year).

Any insight what can be going on?

Re: Split-Plot design replicated in time

SteveDenham — Wed, 31 Jul 2013 12:05:14 GMT

I think what is going on is "small data set". After fitting the variance component for the three-way term, there really isn't any variability left to explain by BI(year) or year. Those zeroes are the result of using REML. Someplace in the output there is probably a statement that the G matrix is not positive definite, which is a note that some variance component has been set to zero in the REML algorithm. Do not despair--the results are correct, and the degrees of freedom probably are as well. However, I would consider adding a Kenward-Rogers adjustment for the denominator degrees of freedom, to get a better estimate of the F probability and of standard errors of the fixed effects, including the lsmeans. Add /ddfm=kenwardrogers to the model statement, and see what effect this has.

Steve Denham

Re: Split-Plot design replicated in time

Ana79 — Wed, 31 Jul 2013 20:05:07 GMT

Thanks Steve,

Yes, The log says the G matrix is not positive definite. I used the ddfm=kenwardroger and the denominator df for the fixed effects change (for water goes from 5 to 10) but it does not change the zero variance for year and Bl(year) (still zeroes).

I tried ddfm=SATTERTHWAITE too and I get same results. So I decided not to adjust by df for now (I know I will need it for the lsmeans).

So, to test the Ho: WP error= 0, or in other word if the water treatment by year variance was substantial, so I followed the steps in Ch2:

1) run the full model

title 'full Yield';

proc mixed data=soy;

class year bl water trt;

model RTO = water trt water*trt;

random year bl(year) Bl*water(year);

run;

year	0
bl(year)	0
Bl*water(year)	27908
Residual	50356

-2 Res Log Likelihood	430.5
AIC (smaller is better)	434.5
AICC (smaller is better)	434.9
BIC (smaller is better)	431.9

2) run the reduced model without the WP error

title 'reduced wp Yield';

proc mixed data=soy3;

class year bl water trt;

model RTO = water trt water*trt;

random year bl(year);

run;

year	0
bl(year)	0
Residual	78272

-2 Res Log Likelihood	433.9
AIC (smaller is better)	435.9
AICC (smaller is better)	436.1
BIC (smaller is better)	434.6

3) run the reduced model without the year and bl(year):

title 'reduced Yield';

proc mixed data=soy3;

class year bl water trt;

model RTO = water trt water*trt;

run;

Residual	78272

-2 Res Log Likelihood	433.9
AIC (smaller is better)	435.9
AICC (smaller is better)	436.1
BIC (smaller is better)	437.3

4) calculate the chi-square statistic using the -2 Res Log Likelihood

for WP error Chi-square= 433.9-430.5= 3.4 so probability with 1 df is 0.0651, and p-value is 0.0325.

for year (and bl(year)) Chi-square= 433.9-433.9=0 so does not have a substantial contribution to the variance (which makes sense from the estimation of variance components table)

So, will be right to say that there is a significant variation water*year, and thus the effect of year on how water trt influenced yield need to be assessed using BLUPs?.

Re: Split-Plot design replicated in time

SteveDenham — Thu, 01 Aug 2013 16:52:49 GMT

Point one: Use the kenward-rogers adjustment in small sample experiments. It computes a Satterthwaite degrees of freedom adjustment, and a shrinkage for the standard errors that reflects redundance (hidden correlation) in the design. It won't change the REML estimates of the variance components.

Point two: Your assertion of significant contribution to variability due to water*year seems supported by the LRT. To assess the effect of year, you will have to look at the blups, using ESTIMATE statements.

Steve Denham

Re: Split-Plot design replicated in time

Ana79 — Thu, 01 Aug 2013 17:59:10 GMT

Thanks Again!