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.
## 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", ...)
## NULLYou 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
?:
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 | \(^\circ\)C | |
wind.2m |
Wind speed (at 2 m) | m/s | |
wind.sqrt |
Square root of wind speed (2 m) | m\(^{1/2}\)/s\(^{1/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:
## [1] '4.2'and
## 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.
## 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.48425409These numbers indicate a primary parameter. So, for example, the
(secondary) parameter air.temp.r1 is used in the
calculation of the primary parameter \(r_1\). 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\(^\circ\)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 bcOur 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.
## 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.r5The 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.
## 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.001560353The 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## 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## [1] TRUEThe 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.
And add the two new columns. Here we are specifying that deep incorporation happens after 0.5 hours.
## ctime TAN.app man.dm air.temp wind.sqrt app.mthd incorp
## 1 168 50 8 20 2 bc deep
## t.incorp
## 1 0.5pred2 <- 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-08Here we see that with incorporation, emission drops to 18% of applied TAN. Shallow incorporation has less of an effect.
## ctime TAN.app man.dm air.temp wind.sqrt app.mthd incorp
## 1 168 50 8 20 2 bc shallow
## t.incorp
## 1 0.5pred3 <- 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.003799257Relative 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.
## ctime
## 1 168## 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## 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-05Notice 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.
## [1] 0.2954223With 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 NAIt 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 NALet’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-03We can see that predicted emission increases substantially as incorporation time goes up.
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 ospred5 <- 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.004506999Dummy variables
Most of the examples given above include “dummy variables”, or
multiple binary variables that together represent a factor. The package
(\(\geq\)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 TRUEFor 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## [1] TRUECalculated 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:
## incorp.deep.r3 incorp.shallow.f4 incorp.deep.f4
## -5.351121 -1.418209 -2.949668The 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 4pred6 <- 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## 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-08Indeed, 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## 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-08Here 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')
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).
Cumulative emission and average interval flux are plotted below.
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.
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), ]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')] <- NAThe 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.
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.
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).
## 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.6887293Confidence 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-01Confidence 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.02411052Note 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## 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.8283200plot(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.
## [1] 100 22It 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.134127In 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## 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## [1] 10000 21The 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():
## 5% 95%
## 0.4517487 0.8148387Compare 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## 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.7597271Data 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.
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')
dat1bpred1b <- alfam2(dat1b, app.name = 'TAN.app', time.name = 'ctime')
pred1b
class(pred1b)
setDT(pred1b)
class(pred1b)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: 3The 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.incorpNegative 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.r5Trying 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.r5Using 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,
## Error: Expect class "character, numeric, integer, NULL" for argument time.incorp but got "logical".## 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 selectedOther 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