rm(list = ls())N-Response OSR, N-Versuch Hohenschulen
1 Prices and other settings
For the economic evaluation of the N rate experiment prices for nitrogen and oilseed rate have to be defined. Additionally some constants are needed to calculate the N demand for OSR according to the german fertiliser ordinance.
#### OSR price in Euro/t
RPr <- 500
##### N-price in Euro/kg N
#NPr <- 0.9
NPr <- 1.2
#ParzF <- 0.9 # Euro/kg N
Sollwert <- 200 # [kg N/ha] alter Sollwert
# reference yield level of DuV [t/ha]
Yref <- 4
# maximum yield level plus for N supply
max_inc <- 40
# Sensitivity of N supply level to av. yield [kg N/0.5 t/ha]
dec_5 <- 15
inc_5 <- 10
Nmin_FJ <- 30 # [kg N/ha] Nmin in Fruehjahr, hier hypothetisch angenommen2 Data input
The experimental data were imported from an Excel file. The data are stored at: “V:/01/Erträge 01.xlsx”
#dat <- read.xlsx("Erträge 01.xlsx", sheetName = "Raps als Prüfkultur")
#dat <- read_excel("V:/01/Erträge 01.xlsx", sheet = "Raps als Prüfkultur")
dat <- read_excel("Erträge 01.xlsx", sheet = "Raps als Prüfkultur")
# some renaming and calculation of additional variables
dat$yield <- dat$`Ertrag (t/ha bei 91% TS)`
dat$`Ertrag (t/ha bei 91% TS)` <- NULL
dat$Oel_91 <- as.numeric(dat$`Öl NPZ (% bei 100% TS)`)*0.91
dat$Oel_100 <- as.numeric(dat$`Öl NPZ (% bei 100% TS)`)
dat$Npr <- as.numeric(dat$`% N im Samen (6.25)`)
dat$HerbstNDgg <- ifelse(dat$NH>0,"NH40" ,"NH0")
dat$NHerbst <- dat$HerbstNDgg
dat$Nentz <- dat$yield*0.91*dat$Npr*10
dat$TKM <- as.numeric(dat$`TKM (g) bei 91% TS`)
dat$`% N im Samen (6.25)` <- NULL
dat$`Öl NPZ (% bei 100% TS)` <- NULL
dat$`Ertrag (t/ha bei 91% TS)` <- NULL
dat$`TKM (g) bei 91% TS` <- NULL
dat$Nges <- dat$N1+dat$N2+dat$N3
dat$Nges2 <- dat$Nges^2
dat$OilYield <- dat$yield*0.91*dat$Oel_100/100
# the price is dependend on oil concentraation
# for every 1% above 40% oil concentration, the price is increased by 1.5% of the base price
dat$price <- RPr+(dat$Oel_91-40)* 1.5/100*RPr
#dat$price <- RPr
dat$Return <- dat$yield*dat$price
dat$NetReturn <- dat$Return-(dat$NH+dat$Nges)*NPr
dat$F_Nmin0_90 <- Nmin_FJ
dat$H_Nha <- 40
dat$key2 <- paste(dat$Jahr, dat$HerbstNDgg,sep = "_")
dat$key <- paste(dat$Jahr, dat$HerbstNDgg, dat$Schlag,sep = "_")
#dat <- dat %>% filter (N1 != 0, N2 != 0)
# save Jahr as numeric variable
dat$nJahr <- dat$Jahr
# convert Jahr and Schlag to factors
dat$Jahr <- as.factor(dat$Jahr)
dat$Schlag <- as.factor(dat$Schlag)3 Data summary
At the Hohenschulen experimental farm, nitrogen (N) rate trials for winter oilseed rape were conducted in the years 2016 to 2024.
In autumn, either no nitrogen or a rate of 40 kg N/ha was applied. In spring, different nitrogen amounts were applied at three dates, with the usual first and second applications following a gradient of 0, 40, 80, or 120 kg N/ha. Additionally, a late nitrogen application of 40 kg N/ha was tested against a control with 0 kg N/ha.
All factor levels of the split applications were fully combined, resulting in a total of 64 nitrogen treatment variants.
This evaluation focuses solely on the response of the target parameters—yield, nitrogen fertilizer cost free econimc net return, oil content, and nitrogen concentration in the seed—as a function of total spring nitrogen application.
No soil mineral nitrogen (Nmin) analysis for the 0–90 cm layer was conducted in spring. Therefore, a default Nmin value of 30 kg N/ha was assumed for this dataset, representing the typical value at the beginning of the vegetation period for winter oilseed rape with relatively low variation.
4 Statistical analysis of the N response data
4.1 Simple summary statistics
The summary statistics of the N-trial data set are shown in Table 1. The data set contains 64 individual N response series, each with a different combination of N application rates in autumn and spring. The target parameters yield, oil concentration at 91% moisture, oil yield, price, and economic net return are summarized.
| Ertrag | NKFL | Oel | Oelertrag | Preis | |
|---|---|---|---|---|---|
| Mean | 3.46 | 1675 | 45.3 | 1.57 | 539.8 |
| Median | 3.32 | 1566 | 45.4 | 1.48 | 540.6 |
| Max | 6.37 | 3197 | 50.2 | 2.95 | 576.7 |
| Min | 0.76 | 357 | 39.9 | 0.34 | 498.9 |
4.2 Simple quadratic model
The experimental data for yield, oil content at 91% moisture, nitrogen concentration of seeds, N offtake, econonomic net return were fitted as quadratic functions of total N spring N supply separately for plots with and without autumn N dressing. Effects of the factor “Schlag” (field) thereby are not included in the model. Thereby the function Quadpars (source file Nopt_OSR_functions.R) is used, which fits a quadratic function to the data and returns the coefficients of the function as well as the R-squared value. The function is applied to each individual N response series separately.
# quadratic fit for seed yield
df.par.q.yield <- QuadPars(dat, IndVar = "Nges", DepVar = "yield", key = "key")
df.par.q.oil <- QuadPars(dat, IndVar = "Nges", DepVar = "Oel_91", key = "key")
df.par.q.Npr <- QuadPars(dat, IndVar = "Nges", DepVar = "Npr", key = "key")
df.par.q.OilYield <- QuadPars(dat, IndVar = "Nges", DepVar = "OilYield", key = "key")
dat$Nentz <- dat$yield*0.91*dat$Npr*10
df.par.q.Nentz <- QuadPars(dat, IndVar = "Nges", DepVar = "Nentz", key = "key")
#df.par.q.Nentz <- QuadPars(dat[!is.na(dat$Nentz),], IndVar = "Nges", DepVar = "Nentz", key = "key")
df.par.q.NetReturn <- QuadPars(dat, IndVar = "Nges", DepVar = "NetReturn", key = "key")
## Zusammenfassung der Koeffizienten in einem data frame
Oilcols <- c("key", "Y0_Oel_91", "slope_Oel_91", "quad_Oel_91", "r.squared_Oel_91")
temp <- df.par.q.oil[,Oilcols]
df.par.q <- merge(x = df.par.q.yield, y = temp, by="key")
Nprcols <- c("key", "Y0_Npr", "slope_Npr", "quad_Npr", "r.squared_Npr")
temp <- df.par.q.Npr[,Nprcols]
df.par.q <- merge(x = df.par.q, y = temp, by="key")
# add columns with codes
#df.par.q <- merge(df.par.q.yield, df.par.q.oil, by = "key")
#df.par.q <- merge(df.par.q, df.par.q.Npr, by = "key")
# Datensaetze ohne Anzeichen fuer abnehmenden Ertragszuwachs werden rausgeschmissen ..
df.par.q <- df.par.q[df.par.q$quad_yield <0,]
df.par.q.oil <- df.par.q.oil[df.par.q.oil$key %in% df.par.q$key,]
df.par.q.Npr <- df.par.q.Npr[df.par.q.Npr$key %in% df.par.q$key,]
df.par.q.OilYield <- df.par.q.OilYield[df.par.q.OilYield$key %in% df.par.q$key,]
df.par.q.NetReturn <- df.par.q.NetReturn[df.par.q.NetReturn$key %in% df.par.q$key,]
# calculation of predicted yields for evaluating the residuals
dat$yield_pred_q <- NA
for (i in 1:length(df.par.q$key)){
key <- df.par.q[i, "key"]
Nges <- dat[dat$key==key, "Nges"]
dat[dat$key==key, "yield_pred_q"] <- YRC (N = Nges, Int = df.par.q[i,"Y0_yield"],
Lin = df.par.q[i,"slope_yield"],
Quad = df.par.q[i, "quad_yield"])
}
dat$resi_q <- dat$yield-dat$yield_pred_q
#rm(df.par.q.oil)
rm(df.par.q.yield)
#df.par.q <- merge(df.par.q, MeanVeg, by = "key", all.x = T)
#rm(MeanVeg)
#df.par.q <- merge(df.par.q, df_regoilYCoefw, by = "key")
#rm(df_regoilYCoefw)
#df.par.q <- merge(df.par.q, df.par.q.NetReturn, by = "key")
#rm(df.par.q.NetReturn)
#df.par.q <- merge(x = df.par.q, y = df.key, by = "key", all.x = T )
df.par.q$model <- "Q"
df.par.q <- as.data.frame(df.par.q)
write.csv(x = df.par.q, file = "dfqpar_Coefficients.csv")4.3 Plot with fitted quadratic yield models
For the prediction of the yield as a function of N supply, a grid of N rates is created. The yield is predicted using the fitted quadratic model coefficients for each individual N response series.
Figure 1 shows the predicted yield as a function of N supply for each year and whether autumn N was applied or not. In general the yield was boosted by autumn N application.
4.4 Mixed effect models for yield
Another approach to analyze the data is to use a mixed effect model where the N rate (and its squared value) are treated as fixed effects and the factors year and block (Schlag) are treated as random effects. For the random effects only their variance will be estimated, which increases the power of the model. Two appproaches were tested, one with a simple random intercept model and one with a random slope model. The latter allows the slope of the relationship between N rate and yield to vary by year.
4.4.1 Simple mixed effect model
The function lme from the package nlme is used to fit a linear mixed effect model to the data. The model includes the fixed effects Nges, NHerbst and Nges2, as well as the interaction of Nges and NHerbst which turned out to be not significant in the ANOVA test. The random effects are Schlag (field) and Jahr (year). The weights argument is used to specify that the variance of the residuals should be allowed to vary by year, which is important for the model to be able to handle heteroscedasticity in the data (Eq. 1).
The model is specified as:
\begin{align*} Y_{ijk} &= \beta_0 + \beta_1 \cdot \text{Nges}_{ijk} + \beta_2 \cdot \text{NHerbst}_{ijk} + \beta_3 \cdot \text{Nges2}_{ijk} \\ &\quad + \beta_4 \cdot (\text{Nges}_{ijk} \cdot \text{NHerbst}_{ijk}) + \beta_5 \cdot (\text{Nges2}_{ijk} \cdot \text{NHerbst}_{ijk}) \\ &\quad + u_j + v_{k(j)} + \varepsilon_{ijk} \end{align*} \tag{1}
Where:
- Y_{ijk}: observed yield for observation i, plot k, year j
- \beta_0, \dots, \beta_5: fixed-effect coefficients
- u_j \sim \mathcal{N}(0, \sigma^2_u): random intercept for year j
- v_{k(j)} \sim \mathcal{N}(0, \sigma^2_v): random intercept for plot k within year j
- \varepsilon_{ijk} \sim \mathcal{N}(0, \sigma^2_j): residual error, with year-specific heteroscedastic variance
# Fit a mixed effect model with fixed effects Nges, NHerbst and Nges2 and random effects Schlag and Jahr
Yield.M2 <- lme(fixed = yield ~ Nges*NHerbst+Nges2*NHerbst,
random = ~1|Jahr/Schlag,
data=dat,
weights = varIdent(form = ~1|Jahr))4.4.1.1 Anova
The anova of the mixed effect model M2 shows that the fixed effects Nges, NHerbst and Nges2 are significant but not their interaction (Table 2).
| numDF | denDF | F-value | p-value | |
|---|---|---|---|---|
| (Intercept) | 1 | 1126 | 140.912 | 0.000 |
| Nges | 1 | 1126 | 2445.794 | 0.000 |
| NHerbst | 1 | 1126 | 147.545 | 0.000 |
| Nges2 | 1 | 1126 | 391.033 | 0.000 |
| Nges:NHerbst | 1 | 1126 | 3.346 | 0.068 |
| NHerbst:Nges2 | 1 | 1126 | 0.965 | 0.326 |
4.4.1.2 Parameter estimates and econmic optimum N rate and yield
The economic optimum is in a first simple version calculated from the yield effect alone without considering the effects of the oil concentration in the seed on the product price.
The estimated effects and parameters are shown in Table 3. There is only a very little difference in the estimated economic N-Optimum (Nopt_q). Yield at economic optimum, however seems to be more affected.
| WithAutumnN | NoAutumnN | |
|---|---|---|
| Nmax | 240.227 | 237.713 |
| Ymax | 4.070 | 3.871 |
| Nopt_q | 205.270 | 206.063 |
| Yopt_q | 4.028 | 3.833 |
| (Intercept) | 2.089 | 1.728 |
| Nges | 0.016 | 0.018 |
| Nges2 | 0.000 | 0.000 |
4.4.2 Mixed effect model with random slope
The model M3 (Eq. 2) is a random slope model that allows the slope of the relationship between Nges and yield to vary by year. This means that the effect of Nges on yield is not constant across years, but varies depending on the year (Jahr). The pdIdent function is used to specify the random effects structure. The model thereby is also fitted using the lme function from the nlme package. Again the weights argument is used to specify that the variance of the residuals should be allowed to vary by year (weights = varIdent(form = ~1|Jahr)), which is important for the model to be able to handle heteroscedasticity in the data. This model does not include the interaction within the random effects structure.
The random statement is Jahr = pdIdent(~1) and adds a random intercept for each year (Jahr).
pdIdent(~1) specifies a diagonal (independent) covariance structure: each level of Jahr gets its own intercept with a variance component, but no covariance with other random effects.
Jahr = pdIdent(~Schlag - 1)
Adds a random effect of Schlag (field or block) within each year.
~Schlag - 1 removes the intercept and includes only the random slopes (i.e., deviations for each Schlag).
Interpreted as: “allow for random deviations between fields (Schlag) nested within years (Jahr)”.
- Jahr = pdIdent(~Nges - 1)
Adds a random effect of nitrogen treatment level (Nges) within each Jahr.
Again, no intercept (-1), only the random effects associated with Nges.
Let Y_{ijkl} be the yield for:
year i,
block/field j,
nitrogen level k,
observation l within that structure.
The model is specified as:
\begin{align*} Y_{ijk} &= \beta_0 + \beta_1 \cdot \text{Nges}_{ijk} + \beta_2 \cdot \text{NHerbst}_{ijk} + \beta_3 \cdot \text{Nges2}_{ijk} \\ &\quad + \beta_4 \cdot (\text{Nges}_{ijk} \cdot \text{NHerbst}_{ijk}) + \beta_5 \cdot (\text{Nges2}_{ijk} \cdot \text{NHerbst}_{ijk}) \\ &\quad + u_j + v_k + w_{ijk}^{(1)} + w_{ijk}^{(2)} + \varepsilon_{ijk} \end{align*} \tag{2} Where:
Where:
- Y_{ijk}: observed yield for observation i, in year j, plot k
- \beta_0, \dots, \beta_5: fixed-effect coefficients
- u_j \sim \mathcal{N}(0, \sigma^2_u): random intercept for year j
- v_k \sim \mathcal{N}(0, \sigma^2_v): random effect of plot k (from
pdIdent(~Schlag - 1))
- w_{ijk}^{(1)} \sim \mathcal{N}(0, \sigma^2_{Nges}): random effect associated with nitrogen level Nges
- w_{ijk}^{(2)} \sim \mathcal{N}(0, \sigma^2_{Nges2}): random effect associated with quadratic nitrogen term Nges2
- \varepsilon_{ijk} \sim \mathcal{N}(0, \sigma^2_j): residual error with year-specific heteroscedastic variance
# the first fit is without the random slope for Nges2
Yield.M3 <- lme(fixed = yield ~ Nges*NHerbst+Nges2*NHerbst,
random = list(Jahr = pdIdent(~1),
Schlag = pdIdent(~Schlag - 1),
Nges = pdIdent(~Nges - 1)#,
#Jahr = pdIdent(~Nges2 -1)#,
),
data=dat,
weights = varIdent(form = ~1|Jahr))#,
# start= start_vals)
# then the model is updated to include the random slope for Nges2
Yield.M3.v2 <- update(Yield.M3,
random = list(Jahr = pdIdent(~1),
Jahr = pdIdent(~Schlag - 1),
Nges = pdIdent(~Nges - 1),
Nges2 = pdIdent(~Nges2 -1)),
control = lmeControl(maxIter = 100, msMaxIter = 100, tolerance = 1e-6))
# save the fitted values of the models
fit.M3 <- as.data.frame(Yield.M3$fitted)
fit.M3.v2 <- as.data.frame(Yield.M3.v2$fitted)
# copy the updated model to Yield.M3
Yield.M3 <- Yield.M3.v2The estimated random effects for the model M3 can be extracted using the ranef function. For the random effect of Nges2 extremely small effects were estimated.
# | eval: false
#| warning: false
#| include: false
#summary(Yield.M3)
#VarCorr(Yield.M3)
head(ranef(Yield.M3)[["Nges2"]]) Nges2
2016/2016/0/0 0.000000e+00
2016/2016/40/1600 3.455160e-12
2016/2016/80/6400 -3.819133e-12
2016/2016/120/14400 -1.071143e-11
2016/2016/160/25600 -7.744293e-12
2016/2016/200/40000 -1.290150e-114.4.2.1 Anova
The anova of this model shows that the fixed effects Nges, NHerbst and Nges2 are significant and in difference to model Yield.M2 also their interaction was significant (Table 4) .
| DF | Sum Sq | Mean Sq | F value | p-value |
|---|---|---|---|---|
| (Intercept) | 1 | 1074 | 101.262 | 0.000 |
| Nges | 1 | 61 | 890.575 | 0.000 |
| NHerbst | 1 | 1074 | 184.350 | 0.000 |
| Nges2 | 1 | 61 | 101.140 | 0.000 |
| Nges:NHerbst | 1 | 1074 | 4.815 | 0.028 |
| NHerbst:Nges2 | 1 | 1074 | 1.489 | 0.223 |
4.4.2.2 Estimated parameters and N optima
Again, the estimates between optimum Nrate in spring for plots with or without N application in autumn were rather small (Table 5).
| WithAutumnN | NoAutumnN | |
|---|---|---|
| Nmax | 239.894 | 236.976 |
| Ymax | 4.103 | 3.900 |
| Nopt_q | 206.681 | 207.132 |
| Yopt_q | 4.063 | 3.864 |
| (Intercept) | 2.024 | 1.642 |
| Nges | 0.017 | 0.019 |
| Nges2 | 0.000 | 0.000 |
4.4.3 Comparison of Mixed Effect Models M2 and M3
The two mixed effect models Yield.M2 and Yield.M3 can be compared using anova. The model M3 allows the slope of the relationship between Nges and yield to vary by year, while M2 does not. The comparison shows that the model M3 is significantly better than model M2, indicating that the random slope model is more appropriate for the data (Table 6).
| call | Model | df | AIC | BIC | logLik | Test | L.Ratio | p-value | |
|---|---|---|---|---|---|---|---|---|---|
| Yield.M2 | lme.formula(fixed = yield ~ Nges * NHerbst + Nges2 * NHerbst, data = dat, random = ~1 | Jahr/Schlag, weights = varIdent(form = ~1 | Jahr)) | 1 | 17 | 995.642 | 1081.346 | -480.821 | NA | NA | |
| Yield.M3 | lme.formula(fixed = yield ~ Nges * NHerbst + Nges2 * NHerbst, data = dat, random = list(Jahr = pdIdent(~1), Jahr = pdIdent(~Schlag - 1), Nges = pdIdent(~Nges - 1), Nges2 = pdIdent(~Nges2 - 1)), weights = varIdent(form = ~1 | Jahr), control = lmeControl(maxIter = 100, msMaxIter = 100, tolerance = 1e-06)) | 2 | 19 | 880.730 | 976.517 | -421.365 | 1 vs 2 | 118.912 | 0 |
4.4.4 Plot of measured vs. predicted for Yield.M3
For plotting of the predicted yield from the mixed effect model M3, a grid of N rates is created. The predictions are made for each year and field (Schlag) combination. The predictions are stored in the data frame grid, which is then used for plotting the predicted yield as a function of N supply.
grid <- expand.grid(Jahr = unique(dat$Jahr), Nges = seq(0, 300, by = 10), NHerbst = unique(dat$NHerbst), Schlag= unique(dat$Schlag))
grid$Nges2 <- grid$Nges^2
grid$Schlag <- as.factor(grid$Schlag)
# the prediction includes the random effects for Schlag and Jahr (level = 1)
grid$Yield.M3.pred <- predict(Yield.M3, newdata = grid, level = 1)The predicted vs. measured yield from the mixed effect model M3 is shown in Figure 2. The predictions are shown for each year and on average of the blocks (Schlag). The points show the measured yield. In comparison to the single fit of the quadratic model it becomes evident, that the random interaction Nges2*Jahr was not predicted accurately (compare Figure 1).
4.4.5 Mixed effects model with lmer
Another approach of analysis is to use the function lmer from the lme4 package to fit a linear mixed effect model to the data. The model includes the fixed effects Nges, NHerbst and Nges2, as well as the interaction of Nges and NHerbst and Nges2 with NHerbst. The random effects are Schlag (field) and Jahr (year).
Because lmer uses a (effects | group) notation the random slope is specified using the (1 + Nges + Nges2 + NHerbst + Nges:NHerbst + Nges2:NHerbst | Jahr) term, which allows the slope of the relationship between Nges and yield to vary by year and Schlag, the latter being nested within Jahr. The term (“effect” | Jahr/Schlag), automatically expands to (“effect” | Jahr) + (1 | Jahr:Schlag).
The equation for the specified lmer model, which includes a “maximal” random effects structure for Jahr and for Schlag nested within Jahr, is as follows.
The model is specified on a single observation level as: \begin{aligned} \text{yield}_{ijk} = & \quad (\beta_0 + b_{0,i} + c_{0,ij}) \\ & + (\beta_1 + b_{1,i} + c_{1,ij}) \cdot \text{Nges}_{ijk} \\ & + (\beta_2 + b_{2,i} + c_{2,ij}) \cdot \text{NHerbst}_{ijk} \\ & + (\beta_3 + b_{3,i} + c_{3,ij}) \cdot \text{Nges2}_{ijk} \\ & + (\beta_4 + b_{4,i} + c_{4,ij}) \cdot (\text{Nges} \cdot \text{NHerbst})_{ijk} \\ & + (\beta_5 + b_{5,i} + c_{5,ij}) \cdot (\text{Nges2} \cdot \text{NHerbst})_{ijk} \\ & + \epsilon_{ijk} \end{aligned} \tag{3}
where: * i is the index for the year (Jahr). * j is the index for the replicate (Schlag) within year i. * k is the index for the individual observation. * \beta_0, \dots, \beta_5 are the fixed-effect coefficients, representing the population average intercept and slopes. * b_{0,i}, \dots, b_{5,i} are the random effects for Jahr i. They represent the deviation of year i’s intercept and slopes from the population average. * c_{0,ij}, \dots, c_{5,ij} are the random effects for Schlag j within Jahr i. They represent the deviation of that specific replicate’s intercept and slopes from the average for year i. * \epsilon_{ijk} is the residual error for each observation.
The random effects and residuals are assumed to follow a normal distribution with a mean of zero: \begin{bmatrix} b_{0,i} \\ \vdots \\ b_{5,i} \end{bmatrix} \sim \mathcal{N} \left( \mathbf{0}, \mathbf{\Sigma_b} \right) \quad , \quad \begin{bmatrix} c_{0,ij} \\ \vdots \\ c_{5,ij} \end{bmatrix} \sim \mathcal{N} \left( \mathbf{0}, \mathbf{\Sigma_c} \right) \quad , \quad \epsilon_{ijk} \sim \mathcal{N}(0, \sigma^2) where \mathbf{\Sigma_b} and \mathbf{\Sigma_c} are the variance-covariance matrices for the random effects at the Jahr and Jahr:Schlag levels, respectively.
library(lme4)
Yield.mod.lmer <- lmer(yield ~ Nges*NHerbst + Nges2*NHerbst +
(1 + Nges + Nges2 + NHerbst + Nges:NHerbst + Nges2:NHerbst | Jahr/Schlag),# +
data = dat
)4.4.5.1 Anova
The Anova of the “Yield.mod.lmer” gives no significant interaction between Nges and NHerbst (Table 7).
| Chisq | Df | Pr(>Chisq) | |
|---|---|---|---|
| (Intercept) | 32.1104534 | 1 | 0.0000000 |
| Nges | 189.6833891 | 1 | 0.0000000 |
| NHerbst | 9.2588144 | 1 | 0.0023436 |
| Nges2 | 60.4207699 | 1 | 0.0000000 |
| Nges:NHerbst | 0.8677199 | 1 | 0.3515873 |
| NHerbst:Nges2 | 0.2852191 | 1 | 0.5933006 |
4.4.5.2 Extracted parameters and economic optimum N rate and yield
Also for the lmer model the fixed effects can be extracted and the economic optimum N rate and yield can be calculated. The results are shown in Table 8. Again the differences in the economic optimum N rate and yield between the two variants with and without autumn N application are rather small.
| WithAutumnN | NoAutumnN | |
|---|---|---|
| Nmax | 246.540 | 244.968 |
| Ymax | 4.126 | 3.911 |
| Nopt_q | 211.633 | 212.967 |
| Yopt_q | 4.084 | 3.872 |
| (Intercept) | 2.036 | 1.660 |
| Nges | 0.017 | 0.018 |
| Nges2 | 0.000 | 0.000 |
4.4.5.3 Prediction plot from the mixed effect model lmer
#unique(as.numeric(row.names(ranef(Yield.mod.lmer)$ID)))
grid <- expand.grid(Jahr = unique(dat$Jahr), Nges = seq(0, 300, by = 10), NHerbst = unique(dat$NHerbst), Schlag= unique(dat$Schlag))
grid$Nges2 <- grid$Nges^2
grid$Yield.pred <- predict(Yield.mod.lmer, newdata = grid, allow.new.levels=TRUE)
grid$yield <- grid$Yield.pred
grid <- grid %>% select(Jahr, Schlag, NHerbst, Nges, yield) %>% group_by(Jahr, NHerbst, Nges) %>% summarise(yield = mean(yield, na.rm = TRUE)) But the predicted yield values are more accurate than the ones from the mixed effect model M3, which is shown in Figure 3.
4.4.6 Alternative implementation of the mixed effect model with random slope with the glmmTMB function from the glmmTMB package
The glmmTMB function from the glmmTMB package is used to fit the same model as above. The model includes the fixed effects Nges, NHerbst and Nges2, as well as the interaction of Nges and NHerbst . The random effects are Schlag (field) and Jahr (year). Note that the interaction of Nges and Herbst had to be removed from the model because of convergence problems. The random slope is specified using the (1 | Jahr:Nges) term, which allows the slope of the relationship between Nges and yield to vary by year. The dispformula argument is used to specify that the variance of the residuals should be allowed to vary by year.
\begin{align*} Y_{ijk} &= \beta_0 + \beta_1 \cdot \text{Nges}_{ijk} + \beta_2 \cdot \text{NHerbst}_{ijk} + \beta_3 \cdot \text{Nges2}_{ijk} \\ &\quad + \beta_4 \cdot (\text{Nges}_{ijk} \cdot \text{NHerbst}_{ijk}) \\ &\quad + b_{0j} + b_{1j} \cdot \text{Nges}_{ijk} + u_{jk} + \varepsilon_{ijk} \end{align*} \tag{4}
Where:
- Y_{ijk}: observed yield for observation i in plot k nested in year j
- \beta_0, \dots, \beta_4: fixed-effect coefficients
- b_{0j} \sim \mathcal{N}(0, \sigma^2_{b_0}): random intercept for year j
- b_{1j} \sim \mathcal{N}(0, \sigma^2_{b_1}): random slope for Nges within year j
- u_{jk} \sim \mathcal{N}(0, \sigma^2_u): random intercept for plot k within year j (i.e. Jahr:Schlag)
- \varepsilon_{ijk} \sim \mathcal{N}(0, \sigma^2_j): residual error with heteroscedastic variance per year (modeled via
dispformula = ~ Jahr)
This model reflects:
A fixed-effects structure with main and interaction terms
A random-intercept and random-slope model for Nges per year
A nested random effect for plots (Schlag) within years
Heteroscedastic residuals modeled by yearlibrary(glmmTMB)
Yield.M3.glmmTMB <- glmmTMB(yield ~ Nges * NHerbst + Nges2 + # fixed effects
(Nges | Jahr) + (1 | Jahr:Schlag), # random effects
dispformula = ~ Jahr, # This models different variances per year
data = dat)The summary of the model Yield.M3.glmmTMB shows the fixed effects estimates and their significance (Table 9). The ANOVA table for the model shows that the fixed effects Nges, NHerbst and Nges2 are significant and in difference to model Yield.M2 also their interaction (Table 9).
Family: gaussian ( identity )
Formula:
yield ~ Nges * NHerbst + Nges2 + (Nges | Jahr) + (1 | Jahr:Schlag)
Dispersion: ~Jahr
Data: dat
AIC BIC logLik -2*log(L) df.resid
727.0 817.9 -345.5 691.0 1131
Random effects:
Conditional model:
Groups Name Variance Std.Dev. Corr
Jahr (Intercept) 7.081e-01 0.841484
Nges 5.388e-06 0.002321 -0.26
Jahr:Schlag (Intercept) 1.638e-02 0.127989
Residual NA NA
Number of obs: 1149, groups: Jahr, 9; Jahr:Schlag, 18
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 2.007e+00 2.849e-01 7.042 1.89e-12 ***
Nges 1.739e-02 9.288e-04 18.728 < 2e-16 ***
NHerbstwithout autumn N -3.167e-01 4.249e-02 -7.453 9.11e-14 ***
Nges2 -3.601e-05 1.701e-06 -21.171 < 2e-16 ***
Nges:NHerbstwithout autumn N 5.563e-04 2.704e-04 2.058 0.0396 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Dispersion model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.29526 0.06431 -20.141 < 2e-16 ***
Jahr2017 -0.12009 0.09204 -1.305 0.191996
Jahr2018 0.08247 0.08975 0.919 0.358164
Jahr2019 0.08327 0.09178 0.907 0.364225
Jahr2020 0.47906 0.09027 5.307 1.12e-07 ***
Jahr2021 0.01674 0.09070 0.185 0.853552
Jahr2022 0.36401 0.09102 3.999 6.35e-05 ***
Jahr2023 0.31478 0.09099 3.460 0.000541 ***
Jahr2024 -0.08089 0.08964 -0.902 0.366851
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 14.4.6.1 Plot the predicted yield from the mixed effect model M3 with glmmTMB
The plot of predicted and observed yield vs. N rate in spring shows a good quality of the predictions. The predictions are made for each year and field (Schlag) combination but the mean value of the Schlag variables are shown (Figure 4).
grid <- expand.grid(Jahr = unique(dat$Jahr), Nges = seq(0, 300, by = 10), NHerbst = unique(dat$NHerbst), Schlag= unique(dat$Schlag))
grid$Nges2 <- grid$Nges^2
grid$yield <- predict(Yield.M3.glmmTMB, newdata = grid, allow.new.levels=TRUE)
grid <- grid %>% select(Jahr, Schlag, NHerbst, Nges, yield) %>% group_by(Jahr, NHerbst, Nges) %>% summarise(yield = mean(yield, na.rm = TRUE))
4.4.6.2 Parameter estimates and economic optimum N rate and yield
With the model Yield.M3.glmmTMB the parameters can be extracted in a similar way as for the model Yield.M3. The economic optimum N rate and yield can be calculated from the fixed effects estimates. The results are shown in @tbl_ExtractParsQuadModelglmmTMBP. Now the differences in the economic optimum N rate and yield between the two variants with and without autumn N application are more clearly visible. The economic optimum N rate is higher for the variant without autumn N application and the yield at economic optimum is higher for the variant with autumn N application.
| WithAutumnN | NoAutumnN | |
|---|---|---|
| Nmax | 241.534 | 249.259 |
| Ymax | 4.107 | 3.927 |
| Nopt_q | 208.208 | 215.933 |
| Yopt_q | 4.067 | 3.887 |
| (Intercept) | 2.007 | 1.690 |
| Nges | 0.017 | 0.018 |
| Nges2 | 0.000 | 0.000 |
5 QP model approach
The former approaches fitted a linear quadratic model to the yield data. In this section we try to fit a quadratic plateau model, which assumes that the yield is not decreasing after reaching its maximum value at a certain N rate:
Y = \begin{cases} a + b \cdot \text{Nges} + c \cdot \text{Nges}^2, & \text{if } \text{Nges} \leq x_0 \\ Y_{\text{max}}, & \text{if } \text{Nges} > x_0 \end{cases} \tag{5}
where:
x_0 = -\frac{b}{2c}, \quad Y_{\text{max}} = a + b \cdot x_0 + c \cdot x_0^2 \tag{6}
As a statistical model with NHerbst a an additional fixed effect and Jahr and Schlag as nested random variables this leads to the following model equation:
Y_{ijk} = \begin{cases} a_{h} + b_{h} \cdot \text{Nges}_{ijk} + c_{h} \cdot \text{Nges}_{ijk}^2 + u_{0,j} + u_{1,j} \cdot \text{Nges}_{ijk}, & \text{if } \text{Nges}_{ijk} \leq x_{0,hj} \\ Y_{\text{max},hj}, & \text{if } \text{Nges}_{ijk} > x_{0,hj} \end{cases} \tag{7}
where:
x_{0,hj} = -\frac{b_{h} + u_{1,j}}{2(c_{h})}, \quad Y_{\text{max},hj} = a_{h} + b_{h} \cdot x_{0,hj} + c_{h} \cdot x_{0,hj}^2 + u_{0,j} + u_{1,j} \cdot x_{0,hj} \tag{8}
where: - i indexes observations, j the combination of year and plot (Jahr/Schlag), h the level of the factor NHerbst. - $ a_h, b_h, c_h$: fixed effect coefficients depending on NHerbst - $ u_{0,j}, u_{1,j} $: random intercept and slope for a and b per Jahr/Schlag - Nges_{ijk} ): nitrogen rate - Residual variances \sigma^2_\epsilon are heterogeneous across years:
\epsilon_{ijk} \sim \mathcal{N}(0, \sigma^2_{Jahr_k})
This model allows:
Fixed effects of parameters a,b,ca,b,c to vary with NHerbst.
Random intercepts and slopes (a, b) at Jahr/Schlag level.
Heteroscedastic residual variances per Jahr.
For fitting this model to the data the packages nlraa provides a self starting function SSquadp3 (https://femiguez.github.io/nlraa-docs/nlraa.html). Also the package nlme is used. In a first step the yield data frame is converted into a grouped data object. In a second step then first estimates of the paramter values are obtained.
library(nlraa)
# creating a grouped data object
tmpG <- groupedData(yield ~ Nges | NHerbst, data = dat[,c("yield", "Nges", "NHerbst","key", "Jahr", "Schlag")])
# fit models to the grouped data
fit3L <- nlsLMList(yield ~ SSquadp3(Nges, a, b, c), data = tmpG)
# adopt the model with random effects
fit3L.nlme <- nlme(fit3L, random==pdDiag(a+b+c~1))
# extract the fixed effects
fxf <- fixef(fit3L.nlme)
# nlme braucht Startwerte für alle Parameter, für die jeweils 3 "zusätzlichen" Effekte der anderen FF ist der Startwert erst mal 0
# nlme needs starting values for the three parameters, obtained from the previous estimated fixed effects, the effect of NHerbst is initially
# set to zero.
mod_base <- nlme(
yield ~ SSquadp3(Nges, a, b, c),
fixed = a + b + c ~ NHerbst,
random = a + b ~ 1 | Jahr/Schlag,
start = c(fxf[1], rep(0,1), fxf[2], rep(0,1), fxf[3], rep(0,1)),
data = tmpG,
control = nlmeControl(maxIter = 500, msMaxIter = 500, tolerance = 1e-6)
)
mod_het <- update(mod_base,
weights = varIdent(form = ~1 | Jahr),
control = nlmeControl(maxIter = 500, msMaxIter = 500, tolerance = 1e-6)
)5.0.0.1 Anova
The model mod_het is fitted with heteroscedasticity, allowing the variance to differ by year. The anova table for this model shows the significance of the fixed effects and their interactions. There is a significant value for the effect of NHerbst on the slope of the quadratic plateau function (?@tbl-AnovaQPmodel).
| numDF | denDF | F-value | p-value | |
|---|---|---|---|---|
| a.(Intercept) | 1 | 1126 | 129.282250 | 0.0000000 |
| a.NHerbst | 1 | 1126 | 180.509353 | 0.0000000 |
| b.(Intercept) | 1 | 1126 | 95.333315 | 0.0000000 |
| b.NHerbst | 1 | 1126 | 6.043529 | 0.0141071 |
| c.(Intercept) | 1 | 1126 | 327.757486 | 0.0000000 |
| c.NHerbst | 1 | 1126 | 1.394455 | 0.2379036 |
5.1 Parameter estimates and economic optimum N rate and yield
The parameter estimates and the optimum N rates and yield levels are shown in Table 12.
| WithAutumN | NoAutumnN | |
|---|---|---|
| Nmax | 241.638 | 237.516 |
| Ymax | 4.135 | 3.984 |
| Nopt_q | 207.150 | 206.101 |
| Yopt_q | 4.094 | 3.946 |
| intercept | 2.104 | 1.829 |
| slope | 0.017 | 0.018 |
| quad | 0.000 | 0.000 |
5.1.1 Plot the predicted yield from the quadratic plateau model with heteroscedasticity
The qp-model fitted the experimental data quite well (Figure 5).
grid <- expand.grid(Jahr = unique(dat$Jahr), Nges = seq(0, 300, by = 10), NHerbst = unique(dat$NHerbst), Schlag= unique(dat$Schlag))
grid$Nges2 <- grid$Nges^2
grid$yield <- predict(mod_het, newdata = grid, allow.new.levels=TRUE)
grid <- grid %>% select(Jahr, Schlag, NHerbst, Nges, yield) %>% group_by(Jahr, NHerbst, Nges) %>% summarise(yield = mean(yield, na.rm = TRUE))