Getting started with the ALFAM2 package

Introduction

The ALFAM2 project is on ammonia volatilization (emission) from field-applied manure, and includes two main products: a database with volatilization measurements and a model for estimating volatilization. The model, which is described in Hafner et al. (2019)1 and Hafner et al. (2024b)2, is the focus of the ALFAM2 R package and this document. The ALFAM2 package is an add-on package for R, which is an environment for statistical computing. With the model, it is possible to predict average volatilization rate and cumulative emission over time, as affected by application method, weather, and manure characteristics. This document provides an introduction to the use of the model, and is aimed for users who are new to R. Readers with some knowledge of R can skip the next section.

Excel or R?

The ALFAM2 model is available in an Excel spreadsheet in addition to the R package that is described in this document. For users who would just like to know cumulative emission for a few scenarios with constant conditions, the Excel model is a good choice. But to work with many different scenarios, or when weather changes over time (e.g., wind or rain), or if you are interested in emission dynamics and not just final cumulative emission, you should use the R package. You can use the ALFAM2 package without much knowledge of R. For anyone planning to use the ALFAM2 model extensively, it may be worthwhile to learn a little bit about R and use the ALFAM2 package, instead of the less efficient Excel spreadsheet model.

Some basics for new R users

The information given in this document should be enough for new R users to install the package and learn enough about R to start using the model (albeit with a lack of understanding about some of the code). For a better understanding, check out the sources mentioned below.

To use R it must be installed on your computer. You can download R and find installation instructions from here: https://www.r-project.org/. And while not required, it is convenient to have a good script editor. The RStudio IDE (integrated development environment) is a popular choice. It can be downloaded from here: https://posit.co/download/rstudio-desktop/.

To use the ALFAM2 package, you will need to install the package (done once, or every time an updated version is needed), load the package (done once in each R session), and then call up the alfam2() function. In R, you will need to be able to install and load packages, call functions, and, ideally, create data frames and import and export data to/from file. For information on these tasks and more, there are many free resources online. CRAN provides various manuals, including a good introduction: https://cran.r-project.org/ (select “Manuals” at the lower left). Posit (RStudio developers) also provides various materials for learning R, although the focus is skewed toward packages from Posit. This book for a course on R has some details: https://www.researchgate.net/publication/325170649_An_Introduction_to_R_for_Beginners. Alternatively, the examples given below may be sufficient.

Installing the ALFAM2 package

Because the ALFAM2 package is on CRAN, it can be installed in the normal way.

install.packages('ALFAM2')

This is the recommended approach.

The latest and earlier versions of the ALFAM2 package are available from a GitHub repository, which includes installation instructions: https://github.com/AU-BCE-EE/ALFAM2.

Every time you open R to use the ALFAM2 package, it must be loaded.

library(ALFAM2)

You can open this vignette with the following call.

vignette("ALFAM2-start")

To check the version of the installed package, use this:

packageVersion("ALFAM2")
## [1] '4.2'

The alfam2() function

The ALFAM2 package includes a single function that is an implementation of the ALFAM2 model: alfam2(). After an explanation of the function, its use is shown in a few examples.

Overview of the function

The alfam2() function can be used for predicting average volatilization rate and cumulative emission over time. The function has several arguments, as shown below.

args(alfam2)
## function (dat, pars = ALFAM2::alfam2pars03, add.pars = NULL, 
##     app.name = "TAN.app", time.name = "ct", time.incorp = NULL, 
##     group = NULL, center = c(app.rate = 40, man.dm = 6, man.tan = 1.2, 
##         man.ph = 7.5, air.temp = 13, wind.2m = 2.7, wind.sqrt = sqrt(2.7), 
##         crop.z = 10), pass.col = NULL, incorp.names = c("incorp", 
##         "deep", "shallow"), prep.dum = TRUE, prep.incorp = TRUE, 
##     add.incorp.rows = FALSE, check = TRUE, warn = TRUE, value = "emis", 
##     conf.int = NULL, pars.ci = ALFAM2::alfam2pars03var, n.ci = NULL, 
##     var.ci = "er", ...) 
## NULL

You can find more details on the arguments (as well as examples) in the help file. As with any R function, you can open the file with ?:

?alfam2

But the most important arguments are described here. Most arguments have default values, and the only one that is required to make predictions is the dat argument, which is a data frame containing some input data, i.e., values of predictor variables and time after slurry application. The dat data frame can contain any number of rows (each representing a time interval), but must contain a column with cumulative time in hours, and the name of this column is indicated with time.name. Typically the data frame will have predictor variables as well, for example, manure dry matter, application method, air temperature, or wind speed. The name of the predictor columns are used to link predictor variables to model parameters, which are set by the pars argument. Usually the default parameter values, based on the measurements in the ALFAM2 database, should be used. Predictor variables used in the available default parameter sets and their default names are given in Table 1 below.

Variable name Description Units Notes
int Intercept terms None
app.mthd.os Open slot application None (logical) Binary variable
app.mthd.cs Closed slot application None (logical) Binary variable
app.mthd.bc Broadcast application None (logical) Binary variable
app.mthd.ts Trailing shoe application None (logical) Binary variable
app.rate Manure application rate t/ha
app.rate.ni Manure app. rate (no inj.) t/ha
man.dm Manure dry matter %
man.ph Manure pH pH units For acidification
man.source.pig Pig manure None (logical) Binary variable
incorp.deep Deep incorporation None (logical) Binary variable
incorp.shallow Shallow incorporation None (logical) Binary variable
air.temp Air temperature ∘C
wind.2m Wind speed (at 2 m) m/s
wind.sqrt Square root of wind speed (2 m) m1/2/s1/2
rain.rate Rainfall rate mm/h
rain.cum Cumulative rain mm
cereal.hght Cereal height cm

Binary variables are dummy variables, and can be created automatically by the alfam2 function (see examples below). Default model parameters and numeric values in the latest set (currently alfam2pars03 object (“Set 3”)) should generally be used. Development of this parameter set is described in Hafner et al. (2024b)3. There are two earlier versions: “Set 1” is available in alfam2pars01, and “Set 2” in alfam2pars02. Set 2 is described in a report on calculation of Danish emission factors4 and a later paper5. And set 1 in the 2019 paper listed below6.

Because values of predicted emission depend on parameter values (strongly depend, in some cases), it is essential that the any application of the model lists the parameter set that was used. Additional, it is good practice to list the version of the package. Before finalizing reports or papers users can check for these details like this:

packageVersion("ALFAM2")
## [1] '4.2'

and

args(alfam2)
## function (dat, pars = ALFAM2::alfam2pars03, add.pars = NULL, 
##     app.name = "TAN.app", time.name = "ct", time.incorp = NULL, 
##     group = NULL, center = c(app.rate = 40, man.dm = 6, man.tan = 1.2, 
##         man.ph = 7.5, air.temp = 13, wind.2m = 2.7, wind.sqrt = sqrt(2.7), 
##         crop.z = 10), pass.col = NULL, incorp.names = c("incorp", 
##         "deep", "shallow"), prep.dum = TRUE, prep.incorp = TRUE, 
##     add.incorp.rows = FALSE, check = TRUE, warn = TRUE, value = "emis", 
##     conf.int = NULL, pars.ci = ALFAM2::alfam2pars03var, n.ci = NULL, 
##     var.ci = "er", ...) 
## NULL

(See the default for the pars argument, alfam2pars03.)

So in this case, we would write “…ALFAM2 R package (version 4.2) was used with the default parameter set 3…”.

Comparing the contents of alfam2pars03 to the variable names given in Table 1, you can see an additional letter and number added to the end of the parameters.

alfam2pars03
##            int.f0    app.mthd.os.f0    app.mthd.cs.f0 
##        0.45305451       -2.89718049       -7.09642528 
## man.source.pig.f0         man.dm.f0            int.r1 
##       -0.95213804        0.49956176       -1.45119862 
##    app.mthd.bc.r1    app.mthd.ts.r1         man.dm.r1 
##        0.73714111       -0.07393662       -0.03300931 
##         man.ph.r1       air.temp.r1      wind.sqrt.r1 
##        0.42121280        0.03321186        0.46104870 
##            int.r2      rain.rate.r2            int.r3 
##       -1.16953266        0.60163865       -2.68829766 
##    app.mthd.cs.r3    incorp.deep.r3         man.ph.r3 
##       -0.38439637       -5.35112099        0.11776977 
## incorp.shallow.f4    incorp.deep.f4            int.r5 
##       -1.41820869       -2.94966810       -1.80000000 
##      rain.rate.r5 
##        0.48425409

These numbers indicate a primary parameter. So, for example, the (secondary) parameter air.temp.r1 is used in the calculation of the primary parameter r1. The most important message here is a simple one: names for predictor variables can be taken from the names given in the default pars argument value, but be sure to omit the last three characters (a “.”, a number, and a letter).

By design, any time a predictor variable is omitted when alfam2() is called, the reference level or value is assumed for that variable (one exception is app.rate.ni). The scenario with reference levels for all predictors is the default scenario, and is the one given in the first row of Table 4 in 7. Predictor values for the default scenario can be found in the center argument. The default application method is trailing hose. The center argument is used for centering predictor variables, and this approach facilities the behavior described above.

Cumulative emission for a single scenario

In this example, let’s assume we are interested in broadcast application of manure that has 8% dry matter (DM), and TAN application of 50 kg/ha. Here wind was 3 m/s, and air temperature 20∘C.

First we need to create a data frame with the input data.

dat1 <- data.frame(ctime = 168, TAN.app = 50, man.dm = 8, 
                   air.temp = 20, wind.sqrt = 2, 
                   app.mthd = 'bc')
print(dat1)
##   ctime TAN.app man.dm air.temp wind.sqrt app.mthd
## 1   168      50      8       20         2       bc

Our predictor variable values are in the columns man.dm and the following ones. The names for the predictor variables must match those names used in the model parameters, which can be seen by checking the parameter object contents (see just above). Here and in most of the example we will let the alfam2 function create our application method dummy variable for us by adding a column app.mthd = 'bc'. Alternatively, we could use app.mthd.bc = TRUE or app.mthd.bc = 1 to add it manually. For this example neither approach is clearly easier, but with combinations that include different application methods, it is easier to let the alfam2 function do the work, as in the first example above. See the section on dummy variables for more details.

Time, in hours after application, is given in the column named ctime here, for cumulative time (although any name can be used).

And now we can call the model function, using default values for most other arguments. We can predict cumulative emission after 7 days (168 hours) with the following call.

pred1 <- alfam2(dat1, app.name = 'TAN.app', time.name = 'ctime')
## Default parameters (Set 3) are being used.
## Warning in alfam2(dat1, app.name = "TAN.app", time.name = "ctime"): Incorporation columns  were dropped 
##     because argument time.incorp is NULL
##     So there is no incorporation.
##     Set check = FALSE to not drop, but then check output.
## Warning in alfam2(dat1, app.name = "TAN.app", time.name = "ctime"): Running with 14 parameters. Dropped 8 with no match.
## These secondary parameters have been dropped:
##   man.source.pig.f0
##   man.ph.r1
##   rain.rate.r2
##   incorp.deep.r3
##   man.ph.r3
##   incorp.shallow.f4
##   incorp.deep.f4
##   rain.rate.r5

The warning messages tell us that the call included some parameters with no associated predictor variables in our data frame given in the dat argument. This is discussed more below. We will turn off the warnings in some examples below.

Let’s look at the predictions.

print(pred1)
##   app.mthd.ts app.mthd.bc app.mthd.os app.mthd.cs ctime  dt
## 1           0           1           0           0   168 168
##              f         s        e       ei         j        er
## 1 2.942948e-34 0.7612382 36.47372 36.47372 0.2171055 0.7294744
##          f0        r1         r2          r3 f4         r5
## 1 0.8103334 0.4139272 0.06768109 0.002049757  1 0.01584893
##         jinst
## 1 0.001560353

The most interesting columns here are called e, which has cumulative emission in the same units as TAN application, and er, which has relative cumulative emission, as a fraction of applied TAN. So in this example, 73% of applied TAN is predicted to be lost by volatilization.

The warning message above is related to an important point: Any excluded predictors are effectively assumed to be at their reference levels.

Application method names

The secondary parameters and dummy variables use two-letter codes for application methods: 'bc' (broadcast), 'ts' (trailing shoe), 'os' (open slot injection), and 'cs' (closed slot injection). Trailing hose is the reference method. These can be used directly as values in a app.mthd column in input data. But there are some aliases that can be used in the app.mthd column: 'broadcast', and 'broadspread' for 'bc'; 'trailing shoe' for 'ts'; 'open slot injection', 'open-slot injection', and 'shallow injection' for 'os'; and 'closed slot injection', 'closed-slot injection', and 'deep injection' for 'cs'. These are not case sensitive. So the following is effectively the same as the previous example.

dat1b <- data.frame(ctime = 168, TAN.app = 50, man.dm = 8, 
                   air.temp = 20, wind.sqrt = 2, 
                   app.mthd = 'Broadcast')
print(dat1b)
##   ctime TAN.app man.dm air.temp wind.sqrt  app.mthd
## 1   168      50      8       20         2 Broadcast
pred1b <- alfam2(dat1, app.name = 'TAN.app', time.name = 'ctime')
## Default parameters (Set 3) are being used.
## Warning in alfam2(dat1, app.name = "TAN.app", time.name = "ctime"): Incorporation columns  were dropped 
##     because argument time.incorp is NULL
##     So there is no incorporation.
##     Set check = FALSE to not drop, but then check output.
## Warning in alfam2(dat1, app.name = "TAN.app", time.name = "ctime"): Running with 14 parameters. Dropped 8 with no match.
## These secondary parameters have been dropped:
##   man.source.pig.f0
##   man.ph.r1
##   rain.rate.r2
##   incorp.deep.r3
##   man.ph.r3
##   incorp.shallow.f4
##   incorp.deep.f4
##   rain.rate.r5
all.equal(pred1, pred1b)
## [1] TRUE

The two-letter codes are probably safer, but in any case the output can be checked for the expected dummy variables.

Adding incorporation

To include incorporation of slurry into the soil after application, we need to add a couple columns to our data frame. First let’s make a new data frame for the example.

dat2 <- dat1

And add the two new columns. Here we are specifying that deep incorporation happens after 0.5 hours.

dat2$incorp <- 'deep'
dat2$t.incorp <- 0.5
print(dat2)
##   ctime TAN.app man.dm air.temp wind.sqrt app.mthd incorp
## 1   168      50      8       20         2       bc   deep
##   t.incorp
## 1      0.5
pred2 <- alfam2(dat2, app.name = "TAN.app", time.name = "ctime", 
                   time.incorp = "t.incorp", warn = FALSE)
print(pred2)
##   app.mthd.ts app.mthd.bc app.mthd.os app.mthd.cs incorp.shallow
## 1           0           1           0           0              0
##   incorp.deep ctime  dt            f        s        e       ei
## 1           1   168 168 1.464181e-35 2.890399 8.824327 8.824327
##            j        er        f0        r1         r2
## 1 0.05252576 0.1764865 0.8103334 0.4139272 0.06768109
##             r3        f4         r5        jinst
## 1 9.132325e-09 0.0497522 0.01584893 2.639606e-08

Here we see that with incorporation, emission drops to 18% of applied TAN. Shallow incorporation has less of an effect.

dat3 <- dat1
dat3$incorp <- 'shallow'
dat3$t.incorp <- 0.5
print(dat3)
##   ctime TAN.app man.dm air.temp wind.sqrt app.mthd  incorp
## 1   168      50      8       20         2       bc shallow
##   t.incorp
## 1      0.5
pred3 <- alfam2(dat3, app.name = "TAN.app", time.name = "ctime", 
                   time.incorp = "t.incorp", warn = FALSE)
print(pred3)
##   app.mthd.ts app.mthd.bc app.mthd.os app.mthd.cs incorp.shallow
## 1           0           1           0           0              1
##   incorp.deep ctime  dt            f        s        e       ei
## 1           0   168 168 5.737058e-35 1.853516 16.83718 16.83718
##           j        er        f0        r1         r2          r3
## 1 0.1002213 0.3367437 0.8103334 0.4139272 0.06768109 0.002049757
##          f4         r5       jinst
## 1 0.1949426 0.01584893 0.003799257

Relative emission only

By default the alfam2() function expects TAN application rate to be in a column named TAN.app. Different names can be specified using the app.name argument. If the specified name is not present in input data frame, the function returns relative emission only (emission as a fraction of applied TAN), with a warning. This can be convenient when only relative (not absolute) emission is needed.

datr <- data.frame(ctime = 168)
print(datr)
##   ctime
## 1   168
predr <- alfam2(datr, app.name = 'TAN.app', time.name = 'ctime')
## Warning in alfam2(datr, app.name = "TAN.app", time.name = "ctime"): Argument app.name is missing or dat is missing column of given name.
##     So function will return relative emission only.
## Default parameters (Set 3) are being used.
## Warning in prepDat(dat, warn = warn): Argument prep.dum = TRUE but there are no variables to convert to dummy variables!
##   Ignoring prep.dum = TRUE.
## Warning in alfam2(datr, app.name = "TAN.app", time.name = "ctime"): Incorporation columns  were dropped 
##     because argument time.incorp is NULL
##     So there is no incorporation.
##     Set check = FALSE to not drop, but then check output.
## Warning in alfam2(datr, app.name = "TAN.app", time.name = "ctime"): Running with 5 parameters. Dropped 17 with no match.
## These secondary parameters have been dropped:
##   app.mthd.os.f0
##   app.mthd.cs.f0
##   man.source.pig.f0
##   man.dm.f0
##   app.mthd.bc.r1
##   app.mthd.ts.r1
##   man.dm.r1
##   man.ph.r1
##   air.temp.r1
##   wind.sqrt.r1
##   rain.rate.r2
##   app.mthd.cs.r3
##   incorp.deep.r3
##   man.ph.r3
##   incorp.shallow.f4
##   incorp.deep.f4
##   rain.rate.r5
predr
##   ctime  dt            f          s         e        ei
## 1   168 168 1.847366e-08 0.04323519 0.2954223 0.2954223
##             j        er        f0         r1         r2
## 1 0.001758466 0.2954223 0.6113652 0.03538355 0.06768109
##            r3 f4         r5        jinst
## 1 0.002049757  1 0.01584893 8.862227e-05

Notice how the column e (emission) is equal to er (relative emission). Even flux is relative, in fraction of applied TAN per hour, so multiplying it by the single interval duration will give total relative emission.

predr$j * predr$dt
## [1] 0.2954223

With this functionality, a minimal example of dat is a data frame with only cumulative time, as in datr above.

Multiple scenarios in a single call

A single function call can be used for multiple scenarios. For example, perhaps we would like to compare 5 different incorporation times. First let’s create a new data frame that contains the inputs. We will need to add a new column with a grouping variable also, to let alfam2() know that each row represents a different scenario.

dat4 <- data.frame(scenario = 1:5, ctime = 168, TAN.app = 50, 
                   man.dm = 8, air.temp = 20, wind.sqrt = 2, 
                   app.mthd = 'bc',
                   incorp = 'deep',
                   t.incorp = c(0.1, 1, 6, 24, NA))
print(dat4)
##   scenario ctime TAN.app man.dm air.temp wind.sqrt app.mthd
## 1        1   168      50      8       20         2       bc
## 2        2   168      50      8       20         2       bc
## 3        3   168      50      8       20         2       bc
## 4        4   168      50      8       20         2       bc
## 5        5   168      50      8       20         2       bc
##   incorp t.incorp
## 1   deep      0.1
## 2   deep      1.0
## 3   deep      6.0
## 4   deep     24.0
## 5   deep       NA

It may seem strange to have a scenario column–isn’t it clear that each row is a different scenario? The answer is no, not when there are multiple time intervals per scenario, for example when one is interested in volatilization rates over time and how they change (see the section on dynamics). Also, it is necessary to submit multiple rows per scenario in order to consider changing weather over time. Note that there is no incorporation for scenario 5–the last row in the data frame. We could also specify this behavior with t.incorp = Inf. To be even more clear, we might use the following, although the alfam2 function will return the same (correct) results for either version of the data frame.

dat4 <- data.frame(scenario = 1:5, ctime = 168, TAN.app = 50, 
                   man.dm = 8, air.temp = 20, wind.sqrt = 2, 
                   app.mthd = 'bc',
                   incorp = c(rep('deep', 4), 'none'),
                   t.incorp = c(0.1, 1, 6, 24, NA))
print(dat4)
##   scenario ctime TAN.app man.dm air.temp wind.sqrt app.mthd
## 1        1   168      50      8       20         2       bc
## 2        2   168      50      8       20         2       bc
## 3        3   168      50      8       20         2       bc
## 4        4   168      50      8       20         2       bc
## 5        5   168      50      8       20         2       bc
##   incorp t.incorp
## 1   deep      0.1
## 2   deep      1.0
## 3   deep      6.0
## 4   deep     24.0
## 5   none       NA

Let’s run the model for these 5 scenarios.

pred4 <- alfam2(dat4, app.name = "TAN.app", time.name = "ctime", 
                   time.incorp = "t.incorp", group = "scenario", 
                   warn = FALSE)
print(pred4)
##   scenario app.mthd.ts app.mthd.bc app.mthd.os app.mthd.cs
## 1        1           0           1           0           0
## 2        2           0           1           0           0
## 3        3           0           1           0           0
## 4        4           0           1           0           0
## 5        5           0           1           0           0
##   incorp.shallow incorp.deep ctime  dt            f         s
## 1              0           1   168 168 1.464181e-35 3.2634462
## 2              0           1   168 168 1.464181e-35 2.5117474
## 3              0           1   168 168 1.464181e-35 1.2012935
## 4              0           1   168 168 1.464181e-35 1.0226386
## 5              0           0   168 168 2.942948e-34 0.7612382
##           e        ei          j        er        f0        r1
## 1  3.290395  3.290395 0.01958568 0.0658079 0.8103334 0.4139272
## 2 14.401664 14.401664 0.08572419 0.2880333 0.8103334 0.4139272
## 3 33.138914 33.138914 0.19725544 0.6627783 0.8103334 0.4139272
## 4 35.413075 35.413075 0.21079211 0.7082615 0.8103334 0.4139272
## 5 36.473720 36.473720 0.21710548 0.7294744 0.8103334 0.4139272
##           r2           r3        f4         r5        jinst
## 1 0.06768109 9.132325e-09 0.0497522 0.01584893 2.980285e-08
## 2 0.06768109 9.132325e-09 0.0497522 0.01584893 2.293809e-08
## 3 0.06768109 9.132325e-09 0.0497522 0.01584893 1.097060e-08
## 4 0.06768109 9.132325e-09 0.0497522 0.01584893 9.339068e-09
## 5 0.06768109 2.049757e-03 1.0000000 0.01584893 1.560353e-03

We can see that predicted emission increases substantially as incorporation time goes up.

barplot(pred4$er, names.arg = paste(dat4$t.incorp), xlab = 't.incorp', ylab = 'er')

And incorporation after 24, or really even 6, hours is not much better than no incorporation.

Scenarios could differ in any way. Below, both temperature and application method vary.

dat5 <- data.frame(scenario = 1:3, ctime = 168, TAN.app = 50, 
                   man.dm = 8, wind.sqrt = 2,
                   air.temp = c(15, 20, 25),
                   app.mthd = c('bc', 'bsth', 'os')
                   )
print(dat5)
##   scenario ctime TAN.app man.dm wind.sqrt air.temp app.mthd
## 1        1   168      50      8         2       15       bc
## 2        2   168      50      8         2       20     bsth
## 3        3   168      50      8         2       25       os
pred5 <- alfam2(dat5, app.name = "TAN.app", time.name = "ctime", 
                   group = "scenario", warn = FALSE)
print(pred5)
##   scenario app.mthd.ts app.mthd.bc app.mthd.os app.mthd.cs ctime
## 1        1           0           1           0           0   168
## 2        2           0           0           0           0   168
## 3        3           0           0           1           0   168
##    dt            f         s        e       ei         j
## 1 168 1.161991e-24 0.8770038 34.56623 34.56623 0.2057514
## 2 168 1.372760e-09 1.5482796 24.50443 24.50443 0.1458597
## 3 168 8.569609e-13 2.1987970 10.72352 10.72352 0.0638305
##          er        f0         r1         r2          r3 f4
## 1 0.6913246 0.8103334 0.28239999 0.06768109 0.002049757  1
## 2 0.4900885 0.8103334 0.07581984 0.06768109 0.002049757  1
## 3 0.2144705 0.1907719 0.11113278 0.06768109 0.002049757  1
##           r5       jinst
## 1 0.01584893 0.001797645
## 2 0.01584893 0.003173597
## 3 0.01584893 0.004506999

Dummy variables

Most of the examples given above include “dummy variables”, or multiple binary variables that together represent a factor. The package (≥v1.2) includes a helper function to calculate these dummy variables internally, making data entry much simpler (and less error-prone) in some cases. The new columns can be seen in the output from the function–for example app.mthd.bc and the similar columns in pred5 above. Dummy variables can still be entered manually though. For example, we could duplicate the behavior of the call above with the following.

dat5b <- data.frame(scenario = 1:3, ctime = 168, TAN.app = 50, 
                   man.dm = 8, wind.sqrt = 2,
                   air.temp = c(15, 20, 25),
                   app.mthd.bc = c(TRUE, FALSE, FALSE),
                   app.mthd.os = c(FALSE, FALSE, TRUE)
                   )
print(dat5b)
##   scenario ctime TAN.app man.dm wind.sqrt air.temp app.mthd.bc
## 1        1   168      50      8         2       15        TRUE
## 2        2   168      50      8         2       20       FALSE
## 3        3   168      50      8         2       25       FALSE
##   app.mthd.os
## 1       FALSE
## 2       FALSE
## 3        TRUE

For scenario 2, application method is not explicitly specified, which means it is the default–trailing hose.

pred5b <- alfam2(dat5b, app.name = "TAN.app", time.name = "ctime", 
                    group = "scenario", warn = FALSE)
print(pred5b)
##   scenario ctime  dt            f         s        e       ei
## 1        1   168 168 1.161991e-24 0.8770038 34.56623 34.56623
## 2        2   168 168 1.372760e-09 1.5482796 24.50443 24.50443
## 3        3   168 168 8.569609e-13 2.1987970 10.72352 10.72352
##           j        er        f0         r1         r2
## 1 0.2057514 0.6913246 0.8103334 0.28239999 0.06768109
## 2 0.1458597 0.4900885 0.8103334 0.07581984 0.06768109
## 3 0.0638305 0.2144705 0.1907719 0.11113278 0.06768109
##            r3 f4         r5       jinst
## 1 0.002049757  1 0.01584893 0.001797645
## 2 0.002049757  1 0.01584893 0.003173597
## 3 0.002049757  1 0.01584893 0.004506999
all.equal(pred5b$e, pred5$e)
## [1] TRUE

Calculated emission is the same as in pred5 above, but in pred5 we also have the dummy variables returned in the output. It is good practice to check these columns to make sure inputs were entered correctly.

This data preparation option invoked by prep.dum = TRUE, which is the default. The alfam2 function will automatically add dummy variables for application method, incorporation, and manure source (currently pig is the only level different from the reference). For this automatic conversion to take place, column names and factor levels must match the relevant part of the parameter names. In the call immediately above, app.mthd matches the beginning of the parameter names app.mthd.bc etc., and the levels bc and os match the final part of these names. The level th has no match–it is (correctly) interpreteted as a reference level.

Here is an example of automatic dummy variable processing with incorporation. Let’s check out the incorporation parameters first:

alfam2pars03[grepl('^incorp', names(alfam2pars03))]
##    incorp.deep.r3 incorp.shallow.f4    incorp.deep.f4 
##         -5.351121         -1.418209         -2.949668

The following example includes both shallow and deep incorporation combined with trailing hose or broadcast application as well as the slurry source. To help ensure that incorporation is indeed applied as intended we set warn = TRUE (because the information entered in the data frame is inconsistent for rows 1 and 4, with incorp set to 'none' but t.incorp set to 4, which will not cause problems, but checking is good practice).

dat6 <- data.frame(scenario = 1:6, ctime = 168, TAN.app = 100, 
                   man.dm = 5, man.ph = 7.2, air.temp = 10, 
                   wind.sqrt = 2, 
                   man.source = c(rep('Cattle', 2), rep('Pig', 4)),
                   app.mthd = rep(c('Broadcast', 'Trailing hose'), 
                                  each = 3),
                   incorp = rep(c('None', 'Shallow', 'Deep'), 2),
                   t.incorp = 4)
print(dat6)
##   scenario ctime TAN.app man.dm man.ph air.temp wind.sqrt
## 1        1   168     100      5    7.2       10         2
## 2        2   168     100      5    7.2       10         2
## 3        3   168     100      5    7.2       10         2
## 4        4   168     100      5    7.2       10         2
## 5        5   168     100      5    7.2       10         2
## 6        6   168     100      5    7.2       10         2
##   man.source      app.mthd  incorp t.incorp
## 1     Cattle     Broadcast    None        4
## 2     Cattle     Broadcast Shallow        4
## 3        Pig     Broadcast    Deep        4
## 4        Pig Trailing hose    None        4
## 5        Pig Trailing hose Shallow        4
## 6        Pig Trailing hose    Deep        4
pred6 <- alfam2(dat6, app.name = "TAN.app", time.name = "ctime", 
                time.incorp = "t.incorp", group = "scenario", 
                warn = TRUE)
## Default parameters (Set 3) are being used.
## Incorporation applied for groups: 2, 3, 5, 6.
## Warning in alfam2(dat6, app.name = "TAN.app", time.name = "ctime", time.incorp = "t.incorp", : Running with 20 parameters. Dropped 2 with no match.
## These secondary parameters have been dropped:
##   rain.rate.r2
##   rain.rate.r5
print(pred6)
##   scenario app.mthd.ts app.mthd.bc app.mthd.os app.mthd.cs
## 1        1           0           1           0           0
## 2        2           0           1           0           0
## 3        3           0           1           0           0
## 4        4           0           0           0           0
## 5        5           0           0           0           0
## 6        6           0           0           0           0
##   incorp.shallow incorp.deep man.source.pig ctime  dt
## 1              0           0              0   168 168
## 2              1           0              0   168 168
## 3              0           1              1   168 168
## 4              0           0              1   168 168
## 5              1           0              1   168 168
## 6              0           1              1   168 168
##              f        s         e        ei          j
## 1 3.556088e-17 3.325689 42.053470 42.053470 0.25031827
## 2 6.932330e-18 3.886253 32.535966 32.535966 0.19366646
## 3 9.752566e-19 6.103600 13.257164 13.257164 0.07891169
## 4 1.186225e-06 4.825469 18.044477 18.044477 0.10740760
## 5 2.312457e-07 4.971811 13.776424 13.776424 0.08200252
## 6 5.901730e-08 6.812750  3.787605  3.787605 0.02254527
##           er        f0         r1         r2           r3
## 1 0.42053470 0.4883753 0.18091287 0.06768109 1.889607e-03
## 2 0.32535966 0.4883753 0.18091287 0.06768109 1.889607e-03
## 3 0.13257164 0.2692079 0.18091287 0.06768109 8.418804e-09
## 4 0.18044477 0.2692079 0.03313816 0.06768109 1.889607e-03
## 5 0.13776424 0.2692079 0.03313816 0.06768109 1.889607e-03
## 6 0.03787605 0.2692079 0.03313816 0.06768109 8.418804e-09
##          f4         r5        jinst
## 1 1.0000000 0.01584893 6.284243e-03
## 2 0.1949426 0.01584893 7.343489e-03
## 3 0.0497522 0.01584893 5.138501e-08
## 4 1.0000000 0.01584893 9.118277e-03
## 5 0.1949426 0.01584893 9.394775e-03
## 6 0.0497522 0.01584893 5.931093e-08

Indeed, the incorporation parameters have been used, and the dummy variables all look correct in the output. We can do a more careful check by viewing the times at which incorporation took place by setting add.incorp.rows = TRUE. Where needed for incorporation calculations, the function will add extra rows which are excluded from the results by default (they can be returned with the results by setting the add.incorp.rows argument to TRUE).

pred6b <- alfam2(dat6, app.name = "TAN.app", time.name = "ctime", 
                 time.incorp = "t.incorp", group = "scenario", 
                 warn = TRUE, add.incorp.rows = TRUE)
## Default parameters (Set 3) are being used.
## Incorporation applied for groups: 2, 3, 5, 6.
## Warning in alfam2(dat6, app.name = "TAN.app", time.name = "ctime", time.incorp = "t.incorp", : Running with 20 parameters. Dropped 2 with no match.
## These secondary parameters have been dropped:
##   rain.rate.r2
##   rain.rate.r5
print(pred6b)
##    scenario ctime  dt            f         s         e
## 1         1   168 168 3.556088e-17  3.325689 42.053470
## 2         2     4   4 1.806765e+01 55.698329 22.801814
## 3         2   168 164 6.932330e-18  3.886253 32.535966
## 4         3     4   4 9.959460e+00 72.505664 12.896523
## 5         3   168 164 9.752566e-19  6.103600 13.257164
## 6         4   168 168 1.186225e-06  4.825469 18.044477
## 7         5     4   4 1.798650e+01 73.849711  3.493427
## 8         5   168 164 2.312457e-07  4.971811 13.776424
## 9         6     4   4 1.798650e+01 73.849711  3.493427
## 10        6   168 164 5.901730e-08  6.812750  3.787605
##            ei           j         er        f0         r1
## 1  42.0534700 0.250318274 0.42053470 0.4883753 0.18091287
## 2  22.8018137 5.700453414 0.22801814 0.4883753 0.18091287
## 3   9.7341520 0.059354585 0.32535966 0.4883753 0.18091287
## 4  12.8965230 3.224130747 0.12896523 0.2692079 0.18091287
## 5   0.3606414 0.002199033 0.13257164 0.2692079 0.18091287
## 6  18.0444771 0.107407602 0.18044477 0.2692079 0.03313816
## 7   3.4934269 0.873356735 0.03493427 0.2692079 0.03313816
## 8  10.2829968 0.062701200 0.13776424 0.2692079 0.03313816
## 9   3.4934269 0.873356735 0.03493427 0.2692079 0.03313816
## 10  0.2941781 0.001793769 0.03787605 0.2692079 0.03313816
##            r2           r3        f4         r5        jinst
## 1  0.06768109 1.889607e-03 1.0000000 0.01584893 6.284243e-03
## 2  0.06768109 1.889607e-03 1.0000000 0.01584893 3.373919e+00
## 3  0.06768109 1.889607e-03 0.1949426 0.01584893 7.343489e-03
## 4  0.06768109 1.889607e-03 1.0000000 0.01584893 1.938802e+00
## 5  0.06768109 8.418804e-09 0.0497522 0.01584893 5.138501e-08
## 6  0.06768109 1.889607e-03 1.0000000 0.01584893 9.118277e-03
## 7  0.06768109 1.889607e-03 1.0000000 0.01584893 7.355865e-01
## 8  0.06768109 1.889607e-03 0.1949426 0.01584893 9.394775e-03
## 9  0.06768109 1.889607e-03 1.0000000 0.01584893 7.355865e-01
## 10 0.06768109 8.418804e-09 0.0497522 0.01584893 5.931093e-08

Here we see the expected effects in primary parameters r3 and f4 at the incorporation times (rows 2, 4, 7 and 9). But note that the f4 mass transfer occurs at the beginning of an interval, and therefore its effects show up in the row with 168 h in the output. This row has results for the interval that starts at the time given in the previous row (or 0, as in scenario 1) and ends at 168 h. For the scenarios with incorporation, the parameter values that reflect incorporation effects are present in any intervals that follow incorporation.

Volatilization dynamics

All the calls above returned results for a single time per scenario (with the exception of some incorporation results, but even they were based on a single input row per scenario). The alfam2 function also predicts dynamics. If your interest is final cumulative emission and changes in weather over time are not important, it is not necessary to look at dynamics. The model uses an analytical expression in each interval, and so results are independent of time step size, as long as conditions (e.g., wind or air temperature) are constant. However, if detailed temporal weather data are available, running the model with multiple intervals will generally improve the accuracy of prediction of final cumulative emission.

Let’s assume we have some high resolution measurements of weather conditions with 2 hour intervals. We’ll create some data to represent this below.

set.seed(1201)
dat7 <- data.frame(ctime = 0:84*2, TAN.app = 100, man.dm = 8, 
                   air.temp = 7 + 7*sin(0:84*2 * 2*pi/24) + 
                              rnorm(85, 0, 2), 
                   wind.sqrt = sqrt(1.5 + 0.4*sin(0:84*2 * 2*pi/24)) + 
                              rnorm(85, 0, 0.12), 
                   app.mthd = 'ts')
plot(air.temp ~ ctime, data = dat7, type = 's', col = 'gray45')

plot(wind.sqrt^2 ~ ctime, data = dat7, type = 's', col = 'blue')

Predictions are made as above. By default, multiple rows in dat are assumed to all belog to the same scenario (same plot, same emission trial) (this is the reason the group argument was needed above).

pred7 <- alfam2(dat7, app.name = 'TAN.app', time.name = 'ctime',
                   warn = FALSE)

Cumulative emission and average interval flux are plotted below.

plot(e ~ ctime, data = pred7, type = 'o', xlab = 'Time (h)', 
     ylab = 'Cumulative emission (kg/ha)')

plot(j ~ ctime, data = pred7, type = 'S', col = 'red', 
     xlab = 'Time (h)', ylab = 'Average flux (kg/ha-h)')

Note that instantaneous flux (jinst) will always be lower than interval average flux when flux monotonically decreases over time.

plot(j ~ ctime, data = pred7, type = 'S', col = 'red', 
     xlab = 'Time (h)', ylab = 'Average flux (kg/ha-h)')
points(jinst ~ ctime, data = pred7, col = 'blue')

Dynamics in the case of incorporation may be interesting. As mentioned in the previous section, extra intervals (rows) will be added for incorporation calculations. These rows are excluded from the results by default and can be returned with the results by setting the add.incorp.rows argument to TRUE.

dat8 <- dat7
dat8$incorp <- "deep"
dat8$t.incorp <- 6.5
pred8 <- alfam2(dat8, app.name = 'TAN.app', time.name = 'ctime',
                time.incorp = 't.incorp', warn = FALSE, 
                add.incorp.rows = TRUE)
plot(e ~ ctime, data = pred8, type = 'o', xlab = 'Time (h)', 
     ylab = 'Cumulative emission (kg/ha)')
abline(v = 6.5, col = 'blue', lty = 2)

plot(j ~ ctime, data = pred8, type = 'S', col = 'red', 
     xlab = 'Time (h)', ylab = 'Average flux (kg/ha-h)')
abline(v = 6.5, col = 'blue', lty = 2)

The drop in flux immediately after incorporation is particularly clear in the flux (second) plot.

Faster evaluation

The alfam2 function is not particularly slow, but its speed can become limiting when applied to many individual scenarios or with very high resolution input data. For example, we can create some fake input data for 1000 different agricultural fields. (These examples are not run during vignette building to avoid excessive check times. But readers can copy/paste them into R to see results.)

set.seed(0812)
dat9 <- expand.grid(field = 1:1000, ct = 1:168, 
                    TAN.app = 100, man.dm = 8, 
                    app.rate.ni = 30, man.source = "pig",
                    man.ph = 7, rain.rate = 0,
                    app.mthd = "bsth")

dat9$air.temp <- 7 + 7*sin(dat9$ct * 2 * pi / 24) + 
                 rnorm(1000, 0, 2)
dat9$wind.sqrt <- sqrt(1.5 + 0.4*sin(dat9$ct * 2 * 2 * pi / 24)) + 
                       rnorm(1000, 0, 0.1)
dat9 <- dat9[order(dat9$field, dat9$ct), ]
head(dat9)
dim(dat9)
system.time(
  pred9 <- alfam2(dat9, app.name = 'TAN.app', time.name = 'ct', 
                  group = 'field', warn = FALSE)
)

Let’s take a small subset for plotting.

pred9sub <- subset(pred9, field %in% 1:100)
pred9sub <- pred9sub[order(pred9sub$field), ]
pred9sub[pred9sub$ct == 168, c('er', 'j')] <- NA

The last line, setting some observations to NA just makes the plot clearer.

plot(j ~ ct, data = pred9sub, type = 'S', col = 'red', 
     xlab = 'Time (h)', ylab = 'Average flux (kg/ha-h)')

There are options for improving function speed, but they introduce complexity and increase the possibility of input errors, and so should be avoided unless speed is critical. The speed can be slightly improved by eliminating some checks using the check argument.

system.time(
  alfam2(dat9, app.name = 'TAN.app', time.name = 'ct', 
         group = 'field', check = FALSE, warn = FALSE)
)

But the improvement is generally trivial. With incorporation, the function call tends to be much slower, and there is more room for improvement.

dat9b <- dat9
dat9b$incorp <- 'shallow'
dat9b$t.incorp <- 4
system.time(
  pred9b <- alfam2(dat9b, app.name = 'TAN.app', time.name = 'ct', 
                   time.incorp = "t.incorp", group = 'field', 
                   warn = FALSE)
)

The issue is preparation of incorporation inputs and determination of parameters. This is not a simple process, because incorporation occurs at a particular time, which may not correspond to a row in the input data, and may differ among locations. Large improvements in speed are possible when these steps are skipped, but doing so requires pre-processing of input data. Pre-processing is done with the same alfam2 function, but the value argument must be set to 'incorp' to get the processed data out.

dat9c <- alfam2(dat9b, app.name = 'TAN.app', time.name = 'ct', 
                time.incorp = 't.incorp', group = 'field', 
                warn = FALSE, value = 'incorp')
head(dat9c)

Note the new, somewhat strange, columns that have been added. And note that this pre-processing takes some time. So separating it only saves time if the actual emission predictions are to be made multiple times. Examples include parameter estimation or cases where other inputs change, e.g., simulations aimed at the effect of changing weather or slurry properties. Running the model function to predict emission is now much quicker, because prep.incorp can be set to FALSE. Because we also prepared dummy variables in the call above that created dat9c, we should set prep.dum = FALSE as well.

system.time(
  pred9c <- alfam2(dat9c, app.name = 'TAN.app', time.name = 'ct', 
                   time.incorp = "t.incorp", group = 'field', 
                   warn = FALSE, prep.dum = FALSE, prep.incorp = FALSE, 
                   check = FALSE)
)
head(pred9b)
head(pred9c)
all.equal(pred9b$e, pred9c$e)

Confidence intervals

The alfam2() function can generate confidence intervals if given multiple parameter sets that themselves represent the distribution of possible parameter values. The package now includes a collection of such parameters that is associated with the default parameter set 3.

set.seed(2609)
dat10 <- data.frame(ctime = 0:84*2, TAN.app = 100, man.dm = 8, 
                   air.temp = 7 + 7*sin(0:84*2 * 2*pi/24) + 
                              rnorm(85, 0, 2), 
                   wind.sqrt = sqrt(1.5 + 0.4*sin(0:84*2 * 2*pi/24)) + 
                              rnorm(85, 0, 0.12), 
                   app.mthd = 'bsth')

Here is a normal call without confidence intervals (CI).

pred10 <- alfam2(dat10, app.name = 'TAN.app', time.name = 'ctime',
                   warn = FALSE)
head(pred10)
##   app.mthd.ts app.mthd.bc app.mthd.os app.mthd.cs ctime dt
## 1           0           0           0           0     0  0
## 2           0           0           0           0     2  2
## 3           0           0           0           0     4  2
## 4           0           0           0           0     6  2
## 5           0           0           0           0     8  2
## 6           0           0           0           0    10  2
##          f        s         e       ei         j         er
## 1 81.03334 18.96666  0.000000 0.000000       NaN 0.00000000
## 2 68.48965 28.21226  2.544995 2.544995 1.2724973 0.02544995
## 3 57.06251 35.53978  5.629435 3.084441 1.5422204 0.05629435
## 4 46.56262 41.15166  9.297690 3.668255 1.8341275 0.09297690
## 5 38.40320 45.33282 11.901947 2.604257 1.3021285 0.11901947
## 6 32.34066 48.42813 13.380717 1.478769 0.7393847 0.13380717
##          f0         r1         r2          r3 f4         r5
## 1 0.8103334 0.01285293 0.06768109 0.002049757  1 0.01584893
## 2 0.8103334 0.01640793 0.06768109 0.002049757  1 0.01584893
## 3 0.8103334 0.02358655 0.06768109 0.002049757  1 0.01584893
## 4 0.8103334 0.03399360 0.06768109 0.002049757  1 0.01584893
## 5 0.8103334 0.02864749 0.06768109 0.002049757  1 0.01584893
## 6 0.8103334 0.01822670 0.06768109 0.002049757  1 0.01584893
##       jinst
## 1 1.0803932
## 2 1.1816016
## 3 1.4187555
## 4 1.6671817
## 5 1.1930767
## 6 0.6887293

Confidence intervals (CI) can be added by specifying a confidence level for the conf.int argument. By default, a set of bootstrap parameters (in alfam2pars03var) are used.

predci1 <- alfam2(dat10, app.name = 'TAN.app', time.name = 'ctime', 
                  warn = FALSE, conf.int = 0.90)
head(predci1)
##    ctime app.mthd.ts app.mthd.bc app.mthd.os app.mthd.cs dt
## 1      0           0           0           0           0  0
## 42     2           0           0           0           0  2
## 53     4           0           0           0           0  2
## 64     6           0           0           0           0  2
## 75     8           0           0           0           0  2
## 2     10           0           0           0           0  2
##           f        s         e       ei         j         er
## 1  81.03334 18.96666  0.000000 0.000000       NaN 0.00000000
## 42 68.48965 28.21226  2.544995 2.544995 1.2724973 0.02544995
## 53 57.06251 35.53978  5.629435 3.084441 1.5422204 0.05629435
## 64 46.56262 41.15166  9.297690 3.668255 1.8341275 0.09297690
## 75 38.40320 45.33282 11.901947 2.604257 1.3021285 0.11901947
## 2  32.34066 48.42813 13.380717 1.478769 0.7393847 0.13380717
##           f0         r1         r2          r3 f4         r5
## 1  0.8103334 0.01285293 0.06768109 0.002049757  1 0.01584893
## 42 0.8103334 0.01640793 0.06768109 0.002049757  1 0.01584893
## 53 0.8103334 0.02358655 0.06768109 0.002049757  1 0.01584893
## 64 0.8103334 0.03399360 0.06768109 0.002049757  1 0.01584893
## 75 0.8103334 0.02864749 0.06768109 0.002049757  1 0.01584893
## 2  0.8103334 0.01822670 0.06768109 0.002049757  1 0.01584893
##        jinst        er.lwr       er.upr
## 1  1.0803932 -8.103778e-18 8.481607e-18
## 42 1.1816016  6.625548e-03 3.997716e-02
## 53 1.4187555  1.775491e-02 7.404345e-02
## 64 1.6671817  3.330055e-02 1.235502e-01
## 75 1.1930767  4.584154e-02 1.585629e-01
## 2  0.6887293  5.262917e-02 1.779778e-01

Confidence limits are given in the output with .lwr and .upr suffixes (for lower and upper). By default CI are only returned for variable er = relative cumulative emission.

plot(er ~ ctime, data = predci1, type = 'l', 
     ylim = c(0, max(predci1$er.upr)))
lines(er.lwr ~ ctime, data = predci1, type = 'l', col = 'blue')
lines(er.upr ~ ctime, data = predci1, type = 'l', col = 'red')

These results come with some warnings. Only uncertainty based on variability in measurements is included; error in model structure or effects of the specific approach used for parameter estimation are not considered. Still, the 90% CI generated with default bootstrap parameter sets is quite wide. Furthermore, estimates are incomplete for any scenarios that include incorporation or closed slot injection, since the parameter estimation observations with these came from a small number of institutions.

We can add any output variables for CI calculation, but internally the quantile function is applied by variable, so the limits returned from different variables may be from different parameter sets and so should be considered separately. (Alternatively, all results from the pars.ci can be returned for external data processing by setting conf.int = 'all'.) Use the var.ci argument to specify variables. Here we request 3 variables.

predci2 <- alfam2(dat10, app.name = 'TAN.app', time.name = 'ctime', 
                  warn = FALSE, conf.int = 0.90, var.ci = c('er', 'j', 'r1'))
head(predci2)
##   ctime app.mthd.ts app.mthd.bc app.mthd.os app.mthd.cs dt
## 1     0           0           0           0           0  0
## 2     2           0           0           0           0  2
## 3     4           0           0           0           0  2
## 4     6           0           0           0           0  2
## 5     8           0           0           0           0  2
## 6    10           0           0           0           0  2
##          f        s         e       ei         j         er
## 1 81.03334 18.96666  0.000000 0.000000       NaN 0.00000000
## 2 68.48965 28.21226  2.544995 2.544995 1.2724973 0.02544995
## 3 57.06251 35.53978  5.629435 3.084441 1.5422204 0.05629435
## 4 46.56262 41.15166  9.297690 3.668255 1.8341275 0.09297690
## 5 38.40320 45.33282 11.901947 2.604257 1.3021285 0.11901947
## 6 32.34066 48.42813 13.380717 1.478769 0.7393847 0.13380717
##          f0         r1         r2          r3 f4         r5
## 1 0.8103334 0.01285293 0.06768109 0.002049757  1 0.01584893
## 2 0.8103334 0.01640793 0.06768109 0.002049757  1 0.01584893
## 3 0.8103334 0.02358655 0.06768109 0.002049757  1 0.01584893
## 4 0.8103334 0.03399360 0.06768109 0.002049757  1 0.01584893
## 5 0.8103334 0.02864749 0.06768109 0.002049757  1 0.01584893
## 6 0.8103334 0.01822670 0.06768109 0.002049757  1 0.01584893
##       jinst      er.lwr     j.lwr      r1.lwr     er.upr
## 1 1.0803932          NA        NA          NA         NA
## 2 1.1816016 0.006625548 0.3312774 0.003886262 0.03997716
## 3 1.4187555 0.017754913 0.5214475 0.006975498 0.07404345
## 4 1.6671817 0.033300553 0.8385077 0.012949531 0.12355023
## 5 1.1930767 0.045841544 0.6216639 0.008565233 0.15856293
## 6 0.6887293 0.052629166 0.3134134 0.004518824 0.17797783
##       j.upr     r1.upr
## 1        NA         NA
## 2 1.9988579 0.02373994
## 3 1.9753496 0.03100817
## 4 2.4640142 0.04595592
## 5 1.6785980 0.03772437
## 6 0.9750093 0.02411052

Note that times with any NaN etc. in one of var.ci columns will be dropped before applying the quantile() function. So here all lwr and upr limits are NA for time = 0 h because j is -Inf.

Confidence intervals can be applied to different groups as well. Here is a demonstration where dry matter varies between groups.

dat11 <- data.frame(ctime = 168, TAN.app = 50, 
                    app.mthd = 'bc', 
                    man.dm = 1:10, air.temp = 20, wind.sqrt = 2) 
predci3 <- alfam2(dat11, app.name = 'TAN.app', time.name = 'ctime', 
                  group = 'man.dm', conf.int = 0.90)
## Default parameters (Set 3) are being used.
## Warning in alfam2(dat = dat, pars = pars, add.pars = add.pars, app.name = app.name, : Incorporation columns  were dropped 
##     because argument time.incorp is NULL
##     So there is no incorporation.
##     Set check = FALSE to not drop, but then check output.
## Warning in alfam2(dat = dat, pars = pars, add.pars = add.pars, app.name = app.name, : Running with 14 parameters. Dropped 8 with no match.
## These secondary parameters have been dropped:
##   man.source.pig.f0
##   man.ph.r1
##   rain.rate.r2
##   incorp.deep.r3
##   man.ph.r3
##   incorp.shallow.f4
##   incorp.deep.f4
##   rain.rate.r5
print(predci3)
##    man.dm ctime app.mthd.ts app.mthd.bc app.mthd.os app.mthd.cs
## 5      -5   168           0           1           0           0
## 4      -4   168           0           1           0           0
## 3      -3   168           0           1           0           0
## 2      -2   168           0           1           0           0
## 1      -1   168           0           1           0           0
## 6       0   168           0           1           0           0
## 7       1   168           0           1           0           0
## 8       2   168           0           1           0           0
## 9       3   168           0           1           0           0
## 10      4   168           0           1           0           0
##     dt            f         s        e       ei          j
## 5  168 2.545770e-56 2.2141944 10.10095 10.10095 0.06012469
## 4  168 2.262953e-52 2.0793418 12.53965 12.53965 0.07464078
## 3  168 1.028902e-48 1.8955697 15.86507 15.86507 0.09443494
## 2  168 2.477581e-45 1.6658229 20.02572 20.02572 0.11920070
## 1  168 3.266919e-42 1.4081777 24.69664 24.69664 0.14700378
## 6  168 2.444469e-39 1.1531110 29.32805 29.32805 0.17457173
## 7  168 1.080170e-36 0.9313427 33.36468 33.36468 0.19859930
## 8  168 2.942948e-34 0.7612382 36.47372 36.47372 0.21710548
## 9  168 5.162215e-32 0.6454527 38.60617 38.60617 0.22979866
## 10 168 6.070504e-30 0.5761587 39.90290 39.90290 0.23751725
##           er        f0        r1         r2          r3 f4
## 5  0.2020190 0.1145835 0.7046755 0.06768109 0.002049757  1
## 4  0.2507930 0.1757817 0.6531002 0.06768109 0.002049757  1
## 3  0.3173014 0.2600650 0.6052998 0.06768109 0.002049757  1
## 2  0.4005143 0.3667769 0.5609979 0.06768109 0.002049757  1
## 1  0.4939327 0.4883753 0.5199384 0.06768109 0.002049757  1
## 6  0.5865610 0.6113652 0.4818841 0.06768109 0.002049757  1
## 7  0.6672936 0.7216410 0.4466150 0.06768109 0.002049757  1
## 8  0.7294744 0.8103334 0.4139272 0.06768109 0.002049757  1
## 9  0.7721235 0.8756362 0.3836319 0.06768109 0.002049757  1
## 10 0.7980580 0.9206566 0.3555538 0.06768109 0.002049757  1
##            r5       jinst    er.lwr    er.upr
## 5  0.01584893 0.004538560 0.1277641 0.2825293
## 4  0.01584893 0.004262145 0.1472534 0.3432488
## 3  0.01584893 0.003885457 0.1933492 0.4139638
## 2  0.01584893 0.003414532 0.2563515 0.4901894
## 1  0.01584893 0.002886422 0.3492118 0.5733555
## 6  0.01584893 0.002363597 0.4626098 0.6416027
## 7  0.01584893 0.001909026 0.5451671 0.7028040
## 8  0.01584893 0.001560353 0.5997174 0.7597271
## 9  0.01584893 0.001323021 0.6080563 0.7992494
## 10 0.01584893 0.001180985 0.6007442 0.8283200
plot(dat11$man.dm, predci3$er, type = 'o', 
     ylim = c(0, max(predci3$er.upr)))
lines(dat11$man.dm, predci3$er.lwr, col = 'blue')
lines(dat11$man.dm, predci3$er.upr, col = 'blue')

By default the model is run with all the different parameter sets provided in pars.ci, which is 100 in the default bootstrap parameter object.

dim(alfam2pars03var)
## [1] 100  22

It is possible to reduce the number used with the n.ci argument.

It may be convenient to combine uncertainty in parameter estimates with uncertainty in model inputs. For this, the dat input data frame must have values that reflect input variable uncertainty. For example, suppose we know that in a particular slurry application event dry matter uncertainty was 2% of fresh mass as a standard deviation.

datuc1 <- data.frame(group = 1:100, ctime = 168, TAN.app = 50, 
                     app.mthd = 'bc', 
                     man.dm = rnorm(100, mean = 8, sd = 2), 
                     air.temp = 20, wind.sqrt = 2)
quantile(datuc1$man.dm)
##        0%       25%       50%       75%      100% 
##  4.361776  6.664609  7.978588  9.508471 13.134127

In this case, we want to get results for all bootstrap parameter sets, not just the quantiles for the confidence intervals. We can get these by setting conf.int = 'all'.

preduc1 <- alfam2(datuc1, app.name = 'TAN.app', time.name = 'ctime', 
                  group = 'group', conf.int = 'all')
## Default parameters (Set 3) are being used.
## Warning in alfam2(dat = dat, pars = pars, add.pars = add.pars, app.name = app.name, : Incorporation columns  were dropped 
##     because argument time.incorp is NULL
##     So there is no incorporation.
##     Set check = FALSE to not drop, but then check output.
## Warning in alfam2(dat = dat, pars = pars, add.pars = add.pars, app.name = app.name, : Running with 14 parameters. Dropped 8 with no match.
## These secondary parameters have been dropped:
##   man.source.pig.f0
##   man.ph.r1
##   rain.rate.r2
##   incorp.deep.r3
##   man.ph.r3
##   incorp.shallow.f4
##   incorp.deep.f4
##   rain.rate.r5
head(preduc1)
##    group ctime app.mthd.ts app.mthd.bc app.mthd.os app.mthd.cs
## 1      1   168           0           1           0           0
## 89     1   168           0           1           0           0
## 96     1   168           0           1           0           0
## 92     1   168           0           1           0           0
## 37     1   168           0           1           0           0
## 16     1   168           0           1           0           0
##     dt            f         s        e       ei         j
## 1  168 5.092709e-35 1.0130663 32.42199 32.42199 0.1929880
## 89 168 2.187692e-10 1.2915088 29.77168 29.77168 0.1772124
## 96 168 1.729811e-14 1.1001513 32.91509 32.91509 0.1959231
## 92 168 3.527354e-31 0.7801358 36.49492 36.49492 0.2172317
## 37 168 6.076613e-34 0.7450051 36.80208 36.80208 0.2190600
## 16 168 3.798794e-13 1.1387640 30.90494 30.90494 0.1839580
##           er        f0         r1         r2          r3 f4
## 1  0.6484397 0.6465281 0.46448352 0.02622215 0.001744226  1
## 89 0.5954337 0.9261497 0.09299761 0.06223041 0.001697912  1
## 96 0.6583018 0.9546963 0.13969538 0.07193463 0.001393751  1
## 92 0.7298985 0.8232866 0.37384233 0.06566458 0.001844284  1
## 37 0.7360417 0.7959515 0.42240658 0.05477944 0.002002155  1
## 16 0.6180987 0.9545427 0.11511617 0.07812444 0.002182948  1
##            r5       jinst par.id
## 1  0.01584893 0.001767016      1
## 89 0.01584893 0.002192868      2
## 96 0.01584893 0.001533337      3
## 92 0.01584893 0.001438792      4
## 37 0.01584893 0.001491616      5
## 16 0.01584893 0.002485863      6
dim(preduc1)
## [1] 10000    21

The output has 10000 rows. And then, for a 90% confidence interval that includes uncertainty in both inputs (only dry matter here) and parameters, we can use quatile():

quantile(preduc1$er, c(0.05, 0.95))
##        5%       95% 
## 0.4517487 0.8148387

Compare this to the confidence interval based on parameter uncertainty only.

datuc2 <- data.frame(group = 1, ctime = 168, TAN.app = 50, 
                     app.mthd = 'bc', 
                     man.dm = 8,
                     air.temp = 20, wind.sqrt = 2)
preduc2 <- alfam2(datuc2, app.name = 'TAN.app', 
                  time.name = 'ctime', group = 'group',
                  conf.int = 0.9)
## Default parameters (Set 3) are being used.
## Warning in alfam2(dat = dat, pars = pars, add.pars = add.pars, app.name = app.name, : Incorporation columns  were dropped 
##     because argument time.incorp is NULL
##     So there is no incorporation.
##     Set check = FALSE to not drop, but then check output.
## Warning in alfam2(dat = dat, pars = pars, add.pars = add.pars, app.name = app.name, : Running with 14 parameters. Dropped 8 with no match.
## These secondary parameters have been dropped:
##   man.source.pig.f0
##   man.ph.r1
##   rain.rate.r2
##   incorp.deep.r3
##   man.ph.r3
##   incorp.shallow.f4
##   incorp.deep.f4
##   rain.rate.r5
print(preduc2)
##   group ctime app.mthd.ts app.mthd.bc app.mthd.os app.mthd.cs
## 1     1   168           0           1           0           0
##    dt            f         s        e       ei         j
## 1 168 2.942948e-34 0.7612382 36.47372 36.47372 0.2171055
##          er        f0        r1         r2          r3 f4
## 1 0.7294744 0.8103334 0.4139272 0.06768109 0.002049757  1
##           r5       jinst    er.lwr    er.upr
## 1 0.01584893 0.001560353 0.5997174 0.7597271

Data import and export

The material in this section is for new R users. Any of the results shown above can be exported as with any data frame in R. The simplest function for this is write.csv(). The CRAN manuals mentioned above include one on import and export: https://cran.r-project.org/ (select “Manuals” at the lower left). The following call will create a comma delimited text file that can be opened with spreadsheet or text editor programs.

write.csv(pred7, 'pred7.csv', row.names = FALSE)

Alternatives include write.csv2, write.table, and for users of the data.table package, fwrite, along with any of the various functions in add-on packages for writing to Excel files.

Except for simple scenarios, it is not very efficient to create a data frame for entering predictor variable values. A more typical approach will be to read data into R from a file, especially when using the model in association with emission measurements. The simplest approach here is to use the read.csv(), fread() or some of the related functions. Alternatively, data can be easily read from Excel files with the read\_xls and related functions in the readxl package.

Note that the alfam2() function can accept data.tables and tibbles. But output is always a data frame, which can of course be changed to a data.table or a tibble after the alfam2 call. You can see this by running the code below. It is not run here in the vignette because not all users will have the necessary packages installed.

library(data.table)
library(ALFAM2)
dat1b <- data.table(ctime = 168, TAN.app = 50, man.dm = 8, 
                   air.temp = 20, wind.sqrt = 2, 
                   app.mthd = 'bc')
dat1b
pred1b <- alfam2(dat1b, app.name = 'TAN.app', time.name = 'ctime')
pred1b
class(pred1b)
setDT(pred1b)
class(pred1b)
library(tibble)
dat1c <- tibble(ctime = 168, TAN.app = 50, man.dm = 8, 
                   air.temp = 20, wind.sqrt = 2, 
                   app.mthd.bc = TRUE)
dat1c
class(dat1c)
pred1c <- alfam2(dat1c, app.name = 'TAN.app', time.name = 'ctime')
class(pred1c)
pred1c <- as_tibble(pred1c)
class(pred1c)

Warnings and error messages

The alfam2() function will return errors or warnings for some cases of problems with input data. The intent of these is to make it easy for users to understand and fix problems with a call or the predictor variable data frame. A few are shown here (but they are not evaluated in order to avoid problems during package checking, so copy and paste the code to see the actual messages). Some checks can be skipped by setting check = FALSE, but doing so increases the risk of getting an unhelpful error message, or worse, overlooking a mistake and proceeding with incorrect output. Some warnings can be surpressed with warn = FALSE in order to avoid cluttering up the console or a log or dynamic report. But users should be careful here also.

The calls above already demonstrate the information on missing predictor variables (which can be supressed by setting warn = FALSE).

Missing values in predictors:

dat6b <- dat6
dat6b[3, 'wind.sqrt'] <- NA

alfam2(dat6b, app.name = "TAN.app", time.name = "ctime", 
       time.incorp = "t.incorp", group = "scenario", 
       check = TRUE, warn = FALSE)
##    app.mthd.os    app.mthd.cs man.source.pig         man.dm 
##              0              0              0              0 
##    app.mthd.bc    app.mthd.ts         man.ph       air.temp 
##              0              0              0              0 
##      wind.sqrt    incorp.deep incorp.shallow 
##              1              0              0
## Error in alfam2(dat6b, app.name = "TAN.app", time.name = "ctime", time.incorp = "t.incorp", : Missing value(s) in predictor variable(s)
##    See above for variables.
##    Check these rows: 3

The printed messages identify which variables and which rows are missing.

Using the wrong names when referring to a column:

pred0 <- alfam2(dat6, app.name = "TAN.app", time.name = "ctim", 
       time.incorp = "t.incorp", group = "scenario")
## Error in alfam2(dat6, app.name = "TAN.app", time.name = "ctim", time.incorp = "t.incorp", : time.name argument you specified (ctim) is not present in dat data frame, which has these columns: scenario, ctime, TAN.app, man.dm, man.ph, air.temp, wind.sqrt, man.source, app.mthd, incorp, t.incorp

Negative time in input data:

dat6c <- dat6
dat6c[1, 'ctime'] <- -10

pred0 <- alfam2(dat6c, app.name = "TAN.app", time.name = "ctime", 
       time.incorp = "t.incorp", group = "scenario", warn = TRUE)
## Warning in alfam2(dat6c, app.name = "TAN.app", time.name =
## "ctime", time.incorp = "t.incorp", : Negative times (variable
## "ctime") found and set to 0.
## Default parameters (Set 3) are being used.
## Incorporation applied for groups: 2, 3, 5, 6.
## Warning in alfam2(dat6c, app.name = "TAN.app", time.name = "ctime", time.incorp = "t.incorp", : Running with 20 parameters. Dropped 2 with no match.
## These secondary parameters have been dropped:
##   rain.rate.r2
##   rain.rate.r5

Trying to use dummy variable prep when it isn’t needed:

dat7 <- data.frame(ctime = 1:10, TAN.app = 100)
pred0 <- alfam2(dat7, app.name = "TAN.app", time.name = "ctime")
## Default parameters (Set 3) are being used.
## Warning in prepDat(dat, warn = warn): Argument prep.dum = TRUE but there are no variables to convert to dummy variables!
##   Ignoring prep.dum = TRUE.
## Warning in alfam2(dat7, app.name = "TAN.app", time.name = "ctime"): Incorporation columns  were dropped 
##     because argument time.incorp is NULL
##     So there is no incorporation.
##     Set check = FALSE to not drop, but then check output.
## Warning in alfam2(dat7, app.name = "TAN.app", time.name = "ctime"): Running with 5 parameters. Dropped 17 with no match.
## These secondary parameters have been dropped:
##   app.mthd.os.f0
##   app.mthd.cs.f0
##   man.source.pig.f0
##   man.dm.f0
##   app.mthd.bc.r1
##   app.mthd.ts.r1
##   man.dm.r1
##   man.ph.r1
##   air.temp.r1
##   wind.sqrt.r1
##   rain.rate.r2
##   app.mthd.cs.r3
##   incorp.deep.r3
##   man.ph.r3
##   incorp.shallow.f4
##   incorp.deep.f4
##   rain.rate.r5

Using an external parameter set (better know what you are doing, not like the example below!):

pred0 <- alfam2(dat6, pars = c(int.f0 = -1, int.r1 = 0), app.name = "TAN.app", 
                time.name = "ctime", time.incorp = "t.incorp", 
                group = "scenario")
## User-supplied parameters are being used.
## Warning in prepIncorp(dat, pars, time.name, time.incorp,
## incorp.names, warn): No incorporation parameters have been
## provided. Skipping incorporation.
## Warning in alfam2(dat6, pars = c(int.f0 = -1, int.r1 = 0), app.name = "TAN.app", : Incorporation columns incorp, t.incorp, incorp.shallow, incorp.deep were dropped 
##     because argument time.incorp is NULL
##     So there is no incorporation.
##     Set check = FALSE to not drop, but then check output.

Note the additional warning about incorporation.

With check = TRUE the alfam2 function checks that all arguments are the right type, and in some cases, contain acceptable values. For example,

alfam2(dat6, app.name = "TAN.app", time.name = "ctim", 
       time.incorp = TRUE, group = "scenario")
## Error: Expect class "character, numeric, integer, NULL" for argument time.incorp but got "logical".
alfam2(dat6, app.name = 95, time.name = "ctim", 
       time.incorp = TRUE, group = "scenario")
## Error: Expect class "character, NULL" for argument app.name but got "numeric".
alfam2(c(ct = 10, man.dm = 2), app.name = 'TAN.app', time.name = "ctim", 
       time.incorp = 6, group = "scenario")
## Error: Expect class "data.frame" for argument dat but got "numeric".

Resulting messages are helpful for identifying problems. Without these checks, the function would typically still return an error, but with a nearly useless message (or worse, incorrect output). For example,

alfam2(dat6, app.name = "TAN.app", time.name = "ctim", 
       time.incorp = TRUE, group = "scenario", check = FALSE)
## Warning in alfam2(dat6, app.name = "TAN.app", time.name =
## "ctim", time.incorp = TRUE, : You set check = FALSE. Be sure to
## verify output!
## Error in `[.data.frame`(dat, , time.name): undefined columns selected

Other examples include the following.

Unavailable options.

pred0 <- alfam2(dat6, app.name = "TAN.app", time.name = "ctim", 
                time.incorp = "t.incorp", group = "scenario", 
                value = "something else")
## Error: Expect one of the following values "emis, incorp" for argument value but got "something else".

Empty input data frame.

datnull <- data.frame()
pred0 <- alfam2(datnull, app.name = "TAN.app", time.name = "ctim", 
                time.incorp = "t.incorp", group = "scenario") 
## Error in alfam2(datnull, app.name = "TAN.app", time.name = "ctim", time.incorp = "t.incorp", : dat has no rows!

Using reserved names.

dat6c <- dat6
dat6c$"__r1" <- 0
pred0 <- alfam2(dat6c, app.name = "TAN.app", 
                time.name = "ctime", time.incorp = "t.incorp", 
                group = "scenario")
## Default parameters (Set 3) are being used.
## Warning in alfam2(dat6c, app.name = "TAN.app", time.name = "ctime", time.incorp = "t.incorp", : dat data frame has some columns with reserved names.
## You can proceed, but there may be problems.
## Better to remove/rename the offending columns: __r1
## Incorporation applied for groups: 2, 3, 5, 6.
## Warning in alfam2(dat6c, app.name = "TAN.app", time.name = "ctime", time.incorp = "t.incorp", : Running with 20 parameters. Dropped 2 with no match.
## These secondary parameters have been dropped:
##   rain.rate.r2
##   rain.rate.r5

Acknowledgements

Christoph Haeni, Johanna Maria Pedersen, Frederik Dalby, and Anders Peter Adamsen provided helpful suggestions and identified errors in earlier drafts of this vignette. Thank you!

References

  1. Hafner, S.D., Pacholski, A., Bittman, S., Carozzi, M., Chantigny, M., Génermont, S., Häni, C., Hansen, M.N., Huijsmans, J., Kupper, T., Misselbrook, T., Neftel, A., Nyord, T., Sommer, S.G. A flexible semi-empirical model for estimating ammonia volatilization from field-applied slurry. Atmospheric Environment. Atmospheric Environment, 199:474-484, 2018. https://doi.org/10.1016/j.atmosenv.2018.11.034↩︎

  2. Hafner S, Pedersen J, Fuss R, Kamp J, Dalby F, Amon B, Pacholski A, Adamsen A, Sommer S., 2024b. Improved tools for estimation of ammonia emission from field-applied animal slurry: refinement of the ALFAM2 model and database. Atmospheric Environment. https://doi.org/10.1016/j.atmosenv.2024.120910↩︎

  3. Hafner S, Pedersen J, Fuss R, Kamp J, Dalby F, Amon B, Pacholski A, Adamsen A, Sommer S., 2024b. Improved tools for estimation of ammonia emission from field-applied animal slurry: refinement of the ALFAM2 model and database. Atmospheric Environment. https://doi.org/10.1016/j.atmosenv.2024.120910↩︎

  4. Hafner, S.D., Nyord, T., Sommer, S.G., Adamsen, A.P.S. 2021. Estimation of Danish emission factors for ammonia from field-applied liquid manure for 1980 to 2019. Danish Centre for Food and Agriculture, Aarhus University, Aarhus, Denmark. Report no. 2021-0251862. https://pure.au.dk/portal/files/223538048/EFreport23092021.pdf↩︎

  5. Hafner, S.D., Kamp, J.N., Pedersen, J. Experimental and model-based comparison of wind tunnel and inverse dispersion model measurement of ammonia emission from field-applied animal slurry. Agricultural and Forest Meteorology, 344, 109790, 2024. https://doi.org/10.1016/j.agrformet.2023.109790↩︎

  6. Hafner, S.D., Pacholski, A., Bittman, S., Carozzi, M., Chantigny, M., Génermont, S., Häni, C., Hansen, M.N., Huijsmans, J., Kupper, T., Misselbrook, T., Neftel, A., Nyord, T., Sommer, S.G. A flexible semi-empirical model for estimating ammonia volatilization from field-applied slurry. Atmospheric Environment. Atmospheric Environment, 199:474-484, 2018. https://doi.org/10.1016/j.atmosenv.2018.11.034↩︎

  7. Hafner, S.D., Pacholski, A., Bittman, S., Carozzi, M., Chantigny, M., Génermont, S., Häni, C., Hansen, M.N., Huijsmans, J., Kupper, T., Misselbrook, T., Neftel, A., Nyord, T., Sommer, S.G. A flexible semi-empirical model for estimating ammonia volatilization from field-applied slurry. Atmospheric Environment. Atmospheric Environment, 199:474-484, 2018. https://doi.org/10.1016/j.atmosenv.2018.11.034↩︎