unbalanced, incomplete split-plot

[palolix]

Hi!

I have a Split-plot with 3 Factors; nematode, treatment and cultivar, plus block.  The design is unbalanced because not all levels of the factor "nematode" have the same number of replicates.  Additionaly, is also an incomplete block because one of the 4 blocks only contain 2 levels of "nematode", and only one level of "treatment".  There is still another limitation; not all the levels of "cultivar" correspond to each of the level of "nematode".  For example, the cultivar Nord was only planted where nematode A and B are, but not where nematode C is.  Cultivar Sang, where only planted where nematode A and C are, but not where B is, and cultivar Beretta was planted where nematode B and C are, but not where A is.  So taken these limitations into account, I would first like to know which design do I have exactly?

I analysed the data using the by statement (by nematode and by cultivar), but the problem is that analysing the data separatly (first by nematode and then by cultivar), I loose degrees of freedom, so better to analyse it toghether in a properly way.

I would greatly appreciate some help!

Caroline

[SteveDenham]

This is a classic case where a mixed model approach is absolutely required, but that you may have to fit a cell means model to get what you want, due to the potential inestimability of main effect means.  The difficult situation arises from the block with only one treatment.  However, I again refer to Walt Stroup's Generalized Linear Mixed Models, and in particular to section 2.4.3 "Fixed or Random?: Tough Calls", where a balanced incomplete block design is described.  This should generalize easily to your situation.

Steve Denham

[palolix]

Hi Steve!

thank you again for your help!  I will certainly read that book.  However, what do you mean with fit a cell means model??  Are there special sas statements for this kind of design?

Thank you!

Caroline

[SteveDenham · Accepted solution]

Another recommended book is Milliken and Johnson's Analysis of Messy Data, both volumes 1 and 3.  In their approach, you essentially set up a one-way analysis of variance, and then compute tests and comparisons using CONTRAST/ESTIMATE/LSMESTIMATE statements.  For instance, you might fit:

proc glimmix;

class nematode cultivar treatment block;

model response=nematode*cultivar*treatment/noint solution;

random block;

lsmeans nematode*cultivar*treatment/e;

run;

From the /e option on the lsmeans statement, you can get the coefficients to construct the tests for main effects and first order interactions.  This is about the only way I know to address the situation you find yourself in.  It is also very appropriate for "wedge" designs.

Steve Denham

[palolix]

Thanks a lot Steve! you are greatly helping me.  I was looking at the Contents of Analysis of Messy Data, and it looks very helpfull for my case, so I will try to get that book.  I also ordered that book from Stroup.  Now my question is; are there specific statements in Proc mixed or in glimmix to construct a model which is unbalanced, has an incomplete block, and also missing treatment combinations, like in my situation?

Thank you very much for that model! what about the noint solution? Should I always omit the intercept, or just when dealing with unbalanced design for example?

Caroline

[SteveDenham]

Noint fits cell means, and solution will return those means as the result.  If you delete noint, the intercept estimate would be the grand mean, and the estimates for each of the combinations of nematode*cultivar*treatment would be deviations from that mean.

As far as specific statements, well, no, you don't really tell MIXED or GLIMMIX that you have some special situation, except by doing things like fitting a cell means model and creating tests using other statements.  One of the great things about a mixed model approach is that you don't need to be concerned about getting an answer when you have unbalance or an incomplete block, since the conditional expectation given the random effects is what you are fitting.  It's the missing treatment combination part that gets tricky.  The highest order interaction will always be estimable, but lower order and main effects may not--hence the use of a cell means model.

Steve Denham

[palolix]

Thank so much Steve!

So by fitting a cell means you mean analysing the data as it were a One-way ANOVA (when actually being a split-plot) in order to solve the problem of missing treatment combinations?

And what about the random error of the split- plot model; blk*trt??  Can I also add it near blk in the random statement?

Thank you!

Caroline

[SteveDenham]

You analyze it as a one-way anova, but specify the necessary random effects, so that tests generated by contrast/lsmestimate statements are correct.  And, yes, you should add blk*trt into the random statement.  I would do it as follows:

random intercept treatment/subject=block;

I would not be surprised if the block by treatment variance component was zero, though.  If that turns out to be the case, then there are two choices:

1. Use the ddfm=kr option (=kr2 in GLIMMIX) to get the shrinkage estimates with proper standard errors.

2. Use method=laplace or method=quad (only available in GLIMMIX).

Note that these two are mutually exclusive.

Steve Denham

[palolix]

Thank you Steve!

Since I have so many random errors (because of my original split-plot model), I do not know if I can still use your new suggestion (random intercept treatment/subject=blk).  So if not, it should look like this?

proc glimmix;

class blk nematode treatment cultivar;

model y= nematode*treatment*cultivar/noint solution;

random blk blk*nematode  blk*nematode*treatment  blk*nematode*cultivar blk*treatment*cultivar;

lsmeans nematode*treatment*cultivar/e;

run;

But what about the main effects and first order interactions from my original split-plot model: nematode  treatment nematode*treatment cultivar nematode*cultivar treatment*cultivar nematode*treatment*cultivar?

You meant to first run the lsmeans statement with the /e option to get the coefficients in order to construct the tests for these main effects and first order interactions?  So doing this I could latter include these in the model?

Thank you very much!

Caroline

[SteveDenham]

I stand corrected, and apologize.  You can include them in the model--you just won't be able to get estimable LSmeans for all of them.  Consequently, we generally construct the tests using contrast/lsmestimate statements, but you don't have to.

As far as the random statement, you will probably have more luck with convergence if you use the subject= option.  This would make the statement look like;

random intercept nematode nematode*treatment nematode*cultivar treatment*cultivar/subject=blk;

The resulting code should then give you the tests and lsmeans that are estimable:

proc glimmix;

class blk nematode treatment cultivar;

model y= nematode|treatment|cultivar/ solution;

random intercept nematode nematode*treatment nematode*cultivar treatment*cultivar/subject=blk;

lsmeans nematode*treatment*cultivar/e;

run;

Steve Denham

Message was edited by: Steve Denham

[palolix]

Thanks a lot Steve for your help!

I run your model, it worked and that is what I got in Log, for the first dependent variable analysis:

NOTE: The model does not contain an intercept. Columns of X are scaled only and not centered.

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

NOTE: Estimated G matrix is not positive definite.

Since I got significant p-values for nematode, nematode*treatment and nematode*cultivar, I run the model for second time, replacing the lsmeans statement with these significant effects and interactions, but unfortunately I only got LSM for nematode*cultivar.  For nematode and nematode*treatment  I got Non-est, no Std error and no p-values.

Caroline

[SteveDenham]

This is what I meant by non-estimable, and why we fit a cell means model.  To get the LSMs for nematode by treatment, you will have to construct them from the available estimates of nematode by treatment by cultivar.  This is why the LSMESTIMATE statement was devised, in my opinion.  Go back to the lsmeans statement with the three-way interaction, add the /e option, and see which combinations need to be made to come up with the LSmeans of interest.

Steve Denham

[palolix]

Hi Steve!

Thanks a lot again!

I run the model again with the lsmeans statement with the three-way interaction, and that is what I got;

In the LSM Table for the nematode*treatment*cultivar interaction, I got all estimates for all combinations, and all p-values significant for all combinations.  In the Solutions for Fixed Effects Table; only for the combinations nematode1*treatment1 and nematode2*treatment1 I got estimates and p-values, for the others combinations corresponding to this interaction, I got estimate 0 and no p-values.

Caroline

[SteveDenham]

The solutions are going to come up zero for the same reason that the LSmeans are non-estimable.  I think I could derive the proper statements if you could share the three-way LSmeans output, including the /e part.  I assume that you have 3 levels of cultivar, 3 levels of nematode, and a reasonable number of levels of treatment.

Steve Denham

[palolix]

Hi Steve, I really appreciate your help!!!!

I have 3 levels of cultivar, 3 of nematode (but not all levels of nematode correspond to each level of cultivar, that is why I have missing combinations), and only 2 levels of treatment.

So, here is the 3-way LSMEANS Table, and next the Solution Table (I am sorry but I do not know a better option to paste them better)

Nordster	1	1	0.6820	0.09826	2	6.94	0.0201
Sang	1	1	0.6051	0.09826	2	6.16	0.0254
Nordster	1	2	0.7213	0.1085	2	6.65	0.0219
Sang	1	2	0.8315	0.1085	2	7.67	0.0166
Beretta	2	1	1.5122	0.1048	2	14.43	0.0048
Nordster	2	1	1.5219	0.1048	2	14.52	0.0047
Beretta	2	2	1.3954	0.1094	2	12.75	0.0061
Nordster	2	2	1.6419	0.1094	2	15.01	0.0044
Beretta	3	1	1.0348	0.09826	2	10.53	0.0089
Sang	3	1	0.9377	0.09826	2	9.54	0.0108
Beretta	3	2	1.3250	0.1085	2	12.22	0.0066
Sang	3	2	1.2005	0.1085	2	11.07	0.0081

nem		1		0.8315	0.1085	5	7.67	0.0006
nem		2		1.2710	0.1305	5	9.74	0.0002
nem		3		1.2005	0.1085	5	11.07	0.0001
trt			1	-0.2628	0.1083	4	-2.43	0.0722
trt			2	0	.	.	.	.
nem*trt		1	1	0.03636	0.08573	6	0.42	0.6862
nem*trt		1	2	0	.	.	.	.
nem*trt		2	1	0.4069	0.08880	6	4.58	0.0038
nem*trt		2	2	0	.	.	.	.
nem*trt		3	1	0	.	.	.	.
nem*trt		3	2	0	.	.	.	.
cultivar	Beretta			0.1245	0.1178	4	1.06	0.3501
cultivar	Nordster			0.3710	0.1376	4	2.70	0.0544
cultivar	Sang			0	.	.	.	.
nem*cultivar	Nordster	1		-0.4811	0.1231	8	-3.91	0.0045
nem*cultivar	Sang	1		0	.	.	.	.
nem*cultivar	Beretta	2		0	.	.	.	.
nem*cultivar	Nordster	2		0	.	.	.	.
nem*cultivar	Beretta	3		0	.	.	.	.
nem*cultivar	Sang	3		0	.	.	.	.
trt*cultivar	Beretta		1	-0.02736	0.1506	4	-0.18	0.8646
trt*cultivar	Nordster		1	-0.2641	0.1771	4	-1.49	0.2101
trt*cultivar	Sang		1	0	.	.	.	.
trt*cultivar	Beretta		2	0	.	.	.	.
trt*cultivar	Nordster		2	0	.	.	.	.
trt*cultivar	Sang		2	0	.	.	.	.
nem*trt*cultivar	Nordster	1	1	0.4513	0.1515	2	2.98	0.0967
nem*trt*cultivar	Sang	1	1	0	.	.	.	.
nem*trt*cultivar	Nordster	1	2	0	.	.	.	.
nem*trt*cultivar	Sang	1	2	0	.	.	.	.
nem*trt*cultivar	Beretta	2	1	0	.	.	.	.
nem*trt*cultivar	Nordster	2	1	0	.	.	.	.
nem*trt*cultivar	Beretta	2	2	0	.	.	.	.
nem*trt*cultivar	Nordster	2	2	0	.	.	.	.
nem*trt*cultivar	Beretta	3	1	0	.	.	.	.
nem*trt*cultivar	Sang	3	1	0	.	.	.	.
nem*trt*cultivar	Beretta	3	2	0	.	.	.	.
nem*trt*cultivar	Sang	3	2	0	.	.	.	.