1) Baseline random intercept model

We start with a simple random intercept model for stand count $y_{ij}$, where $j$ indexes subsamples within plot $i$:

$$y_{ij}=\mu +b_i+{\epsilon }_{ij}$$

where

• $\mu$: overall mean (fixed intercept),
• $b_i\sim \mathcal{N}(0,{\sigma }_{\text{plot}}^2)$: random effect for plot,
• ${\epsilon }_{ij}\sim \mathcal{N}(0,{\sigma }_{\text{residual}}^2)$ residual variation at the subsample level.

mod1 = lmer(stand ~ 1 + (1|plot_id), data = dat_long, REML = TRUE)
VarCorr(mod1)

##  Groups   Name        Std.Dev.
##  plot_id  (Intercept) 3.6490  
##  Residual             2.9687

The between-plot variability ${\sigma }_{\text{plot}}^2$ (variation in average stand across plots) is ${3.649}^2=13.31$. This measures how different the average stand counts are between plots. A higher value indicates more variability from plot to plot. The residual variance ${\sigma }_{\text{residual}}^2={2.9687}^2=8.82$ captures how much individual subsamples differ from their plot average. It includes noise, measurement error, and row/subsample-level micro-variation.

Total variance of observation $y_{ij}$ is

$$\text{Var}(y_{ij})={\sigma }_{\text{plot}}^2+{\sigma }_{\text{residual}}^2=22.13$$

vc = as.data.frame(VarCorr(mod1))
var_plot     = vc[vc$grp == "plot_id", "vcov"]
var_resid    = attr(VarCorr(mod1), "sc")^2
var_total = var_plot + var_resid
icc = var_plot / var_total;icc

## [1] 0.6017259

2) Row side as a fixed effect

Including row_side as a fixed effect allows us to test whether there is a systematic mean difference in stand between the left and right rows across plots. However, there’s an important caveat: the labels “left” and “right” may not consistently reflect the same physical orientation in the field. For example, if the planter always moves in the same direction, left and right consistently refer to the same sides of the planter. But if planting follows a serpentine pattern (alternating directions), the same side of the planter might be labeled “left” in one pass and “right” in the next. In that case, the label reflects position relative to the viewer, not field machinery.

This limits the interpretability of row_side as a biologically or operationally meaningful factor. Still, we fit the model to assess whether there is a consistent average difference in stand between rows.

$$y_{ijk}=\mu +{\alpha }_k+b_i+{\epsilon }_{ijk}$$

where

• $\mu$: overall mean (fixed intercept),
• ${\alpha }_k$: fixed effect of side $k\in$ {left, right},
• $b_i\sim \mathcal{N}(0,{\sigma }_{\text{plot}}^2)$: random effect for plot $i$,
• ${\epsilon }_{ijk}\sim \mathcal{N}(0,{\sigma }_{\text{residual}}^2)$: residual error.

mod2 = lmer(stand ~ 1 + row_side +(1|plot_id), data = dat_long, REML = TRUE)
VarCorr(mod2)

##  Groups   Name        Std.Dev.
##  plot_id  (Intercept) 3.6492  
##  Residual             2.9677

The estimated row side effect (${\alpha }_k$) is –0.186 (not shown), suggesting slightly lower stand counts on the right row on average. However, the effect is not statistically significant. The variance components are nearly unchanged from the baseline model.

We can explore this visually by computing left–right differences in average stand per plot:

Fig. 2. Distribution of within-plot differences in stand counts between left and right row sides. Histogram displays the difference in mean stand counts per plot (left – right), based on all subsamples within each row side. The dashed vertical line at zero indicates perfect symmetry.

Another way to formally test whether mean stand differs by row side, we use a paired t-test:

t.test(stand_side_summary$left, stand_side_summary$right, paired = TRUE)

## 
## 	Paired t-test
## 
## data:  stand_side_summary$left and stand_side_summary$right
## t = 1.2551, df = 402, p-value = 0.2102
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
##  -0.1056323  0.4786678
## sample estimates:
## mean difference 
##       0.1865178

As mentioned before, while the estimated difference in stand between left and right rows is small and not statistically significant, this result should be interpreted with caution. The row_side labels may not consistently represent the same physical side of the planter or field due to alternating planting directions (serpentine patterns). As a result, any real row-to-row differences in emergence—caused by planter variability, tire compaction, or edge effects—may be masked by the inconsistent labeling.

Thus, the lack of significance does not imply row sides are truly equivalent, only that the current labeling scheme may not reliably capture physical asymmetries across plots.

3) Row side as a random effect nested within plot

If we suspect that stand counts systematically differ between the two rows within a plot—but that the meaning of “left” and “right” is inconsistent across plots—then treating row_side as a random effect nested within plot is a better approach than using it as a fixed effect.

This specification captures within-plot row-level heterogeneity without relying on the labels “left” and “right” to carry any consistent directional meaning. What matters is that each plot has two rows, and those rows may differ in stand due.

In this model, we treat row_side as a random effect specific to each plot:

$$y_{ijk}=\mu +{\alpha }_{k(i)}+b_i+{\epsilon }_{ijk}$$

• $\mu$: overall mean (fixed intercept),
• ${\alpha }_{k(i)}\sim \mathcal{N}(0,{\sigma }_{\text{row:plot}}^2)$: random effect for row side $k$ nested within plot $i$,
• $b_i\sim \mathcal{N}(0,{\sigma }_{\text{plot}}^2)$: random plot effect,
• ${\epsilon }_{ijk}\sim \mathcal{N}(0,{\sigma }_{\text{residual}}^2)$: subsample-level residual error (variation between subsamples within row).

mod3 = lmer(stand ~ 1 + (1|plot_id/row_side), data = dat_long, REML = TRUE)
VarCorr(mod3)

##  Groups           Name        Std.Dev.
##  row_side:plot_id (Intercept) 1.3785  
##  plot_id          (Intercept) 3.5433  
##  Residual                     2.7700

This model includes an additional variance component representing variation between row sides within each plot. The estimated variance ${\sigma }_{\text{plot:row}}^2={1.3785}^2=1.90$ shows that this within-plot variation is not huge but it is not zero either. On the scale of plants per acre, this corresponds to a standard deviation of approximately 10,000 plants per acre between row sides.

vc = as.data.frame(VarCorr(mod3))
var_row   = vc[vc$grp == "row_side:plot_id", "vcov"]
var_plot  = vc[vc$grp == "plot_id", "vcov"]
var_resid = attr(VarCorr(mod3), "sc")^2
var_total = var_plot + var_row + var_resid
icc_plot = var_plot / var_total;icc_plot

## [1] 0.5673863

icc_row = var_row / var_total;icc_row

## [1] 0.08587233

iic_sub= var_resid / var_total;iic_sub

## [1] 0.3467414

mod2a = lmer(stand ~ 1 + (1 | plot_id), data = dat_long, REML = FALSE)
mod3a = lmer(stand ~ 1 + (1 | plot_id/row_side), data = dat_long, REML = FALSE)
anova(mod2a, mod3a)

## Data: dat_long
## Models:
## mod2a: stand ~ 1 + (1 | plot_id)
## mod3a: stand ~ 1 + (1 | plot_id/row_side)
##       npar   AIC   BIC  logLik deviance Chisq Df Pr(>Chisq)    
## mod2a    3 13054 13072 -6524.3    13048                        
## mod3a    4 13001 13024 -6496.6    12993 55.25  1  1.061e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Theoretical standard error calculation

I use variance components from a linear mixed-effects model to derive the theoretical standard error (SE) of the mean for a given plot. Each plot contains $J=2$ row sides and $K$ subsamples per row side, for a total of $J\cdot K$ subsamples. For our case, the mean for a given plot is:

$${\overline{y}}_i=\frac{1}{J\cdot K}\sum _{j=1}^J\sum _{k=1}^K{y}_{ijk}$$

The variance of ${\overline{y}}_i$ is:

$$\mathrm{Var}({\overline{y}}_i)=\mathrm{Var}(\frac{1}{J\cdot K}\sum _{j=1}^J\sum _{k=1}^K(\mu +{\alpha }_i+{\beta }_{j(i)}+{\epsilon }_{k(ij)}))$$

Since $\mu$ is constant, and ${\alpha }_i$ is constant within plot $i$, the variance simplifies to:

$$\mathrm{Var}({\overline{y}}_i)=\mathrm{Var}({\alpha }_i+\frac{1}{J}\sum _{j=1}^J{\beta }_{j(i)}+\frac{1}{J\cdot K}\sum _{j=1}^J\sum _{k=1}^K{\epsilon }_{k(ij)})$$

Assuming independence and homoscedasticity, the variance becomes:

$$\begin{array}{rl}\mathrm{Var}({\overline{y}}_i)& ={\sigma }_{\text{plot}}^2+\frac{1}{J^2}\cdot J\cdot {\sigma }_{\text{row:plot}}^2+\frac{1}{(J\cdot K{)}^2}\cdot J\cdot K\cdot {\sigma }_{\text{residual}}^2\\ & ={\sigma }_{\text{plot}}^2+\frac{{\sigma }_{\text{row:plot}}^2}{J}+\frac{{\sigma }_{\text{residual}}^2}{J\cdot K}\end{array}$$

However, when assessing the precision of the mean estimate within a single plot (i.e., ignoring variation among plots), the between-plot variance ${\sigma }_{\text{plot}}^2$is excluded. Thus, the relevant variance is

$$\begin{array}{rl}\mathrm{Var}({\overline{y}}_i)& =\frac{{\sigma }_{\text{row:plot}}^2}{J}+\frac{{\sigma }_{\text{residual}}^2}{J\cdot K}\end{array}$$

which quantifies the sampling variance of the plot mean estimator given the nested sampling design.

Formulas are boring but they clearly show you that as $K$ increases, $\frac{{\sigma }_{\text{residual}}^2}{J\cdot K}$ decreases, improving precision. Another thing it tells you is that the $\frac{{\sigma }_{\text{row:plot}}^2}{J}$ term is constant with respect to $K$ (it does not depend on it). Thus, there is diminishing return as $K$ grows because the only term decreasing with $K$ is residual variance contribution.

Theoretical SE calculation

To find the point of diminishing returns, we calculated the theoretical SE for $K=1to10$. Even with only three subsamples per row, stats lets us extrapolate and optimize sampling—one reason it’s worth knowing the math.

calculate_se_nested = function(k, j, sigma_row, sigma_resid) {
  sqrt((sigma_row^2) / j + (sigma_resid^2) / (j * k))
}

j = 2  # row sides per plot

k_values = 1:10
se_curve_nested = tibble(
  k = k_values,
  se = round(calculate_se_nested(k_values, j = j, sigma_row = sigma_row, sigma_resid = sigma_resid),2)
)

se_curve_nested %>%
  mutate(
    se_reduction = round(lag(se) - se, 2),
    relative_reduction = round(se_reduction / lag(se), 2),
    percent_reduction = round(relative_reduction * 100, 2)
  )  %>%
  select(
    `Number of Subsamples (K)` = k,
    `Standard Error` = se,
    `SE Reduction` = se_reduction,
    `Percent Reduction (%)` = percent_reduction
  ) %>%
  kableExtra::kable(align = "c", format = "html") %>%
  kableExtra::kable_styling(full_width = FALSE, position = "center")

The theoretical standard error calculations illustrate three distinct phases in precision gains associated with increasing subsamples per row side. In the initial rapid improvement phase $(K=1\text{–}3)$, each additional subsample yields large reductions in SE (11.8–22.6%), offering a high return on sampling effort. Notably, increasing from 1 to 2 subsamples provides the largest gain (~22.6%) with minimal added effort, making 2 subsamples a practical and efficient choice in many contexts. Between $K=3\text{–}6$, diminishing returns begin: SE reductions moderate to 5–15% per additional subsample, so further increases in $K$ should be carefully balanced against labor costs. Beyond $K=6$, an asymptotic phase is reached where precision gains are not offset by additional subsampling.

I believe choosing $K=2$ does the job, on average.

Empirical validation

We use simulation to estimate the standard error of plot means across subsample sizes $K$. The goal is to identify the “elbow” in the SE curve — the point beyond which increasing $K$ yields only marginal reductions in standard error. Beyond this point, additional sampling effort does not produce proportional gains in precision.

Fig. 3. Theoretical and empirical standard errors of plot-level stand estimates as a function of subsamples per row side ($K$). Theoretical SE (red dashed line) are calculated from variance components of the nested linear mixed-effects model, while empirical SE (black solid line) are derived from simulation with repeated random subsampling. Both curves demonstrate diminishing returns in precision gains as $K$ increases, validating the theoretical model and supporting an efficient sampling intensity. $K=2$ is highlighted to show marginal gains in precision beyond that point.

Conclusion

This analysis supports collecting two 1-meter subsamples per row side when estimating stand counts in soybean trials. This sampling effort captures within-plot variability efficiently, offering a balance between precision and labor. For trials not focused on seed treatment effects, even a single subsample may suffice (Yes, you can improve your fungicide model by including stand count as a covariate. Since stand is a non-confounding effect, adjusting for it helps control variability without biasing treatment effects)