Summarizing unexplained variability within groups using fitted random effects or residuals I have repeated observations (count data) conducted in different plots (say, 5 times each plot). I would like to conduct a regression analysis of the observation data at the plot level with plot-level independent variables. I'd like to consider approaches besides putting the plot-level independent variables into an observation-level model because I plan to build a structural equation model that will relate to other datasets at the plot level.
I could average the repeated observations for each plot, but each observation has individual sources of variability (e.g. weather conditions, time of day, surveyor ID, etc.) that could be systematically affecting the data, so I would like to factor them out.
Would it be possible to remove these nuisance sources of variability by either:


*

*Making an observation-level mixed-effects model with plot as a random effect grouping and running the plot-level analysis on the fitted random effects coefficients.


or


*Fitting an observation-level model without mixed effects and averaging the residuals of this model by plot groupings and running the plot-level analysis on the averaged residuals.


Are either of these an appropriate approach? Or is there a more recommended or standard approach that I'm naively not aware of?
 A: The multilevel model allows you to simultaneously model within- and between-plot variation. Consider a simplified multilevel model, with the Xs representing  within-plot variables and W representing a between-plot variable. 
The within-plot model:
$y_{ij}$ = $\beta0_j$ + $\beta1_jX_{1ij}$ + $\beta2X_{2ij}...$ + $e_{ij}$ 
Here you can include as many X variables as are necessary to explain within-plot variation in your outcome. Recognize that X variables, because they are measured within-plots, also have between-plot variation, and thus you might also consider including their plot mean in the model for the between-plot mean outcome (intercept):
$\beta0_j$  =  $\gamma_{00}$ + $\gamma_{01}W_j$ + $u_{0j}$ 
Where W could include plot means of all X variables as well as those plot-level variables you are interested in using as predictors of between-plot variance.
This model could be further extended by modeling between-plot variation in a within-plot slope. For example, the j in $\beta_{1j}$:
$\beta1_j$ (cluster slope) = $\gamma_{10}$ + $\gamma_{11}W_j$ + $u_{1j}$ 
Note that for both $\beta0_j$ and $\beta1_j$, the model assumes that these latent variables have a Normal distribution, with a mean of 0 (relative to the fixed effect parameter estimate) and an estimated variance ($\sigma^2_0$ and $\sigma^2_1$ for the intercept and slope, respectively).
Unless I am missing something, it seems that the multilevel model can accommodate both your desire to control for within-plot differences that might affect between-plot variance in the outcome and allow you to model between-plot variance in the mean level of the outcome. No need for a two-step approach where you mess with residuals from the first step.
