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?