# How would you model this random effects structure?

I have a sort of weird and complicated model design, and I'd like to get your opinion on how best to model the error structure.

I have 100 sites, with each site falling into 1 of 4 different forest type categories (so 25 of each forest type). Within each site, I have 4 plots, each one with a different manipulative treatment. The outcome is number of new seedlings within each plot. I'm interested in how forest type, treatment, and their interaction affect seedling growth.

So, the most basic model I can think of is (using nlme):

lme(seedling ~ forest.type + treatment + forest.type*treatment, random=~1|site)

This sort of seems to be right, since each site probably has some random effect on the number of seedlings. But it's also doesn't seem to be totally correct, since treatment is nested within site, and the random intercept might differ across treatments. So another model I've thought of is:

lme(seedling ~ forest.type + treatment + forest.type*treatment, random=~1+treatment|site)

And, while model fits (though lmer refuses to even try), it also doesn't quite seem right, since I don't have any replication of treatment within site.

I know the geographic locations of all of the plots, so I've also tried some models that just use a spatially correlated error structure in place of a random effect, but I have no way to know if treatment affects the correlations, so I don't feel totally comfortable with this approach.

Do either of these above models seems appropriate? Or is there a different model, or a different approach that you'd suggest?

Thanks.

This is very similar to a classic split plot design.

Forest type is a fixed effect with four levels. It is a bit of a stretch to think of forest type being randomly assigned to sites. Instead, sites are a random effect nested in forest type.

Treatment is a fixed effect with three levels.

Plots are a random effect nested in sites (nested in forest type). Plots are uniquely identified in this case by site by treatment combinations.

Formally, the site by treatment interaction would be a random effect also. In this case, it is confounded with plot-to-plot variability and error variance.

The test for forest type would use the site term. The estimate for site would include the site-to-site variance $\sigma^2_{\textrm{site}}$ as well as any variability introduced in preparing each site.

The tests for treatment and for the treatment by forest type interaction would use the site by treatment term. The estimate for site by treatment term would include plot-to-plot variance $\sigma^2_{\textrm{plot}}$, site by treatment variance $\sigma^2_{\textrm{site}\times\textrm{treatment}}$, and error variance $\sigma^2_{\textrm{error}}$.

I don't have a good way yet to get correct expected mean squares and degrees of freedom for these sorts of statistical analysis in general. For this design maybe it's feasible using packages for calculating stratified error variances.

My take is that forest.type should be tested with 3 and 96 degrees of freedom, treatment should be tested with 3 and 288 degrees of freedom, and forest.type:treatment with 9 and 288 degrees of freedom.

It looks like nlme would produce the correct analysis for the fixed effects:

lme(seedling ~ forest.type + treatment + forest.type*treatment, random=~1/forest.type/site)