I am currently working on some analyses where I examine the effects of a variety of landscape-scale variables on bee abundance.

Bees were collected in pan traps in a nested sampling design as follows:

- 97 landscapes of interest in an agricultural environment.

- Within each landscape, bees were collected in 3 crop fields.

- For each field, 6 traps were deployed along a transect at the edge of the field, and 6 traps were deployed along a transect actually in the field.

- Each transect in each landscape was sampled twice over the course of the summer.

This being an ecological field study, things did not go perfectly as planned... Sometimes some or all of the pan traps in a transect in one of the rounds of sampling were lost/destroyed.

I want to report all of my data at a landscape scale (i.e. pooling my bee abundance data across sampling rounds, transects, and fields within each landscape, so that I can easily assess the effects of my landscape predictor variables). I know that I can partially deal with the problem of differential trap recovery by using bee density (number of bees per successfully recovered trap) as a response variable, rather than total bee counts. However, this does not fully account for the effects of field, transect type (in-field or field border), or sampling round (1 or 2). For example, if bees naturally vary in abundance based on sampling round, a landscape that had more traps recovered in Round 1 than in Round 2 will skew my final abundance estimates.

Can anyone suggest a way that I can account for differential sampling effort across these hidden, nested variables in a linear model of the form, "landscape-scale bee abundance ~ landscape variable 1 + landcape variable 2 + ... + landscape variable n"?

Solutions that would work in R would be ideal.

Thank you so much!