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", ...)
## 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 ?:
?alfam2But 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 will be described in a forthcoming paper. 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 factors2 and a later paper3. And set 1 in the 2019 paper listed below4.
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.1.3'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.1.3) 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.43613933 -2.93492578 -7.80196997
## man.source.pig.f0 man.dm.f0 int.r1
## -0.85171386 0.49659337 -1.46760800
## app.mthd.bc.r1 app.mthd.ts.r1 man.dm.r1
## 0.71991146 -0.09333684 -0.02843126
## man.ph.r1 air.temp.r1 wind.sqrt.r1
## 0.44886708 0.03454900 0.46628989
## int.r2 rain.rate.r2 int.r3
## -1.20493824 0.62051420 -2.71593590
## app.mthd.cs.r3 incorp.deep.r3 man.ph.r3
## -0.34883867 -1.96259695 0.03557064
## incorp.shallow.f4 incorp.deep.f4 int.r5
## -1.37979544 -3.26822034 -1.80000000
## rain.rate.r5
## 0.34944126These 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 5. 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.
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"): 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.
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 5.76975e-33 0.7727432 36.46375 36.46375 0.2170461 0.729275
## f0 r1 r2 r3 f4 r5
## 1 0.8067963 0.4014868 0.06238235 0.001923376 1 0.01584893
## jinst
## 1 0.001486275The 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 Broadcastpred1b <- 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"): 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.r5all.equal(pred1, pred1b)## [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.
dat2 <- dat1And 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.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 2.116209e-34 2.920113 8.303934 8.303934
## j er f0 r1 r2
## 1 0.04942818 0.1660787 0.8067963 0.4014868 0.06238235
## r3 f4 r5 jinst
## 1 2.096366e-05 0.03667766 0.01584893 6.121627e-05Here we see that with incorporation, emission drops to 17% 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.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 1.159961e-33 1.892889 16.61548 16.61548
## j er f0 r1 r2
## 1 0.09890167 0.3323096 0.8067963 0.4014868 0.06238235
## r3 f4 r5 jinst
## 1 0.001923376 0.2010419 0.01584893 0.003640736Multiple 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 2.116209e-34 3.2854505
## 2 0 1 168 168 2.116209e-34 2.5463555
## 3 0 1 168 168 2.116209e-34 1.2097252
## 4 0 1 168 168 2.116209e-34 1.0163195
## 5 0 0 168 168 5.769750e-33 0.7727432
## e ei j er f0 r1
## 1 2.865008 2.865008 0.01705362 0.05730016 0.8067963 0.4014868
## 2 13.828814 13.828814 0.08231437 0.27657629 0.8067963 0.4014868
## 3 32.999739 32.999739 0.19642702 0.65999477 0.8067963 0.4014868
## 4 35.477783 35.477783 0.21117728 0.70955566 0.8067963 0.4014868
## 5 36.463748 36.463748 0.21704612 0.72927495 0.8067963 0.4014868
## r2 r3 f4 r5 jinst
## 1 0.06238235 2.096366e-05 0.03667766 0.01584893 6.887508e-05
## 2 0.06238235 2.096366e-05 0.03667766 0.01584893 5.338094e-05
## 3 0.06238235 2.096366e-05 0.03667766 0.01584893 2.536027e-05
## 4 0.06238235 2.096366e-05 0.03667766 0.01584893 2.130578e-05
## 5 0.06238235 1.923376e-03 1.00000000 0.01584893 1.486275e-03We 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 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 2.368445e-23 0.8921571 34.53148 34.53148 0.20554452
## 2 168 2.960794e-09 1.5370593 25.06230 25.06230 0.14918037
## 3 168 1.249129e-12 2.2470157 10.39953 10.39953 0.06190198
## er f0 r1 r2 r3 f4
## 1 0.6906296 0.8067963 0.26972814 0.06238235 0.001923376 1
## 2 0.5012460 0.8067963 0.07651733 0.06238235 0.001923376 1
## 3 0.2079907 0.1815918 0.11389505 0.06238235 0.001923376 1
## r5 jinst
## 1 0.01584893 0.001715953
## 2 0.01584893 0.002956343
## 3 0.01584893 0.004321855Dummy 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 2.368445e-23 0.8921571 34.53148 34.53148
## 2 2 168 168 2.960794e-09 1.5370593 25.06230 25.06230
## 3 3 168 168 1.249129e-12 2.2470157 10.39953 10.39953
## j er f0 r1 r2
## 1 0.20554452 0.6906296 0.8067963 0.26972814 0.06238235
## 2 0.14918037 0.5012460 0.8067963 0.07651733 0.06238235
## 3 0.06190198 0.2079907 0.1815918 0.11389505 0.06238235
## r3 f4 r5 jinst
## 1 0.001923376 1 0.01584893 0.001715953
## 2 0.001923376 1 0.01584893 0.002956343
## 3 0.001923376 1 0.01584893 0.004321855all.equal(pred5b$e, pred5$e)## [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:
alfam2pars03[grepl('^incorp', names(alfam2pars03))]## incorp.deep.r3 incorp.shallow.f4 incorp.deep.f4
## -1.962597 -1.379795 -3.268220The 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.r5print(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 2.154847e-15 3.367835 41.518626 41.518626 0.24713468
## 2 4.332145e-16 3.970446 31.255779 31.255779 0.18604631
## 3 4.670722e-17 6.103392 13.186915 13.186915 0.07849354
## 4 4.535890e-06 4.836810 18.548833 18.548833 0.11040972
## 5 9.119037e-07 4.986344 13.869975 13.869975 0.08255937
## 6 1.663658e-07 6.801341 3.838452 3.838452 0.02284793
## er f0 r1 r2 r3
## 1 0.41518626 0.4848911 0.16173905 0.06238235 1.876692e-03
## 2 0.31255779 0.4848911 0.16173905 0.06238235 1.876692e-03
## 3 0.13186915 0.2865564 0.16173905 0.06238235 2.045483e-05
## 4 0.18548833 0.2865564 0.03082502 0.06238235 1.876692e-03
## 5 0.13869975 0.2865564 0.03082502 0.06238235 1.876692e-03
## 6 0.03838452 0.2865564 0.03082502 0.06238235 2.045483e-05
## f4 r5 jinst
## 1 1.00000000 0.01584893 0.0063203882
## 2 0.20104186 0.01584893 0.0074513020
## 3 0.03667766 0.01584893 0.0001248439
## 4 1.00000000 0.01584893 0.0090773404
## 5 0.20104186 0.01584893 0.0093578576
## 6 0.03667766 0.01584893 0.0001391254Indeed, 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.r5print(pred6b)## scenario ctime dt f s e
## 1 1 168 168 2.154847e-15 3.367835 41.518626
## 2 2 4 4 1.978361e+01 55.658102 21.122409
## 3 2 168 164 4.332145e-16 3.970446 31.255779
## 4 3 4 4 1.169153e+01 70.995544 12.779141
## 5 3 168 164 4.670722e-17 6.103392 13.186915
## 6 4 168 168 4.535890e-06 4.836810 18.548833
## 7 5 4 4 1.973751e+01 72.210464 3.489598
## 8 5 168 164 9.119037e-07 4.986344 13.869975
## 9 6 4 4 1.973751e+01 72.210464 3.489598
## 10 6 168 164 1.663658e-07 6.801341 3.838452
## ei j er f0 r1
## 1 41.5186256 0.247134676 0.41518626 0.4848911 0.16173905
## 2 21.1224091 5.280602276 0.21122409 0.4848911 0.16173905
## 3 10.1333704 0.061788844 0.31255779 0.4848911 0.16173905
## 4 12.7791409 3.194785229 0.12779141 0.2865564 0.16173905
## 5 0.4077736 0.002486424 0.13186915 0.2865564 0.16173905
## 6 18.5488327 0.110409718 0.18548833 0.2865564 0.03082502
## 7 3.4895980 0.872399503 0.03489598 0.2865564 0.03082502
## 8 10.3803767 0.063294980 0.13869975 0.2865564 0.03082502
## 9 3.4895980 0.872399503 0.03489598 0.2865564 0.03082502
## 10 0.3488537 0.002127157 0.03838452 0.2865564 0.03082502
## r2 r3 f4 r5 jinst
## 1 0.06238235 1.876692e-03 1.00000000 0.01584893 0.0063203882
## 2 0.06238235 1.876692e-03 1.00000000 0.01584893 3.3042346857
## 3 0.06238235 1.876692e-03 0.20104186 0.01584893 0.0074513020
## 4 0.06238235 1.876692e-03 1.00000000 0.01584893 2.0242138725
## 5 0.06238235 2.045483e-05 0.03667766 0.01584893 0.0001248439
## 6 0.06238235 1.876692e-03 1.00000000 0.01584893 0.0090773404
## 7 0.06238235 1.876692e-03 1.00000000 0.01584893 0.7439259695
## 8 0.06238235 1.876692e-03 0.20104186 0.01584893 0.0093578576
## 9 0.06238235 1.876692e-03 1.00000000 0.01584893 0.7439259695
## 10 0.06238235 2.045483e-05 0.03667766 0.01584893 0.0001391254Here 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.5pred8 <- 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')] <- 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.
dat9b <- dat9
dat9b$incorp <- 'shallow'
dat9b$t.incorp <- 4system.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
## 1 80.67963 19.32037 7.689719e-16 7.689719e-16 Inf
## 2 68.97939 27.79484 2.474397e+00 2.474397e+00 1.2371987
## 3 58.13839 34.59201 5.525259e+00 3.050861e+00 1.5254307
## 4 47.97373 39.86368 9.234365e+00 3.709106e+00 1.8545532
## 5 40.02169 43.84697 1.187339e+01 2.639022e+00 1.3195110
## 6 34.09765 46.84651 1.335826e+01 1.484874e+00 0.7424372
## er f0 r1 r2 r3 f4
## 1 7.689719e-18 0.8067963 0.01234782 0.06238235 0.001923376 1
## 2 2.474397e-02 0.8067963 0.01595685 0.06238235 0.001923376 1
## 3 5.525259e-02 0.8067963 0.02310844 0.06238235 0.001923376 1
## 4 9.234365e-02 0.8067963 0.03370397 0.06238235 0.001923376 1
## 5 1.187339e-01 0.8067963 0.02823363 0.06238235 0.001923376 1
## 6 1.335826e-01 0.8067963 0.01771422 0.06238235 0.001923376 1
## r5 jinst
## 1 0.01584893 1.0333782
## 2 0.01584893 1.1541535
## 3 0.01584893 1.4100206
## 4 0.01584893 1.6935780
## 5 0.01584893 1.2142916
## 6 0.01584893 0.6941167Confidence 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
## 1 80.67963 19.32037 7.689719e-16 7.689719e-16 Inf
## 42 68.97939 27.79484 2.474397e+00 2.474397e+00 1.2371987
## 53 58.13839 34.59201 5.525259e+00 3.050861e+00 1.5254307
## 64 47.97373 39.86368 9.234365e+00 3.709106e+00 1.8545532
## 75 40.02169 43.84697 1.187339e+01 2.639022e+00 1.3195110
## 2 34.09765 46.84651 1.335826e+01 1.484874e+00 0.7424372
## er f0 r1 r2 r3 f4
## 1 7.689719e-18 0.8067963 0.01234782 0.06238235 0.001923376 1
## 42 2.474397e-02 0.8067963 0.01595685 0.06238235 0.001923376 1
## 53 5.525259e-02 0.8067963 0.02310844 0.06238235 0.001923376 1
## 64 9.234365e-02 0.8067963 0.03370397 0.06238235 0.001923376 1
## 75 1.187339e-01 0.8067963 0.02823363 0.06238235 0.001923376 1
## 2 1.335826e-01 0.8067963 0.01771422 0.06238235 0.001923376 1
## r5 jinst er.lwr er.upr
## 1 0.01584893 1.0333782 -7.623892e-18 7.970631e-18
## 42 0.01584893 1.1541535 6.503375e-03 3.944416e-02
## 53 0.01584893 1.4100206 1.860650e-02 8.146614e-02
## 64 0.01584893 1.6935780 3.901266e-02 1.316395e-01
## 75 0.01584893 1.2142916 5.311177e-02 1.658405e-01
## 2 0.01584893 0.6941167 5.968924e-02 1.845955e-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
## 1 80.67963 19.32037 7.689719e-16 7.689719e-16 Inf
## 2 68.97939 27.79484 2.474397e+00 2.474397e+00 1.2371987
## 3 58.13839 34.59201 5.525259e+00 3.050861e+00 1.5254307
## 4 47.97373 39.86368 9.234365e+00 3.709106e+00 1.8545532
## 5 40.02169 43.84697 1.187339e+01 2.639022e+00 1.3195110
## 6 34.09765 46.84651 1.335826e+01 1.484874e+00 0.7424372
## er f0 r1 r2 r3 f4
## 1 7.689719e-18 0.8067963 0.01234782 0.06238235 0.001923376 1
## 2 2.474397e-02 0.8067963 0.01595685 0.06238235 0.001923376 1
## 3 5.525259e-02 0.8067963 0.02310844 0.06238235 0.001923376 1
## 4 9.234365e-02 0.8067963 0.03370397 0.06238235 0.001923376 1
## 5 1.187339e-01 0.8067963 0.02823363 0.06238235 0.001923376 1
## 6 1.335826e-01 0.8067963 0.01771422 0.06238235 0.001923376 1
## r5 jinst er.lwr j.lwr r1.lwr
## 1 0.01584893 1.0333782 NA NA NA
## 2 0.01584893 1.1541535 0.006503375 0.3251688 0.003883792
## 3 0.01584893 1.4100206 0.018606498 0.5823830 0.007295156
## 4 0.01584893 1.6935780 0.039012658 0.8729245 0.012036472
## 5 0.01584893 1.2142916 0.053111766 0.6494888 0.009616745
## 6 0.01584893 0.6941167 0.059689241 0.3076063 0.004538138
## er.upr j.upr r1.upr
## 1 NA NA NA
## 2 0.03944416 1.972208 0.02502532
## 3 0.08146614 2.184737 0.03507628
## 4 0.13163948 2.516113 0.05096457
## 5 0.16584050 1.710051 0.04281883
## 6 0.18459550 0.967844 0.02699123Note 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, : 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.r5print(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 7.658989e-51 2.2627983 9.80978 9.80978 0.05839155
## 4 168 1.010062e-47 2.1261991 12.24535 12.24535 0.07288899
## 3 168 8.375497e-45 1.9399791 15.56754 15.56754 0.09266391
## 2 168 4.435953e-42 1.7066892 19.73244 19.73244 0.11745498
## 1 168 1.524679e-39 1.4439375 24.42778 24.42778 0.14540343
## 6 168 3.467761e-37 1.1819796 29.11523 29.11523 0.17330495
## 7 168 5.351155e-35 0.9518677 33.24124 33.24124 0.19786453
## 8 168 5.769750e-33 0.7727432 36.46375 36.46375 0.21704612
## 9 168 4.482966e-31 0.6480308 38.72067 38.72067 0.23048017
## 10 168 2.584492e-29 0.5703355 40.14309 40.14309 0.23894696
## er f0 r1 r2 r3 f4
## 5 0.1961956 0.1143733 0.6348777 0.06238235 0.001923376 1
## 4 0.2449070 0.1750525 0.5946465 0.06238235 0.001923376 1
## 3 0.3113507 0.2585266 0.5569646 0.06238235 0.001923376 1
## 2 0.3946487 0.3642309 0.5216705 0.06238235 0.001923376 1
## 1 0.4885555 0.4848911 0.4886130 0.06238235 0.001923376 1
## 6 0.5823046 0.6073387 0.4576503 0.06238235 0.001923376 1
## 7 0.6648248 0.7176294 0.4286497 0.06238235 0.001923376 1
## 8 0.7292750 0.8067963 0.4014868 0.06238235 0.001923376 1
## 9 0.7744134 0.8727971 0.3760452 0.06238235 0.001923376 1
## 10 0.8028618 0.9185280 0.3522157 0.06238235 0.001923376 1
## r5 jinst er.lwr er.upr
## 5 0.01584893 0.004352211 0.1219644 0.3412941
## 4 0.01584893 0.004089479 0.1478557 0.4070369
## 3 0.01584893 0.003731308 0.1793219 0.4768230
## 2 0.01584893 0.003282604 0.2338195 0.5359166
## 1 0.01584893 0.002777234 0.3204419 0.5895909
## 6 0.01584893 0.002273391 0.4173289 0.6421670
## 7 0.01584893 0.001830799 0.5289427 0.7026018
## 8 0.01584893 0.001486275 0.6190506 0.7618915
## 9 0.01584893 0.001246407 0.6313494 0.8054500
## 10 0.01584893 0.001096969 0.6360032 0.8379472plot(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 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, : 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.r5head(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 1.639000e-29 1.0926849 31.60240 31.60240 0.1881095
## 89 168 1.203993e-31 0.7554343 37.03299 37.03299 0.2204345
## 96 168 1.500723e-13 1.0414376 33.26095 33.26095 0.1979819
## 92 168 8.513323e-13 1.0775028 33.14910 33.14910 0.1973161
## 37 168 3.400932e-34 0.8590680 35.16440 35.16440 0.2093119
## 16 168 4.381991e-23 0.8250588 36.00961 36.00961 0.2143429
## er f0 r1 r2 r3 f4
## 1 0.6320479 0.6313041 0.3943238 0.02075317 0.001479486 1
## 89 0.7406599 0.8787759 0.3618556 0.08443783 0.001826523 1
## 96 0.6652191 0.7818514 0.1612706 0.03631033 0.001365639 1
## 92 0.6629820 0.8341612 0.1434809 0.04415406 0.001261752 1
## 37 0.7032879 0.8491639 0.3798781 0.10114791 0.001853033 1
## 16 0.7201922 0.8490724 0.2675120 0.06124022 0.001762766 1
## r5 jinst par.id
## 1 0.01584893 0.001616612 1
## 89 0.01584893 0.001379818 2
## 96 0.01584893 0.001422228 3
## 92 0.01584893 0.001359542 4
## 37 0.01584893 0.001591881 5
## 16 0.01584893 0.001454386 6dim(preduc1)## [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():
quantile(preduc1$er, c(0.05, 0.95))## 5% 95%
## 0.4580098 0.8217002Compare 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, : 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.r5print(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 er
## 1 168 5.76975e-33 0.7727432 36.46375 36.46375 0.2170461 0.729275
## f0 r1 r2 r3 f4 r5
## 1 0.8067963 0.4014868 0.06238235 0.001923376 1 0.01584893
## jinst er.lwr er.upr
## 1 0.001486275 0.6190506 0.7618915Data 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')
dat1bpred1b <- 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: 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"): 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.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" 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