lowest_AIC <- mixed_model(yield ~ stock,
                    random = ~  1 | groups,
                    data = trees,
                    family = hurdle.lognormal(), n_phis = 1,
                    zi_fixed = ~ stock
)

> AIC(lowest_AIC)
[1] 2784.14
> mean(residuals(lowest_AIC)^2)
[1] 0.1474714
> summary(lowest_AIC)

Call:
mixed_model(fixed = yield ~ stock, random = ~1 | groups, data = trees,
    family = hurdle.lognormal(), zi_fixed = ~stock, n_phis = 1)

Data Descriptives:
Number of Observations: 20583
Number of Groups: 2

Model:
 family: hurdle log-normal
 link: identity

Fit statistics:
  log.Lik     AIC      BIC
 -1386.07 2784.14 2776.299

Random effects covariance matrix:
              StdDev
(Intercept) 1.396311

Fixed effects:
            Estimate Std.Err  z-value p-value
(Intercept)   0.0004  0.0061   0.0731 0.94177
stock        -0.0907  0.0062 -14.7410 < 1e-04

Zero-part coefficients:
            Estimate Std.Err  z-value p-value
(Intercept)  -2.4295  0.0257 -94.3862 < 1e-04
stock         0.1934  0.0238   8.1089 < 1e-04

log(residual std. dev.):
  Estimate Std.Err
   -0.1746  0.0051

Integration:
method: adaptive Gauss-Hermite quadrature rule
quadrature points: 11

Optimization:
method: EM
converged: TRUE

best <- mixed_model(yield ~ stock + slope + precipitation + n_trees + soil_depth + rocks  + position,
                    random = ~  1 | groups,
                    data = trees,
                    family = hurdle.lognormal(), n_phis = 1,
                    zi_fixed = ~ stock + slope + precipitation + rocks + position,
                    zi_random = ~  1 | groups

)

> AIC(best)
[1] 2839.902
> mean(residuals(best)^2)
[1] 0.1890385
> summary(best)

Call:
mixed_model(fixed = yield ~ stock + slope + precipitation + n_trees +
    soil_depth + rocks + position, random = ~1 | groups, data = trees,
    family = hurdle.lognormal(), zi_fixed = ~stock + slope +
        precipitation + rocks + position, zi_random = ~1 | groups,
    n_phis = 1)

Data Descriptives:
Number of Observations: 20528
Number of Groups: 2

Model:
 family: hurdle log-normal
 link: identity

Fit statistics:
   log.Lik      AIC      BIC
 -1387.951 2839.902 2798.083

Random effects covariance matrix:
                StdDev    Corr
(Intercept)     1.4117
zi_(Intercept)  0.1331  0.0477

Fixed effects:
                           Estimate Std.Err  z-value p-value
(Intercept)                  0.4576  0.0176  25.9341 < 1e-04
stock                       -0.0644  0.0059 -10.8639 < 1e-04
slope                       -0.0593  0.0059 -10.1307 < 1e-04
precipitation               -0.0858  0.0061 -14.1035 < 1e-04
n_trees                     -0.0630  0.0060 -10.5368 < 1e-04
soil_depth                   0.0590  0.0061   9.7462 < 1e-04
rocks2                      -0.0847  0.0155  -5.4514 < 1e-04
rocks3                      -0.0942  0.0178  -5.2897 < 1e-04
rocks4                      -0.1437  0.0184  -7.8039 < 1e-04
rocks5                      -0.1933  0.0258  -7.5044 < 1e-04
rocks6                      -0.4962  0.0574  -8.6376 < 1e-04
position-0.561514296868636  -0.2333  0.0180 -12.9644 < 1e-04
position0.381120787744946   -0.4953  0.0193 -25.6413 < 1e-04
position1.32375587235853    -0.8913  0.0234 -38.0653 < 1e-04
position2.26639095697211    -1.0993  0.0276 -39.8970 < 1e-04

Zero-part coefficients:
                           Estimate Std.Err  z-value  p-value
(Intercept)                 -3.5587  0.1180 -30.1541  < 1e-04
stock                        0.1724  0.0255   6.7729  < 1e-04
slope                        0.2517  0.0272   9.2508  < 1e-04
precipitation               -0.1457  0.0296  -4.9197  < 1e-04
rocks2                       0.3587  0.0719   4.9900  < 1e-04
rocks3                       0.3497  0.0823   4.2506  < 1e-04
rocks4                       0.3951  0.0845   4.6742  < 1e-04
rocks5                       0.5387  0.1110   4.8507  < 1e-04
rocks6                       0.5434  0.2539   2.1404 0.032319
position-0.561514296868636   0.3085  0.1215   2.5384 0.011136
position0.381120787744946    0.8707  0.1219   7.1406  < 1e-04
position1.32375587235853     1.5953  0.1251  12.7515  < 1e-04
position2.26639095697211     2.0707  0.1270  16.3061  < 1e-04

log(residual std. dev.):
  Estimate Std.Err
   -0.2856  0.0052

Integration:
method: adaptive Gauss-Hermite quadrature rule
quadrature points: 11

Optimization:
method: hybrid EM and quasi-Newton
converged: TRUE