Erratic behavior of split-plot model: Contaminated data, script error, or bug?

Question

Problem: We have observed erratic behavior of a split-plot model when run with different subsets of our data. For two of six response-variable columns, the script produces incorrect outputs for the random factor in the model. Specifically, the random factor accounts for 0% of the total variance in the REML Variance Component Estimates table, and the table of Least Squares Means displays a column of means but not of standard errors.

We have checked the column properties and replaced the data in the problematic columns, to no avail. Two lines of evidence suggest that this behavior is not due to a script error:

* The script works properly on the other four of response variable columns 
* A reduced model—including only the random factor and the problematic data columns—produces incorrect output for one of the response variables while for the other variable the output appears to be correct

We include a brief summary of our experimental design and approach to analysis, as well as examples of the problematic data, the scripts used, and model output. We also attach a pdf showing the script and output for the reduced models described above.

Thanks for any assistance in solving this puzzle!

Experimental design: An ecological field experiment consisting of four (or three) spatially defined blocks (replicates) to which 3 factors were randomly assigned in a split-plot design. At the whole plot scale,
two two-level factors—Limpet removal (+/-) and Moisture addition (+/-), both relatively hard to change—were randomly assigned to four plots within each block. At the sub-plot scale, Blade removal (Blades +/-) was randomly assigned to two sub-plots within each whole plot (see sample layout, below).

Sample dataset

Blk#

Limp

Moist

Blades

Ja88Ap88

Mr89Nv89

1

+L

+M

+B

0

-15

1

-L

+M

+B

-5.88235

-2.94118

1

+L

-M

+B

5.88235

-11.7647

1

-L

-M

+B

-14.4118

-38.2353

2

+L

+M

+B

-2.94118

2.94118

2

-L

+M

+B

0

-20.5882

2

+L

-M

+B

0

-15

2

-L

-M

+B

0.29412

-2.94118

3

+L

+M

+B

0

.

3

-L

+M

+B

-32.3529

.

3

+L

-M

+B

-17.6471

.

3

-L

-M

+B

-15

.

4

+L

+M

+B

0

-2.64706

4

-L

+M

+B

-5.88235

-14.4118

4

+L

-M

+B

-5.58824

-43.8235

4

-L

-M

+B

-5.88235

-35.2941

1

+L

+M

-B

-5.88235

11.47059

1

-L

+M

-B

3.23529

5.29412

1

+L

-M

-B

-2.94118

-8.82353

1

-L

-M

-B

-2.94118

-11.7647

2

+L

+M

-B

-5.88235

6.17647

2

-L

+M

-B

0

-6.17647

2

+L

-M

-B

2.94118

-26.4706

2

-L

-M

-B

3.23529

-20.2941

3

+L

+M

-B

-14.7059

.

3

-L

+M

-B

14.70588

.

3

+L

-M

-B

-0.29412

.

3

-L

-M

-B

-5.88235

.

4

+L

+M

-B

-8.82353

-26.4706

4

-L

+M

-B

-26.4706

-20.5882

4

+L

-M

-B

-14.4118

-44.1177

4

-L

-M

-B

2.94118

-47.3529

Analysis. To simplify the analysis, we chose to model the four whole-plot treatment combinations as 4 levels of a single factor (Concatenate[Limp,Moist] in the script, below). We plan to use post-hoc tests to distinguish among treatment means. The random effect (error term for the whole plot) is a concatenation of the block, limpet, and moisture codes (i.e. :Name( "Concatenate[Blk#,Limp,Moist]" )[:Blk#], below).

Sample script for full model.

Fit Model(
Y( :Ja88Ap88 ),
Effects(
:Blk#,
:Name( "Concatenate[Limp,Moist]" ),
:Blades,
:Blades * :Name( "Concatenate[Limp,Moist]" )
),
Random Effects( :Name( "Concatenate[Blk#,Limp,Moist]" )[:Blk#] ),
NoBounds( 1 ),
Personality( "Standard Least Squares" ),
Method( "REML" ),
Emphasis( "Minimal Report" ),
Run(
)
),

Partial output for full model.

REML Variance component estimates For Ja88Apr88 (n = 4 blocks)

Random Effect

Var Ratio

Var Component

Std Error

95% Lower

95% Upper

Wald p-Value

Pct of Total

Concatenate[Blk#,Limp,Moist][Blk#]

-0.324822

-40.38441

27.376421

-94.04121

13.272392

0.1402

0.000

Residual

124.32778

50.756605

63.930876

338.78407

100.000

Total

124.32778

50.756605

63.930876

338.78407

100.000

Fixed effect tests For Ja88Apr88 Note denominator degrees of freedom appear to be correct.

Source

Nparm

DF

DFDen

F Ratio

Prob > F

Blk#

3

9

3.1201

0.0808

Concatenate[Limp,Moist]

3

9

0.2215

0.8791

Blades

1

12

0.3675

0.5557

Blades*Concatenate[Limp,Moist]

3

12

1.0059

0.4238

Least Squares Means Table For Ja88Apr88 (n = 4 blocks)

Level

Least Sq Mean

Std Error

[1]1, +L, +M

-2.28820

.

[1]1, +L, -M

-9.03709

.

[1]1, -L, +M

-8.62806

.

[1]1, -L, -M

8.48276

.

[2]2, +L, +M

7.78410

.

[2]2, +L, -M

-1.69162

.

[2]2, -L, +M

-3.73673

.

[2]2, -L, -M

-3.53222

.

[3]3, +L, +M

-11.31712

.

[3]3, +L, -M

-6.88603

.

[3]3, -L, +M

-11.93065

.

[3]3, -L, -M

-5.45444

.

[4]4, +L, +M

-14.25232

.

[4]4, +L, -M

-2.45880

.

[4]4, -L, +M

4.22192

.

[4]4, -L, -M

-19.56963

.

Erratic behavior of split-plot models_reduced model script and output.pdf

Sample layout and data.jpg

cbhalpern · Answer

Thank you to markbailey and statman for your thoughtful replies. We will attempt to clarify the nature of the design and to answer your data questions. We hope we will have provided enough information that you can address our initial question regarding the erratic performance of essentially identical models.

The goal of the experiment was to determine the separate and joint effects of three factors on a perennial seaweed:

* limpet herbivores (tested by retaining or removing limpets, +/-L)
* natural desiccation (tested by adding moisture, +/-M)
* physical disturbance (tested by removing blades of the seaweed, +/-B)

We are particularly interested in the seaweed’s response to and recovery from disturbance, thus the response variable is the change in seaweed abundance (cover) between sampling dates (cover[t] minus cover[t-1]). Seaweed abundance was measured in each of the 32 sampling units (sub-plots) at each sampling date.

[Statman (#4) is correct in noting that the data are not strictly continuous, but derive from point counts converted to “cover (%)” (i.e. seaweed presence/absence under a grid of 34 points, yielding a cover range of 0-100%). Because the change in cover from one time to the next can take on negative values, we cannot assume a Poisson distribution.]

Our field design is typical of ecological experiments. Treatments were randomized within spatially arrayed blocks to account for unmeasured environmental variation. Treatment assignments to sampling units remained fixed throughout the three-year experiment. To further clarify, we've attached a diagram showing

* the spatial layout of the blocks (blocks are shown as contiguous but are, in fact, separated by varying distances)
* the random assignment of L and M treatments to whole plots (factorial combinations of the two)
* the random assignment of the B treatment to sub-plots within each whole plot.

Note that the spatial scales of the treatments differ: L and M are larger (harder to manipulate) and B is smaller (easier to manipulate). Thus, randomization of the B treatment was restricted, unlike in a randomized complete block, three-factor design (as suggested by markbailey).

Based on the physical layout and randomization of treatments, it seems appropriate to conduct a split-plot analysis, testing the effects of L and M with a different error term (variation associated with whole plots) than that used for model terms that include B (variation associated with the sub-plots). However, for simplicity, we treated the combination of L and M as a single factor with 4 levels. Our plan is to explore the individual and interactive effects of L and M through post-hoc contrasts. Given low replication, more complex models that include all higher-level (2- and 3-way) interactions seem susceptible to overfitting (hopefully this answers statman #3).

The data included in our initial posting represent the two response intervals — Ja88 to Ap88 and Mr89 to Nv89 — for which we received the puzzling model output also included in our initial post. Specifically, the random factor (representing the whole-plot error term, reflecting variation among the 16 Block x L x M experimental units) accounted for 0% of the total variance in the REML Variance Component Estimates table, and the table of Least Squares Means displayed a column of means but not of standard errors. For the other four intervals, the model output appeared fine, including proper accounting of numerator and denominator dfs for the whole plots and subplots. We have tried replacing the data in the problematic columns and verified that the column properties match those of the other four columns.

Do you now have enough information to help us pinpoint a problem that occurs with some, but not all instances of running essentially identical models, replacing only the response variable (sampling interval)?

Assignment of treatments to blocks and plots.pdf

Mark_Bailey · Answer

Please read @statman's reply carefully.

I am also confused by your post:

* You have concatenated the Limpet and Moisture levels and then treated their joint changes as a random effect on the response? Wouldn't you expect these factors to produce fixed effects?
* What you call Blk# is actually the whole plot? That is, these runs reset the hard to change factor once? Then it should probably be the only random effect in the model.
* It appears that Blades is the whole plot (hard to change) factor, but then how is Blades set to both -B and +B in the same whole plot?
* Your script indicates that you are using the same factors for both the fixed and random effects. Is that intentional?
* You have lost all of whole plot = 3?

Perhaps this model will provide the information that you want as more stable estimates:

Fit Model(
	Y( :Ja88Ap88 ),
	Effects( :Limp, :Moist, :Blades, :Blades * :Limp, :Blades * :Moist ),
	Random Effects( :Blk# ),
	No Bounds( 1 ),
	Personality( "Standard Least Squares" ),
	Method( "REML" ),
	Emphasis( "Minimal Report" ),
	Run()
);

statman · Answer

Let me preface my comments by I don't understand the situation you are describing from a scientific or engineering point of view so my comments may be completely off base.  It does make sense that changing limpets and moisture would be difficult to do randomly. I don't know what your hypotheses are nor do I understand the response variable or the factors for that matter? I have some questions regarding your experiment. The difficulty I am having is understanding your design plan (based on the sample data you provided):

1. This statement makes no sense to me: "to which 3 factors were randomly assigned in a split-plot design".  Split-plot designs purposely restrict randomization, so I don't understand what you are doing.  In the sample data set, it looks like your "hard to change factors" are changing a lot!  The subplot factor doesn't appear to change nearly as much as your whole plot factors?  Do you have the run order?

2. I also don't understand your use of the block?   I presume you are using the blocks to get an estimate of the whole plot errors? Let me see if I understand what you are doing...You are running a 2^3 factorial experiment in 4 different "locations/spaces" (blocks).  Within each block you are intending to restrict randomization on 2 of 3 factors (L and M) in a factorial matrix.  This sounds more like a split-split-plot.  I added just the design in a jmp file to see if this is what your are running.

3. Why would you combine L & M?  If you do this you lose the ability to quantify the L*M interaction and any other 2nd order effect with L and M?  Instead you can estimate quadratic and cubic effects???

4. Your measurement system look to have a discrimination issue (it doesn't appear continuous).  There are only some numbers available looking at the sample data set?  Is this the real data?

Food for thought (pardon my rambling)... the reason you need to be aware and take into account restrictions on randomization is due to the potential for unanticipated changes in the noise over the course of the experiment which in turn impacts the experimental errors and subsequently the statistical tests.  The thinking is the whole plot errors will be different than the sub-plot errors and therefore a different test of significance should be performed for each.  If in planning your experiment, if you identify the noise variables and hypothesize about the affects of those noise variables, you can take that into account in the designing of your experiment (e.g., what is changing within block and what is changing between block?).  If there is very little difference between the noise changing at the whole plot level (between block) and the noise changing at the sub-pot level, then taking into account the restriction may be of no consequence.  If, however, there are significant differences in the noise (and therefore the estimated errors) at the whole plot level and the noise at the sub-plot level, then tests of significance should reflect that.

split-split-plot.jmp