Fit seasonal wave and continuous ancillary data for trend analysis

Share:

Description

Function to prepare data and fit the seawaveQ model.

Usage

1
2
3
4
  fitswavecav(cdat, cavdat, tanm = "trend1", pnames,
    yrstart = 0, yrend = 0, tndbeg = 0, tndend = 0,
    iwcav = c("none"), dcol = "dates",
    qwcols = c("R", "P"), mclass = 1)

Arguments

cdat

is the concentration data

cavdat

is the continuous (daily) ancillary data

tanm

is a character identifier that names the trend analysis run. It is used to label output files.

pnames

are the parameters (water-quality constituents) to analyze (omit the the starting character, for example for sulfate data indicated by P00945, enter "00945").

yrstart

is the starting year of the analysis (treated as January 1 of that year). Zero means the start date will be determined by the start date of cavdat, the continuous ancillary data.

yrend

is the ending year of the analysis (treated as December 31 of that year). Zero means the end date will be determined by the end date of cavdat, the continuous ancillary data.

tndbeg

is the beginning (in whole or decimal years) of the trend period. Zero means the begin date will be the beginning of the concentration data, cdat.

tndend

is the end of the trend (treated as December 31 of that year). Zero means the end date will be the end of the concentration data, cdat.

iwcav

is a character vector indicating which continuous ancillary variables to include, if none are used for analysis, use iwcav=c("none").

dcol

is the column name for the dates, should be the same for both cdat and cavdat

qwcols

is a character vector with the beginning of the column headers for remarks code (default is R), and beginning of column headers for concentration data (default is P for parameter).

mclass

has not been implemented yet but will provide additional model options.

Format

The data frame returned has one row for each parameter analyzed and the number of columns depend on the number of continuous ancillary variables used. The general format is as follows:

pname character Parameter analyzed
mclass numeric Currently a value of 1
jmod numeric The choice of pulse input function, an integer 1--14.
hlife numeric the model half-life in months, an integer, 1 to 4 months
cmaxt numeric the decimal season of maximum concentration
scl numeric the scale factor from the survreg.object
loglik numeric the log-likelihood for the model
cint numeric coefficient for model intercept
cwave numeric coefficient for the seasonal wave
ctnd numeric coefficient for the trend component of model
c[alphanumeric] numeric 0 or more coefficients for the continuous ancillary variables
seint numeric standard error for the intercept
sewave numeric standard error for the seasonal wave
setnd numeric standard error for the trend
se[alphanumeric] numeric 0 or more standard errors for the continuous ancillary variables
pvaltnd numeric the p-value for the trend line

Details

Fits the seawaveQ model (Vecchia and others, 2008) using a seasonal wave and continuous ancillary variables (streamflow anomalies and other continuous variables such as conductivity or sediment) to model water quality.

Value

a pdf file containing plots of the data and modeled concentration, a text file containing a summary of the survival regression call for each model selected, and a list. The first element of the list is a data frame described under format. The second element of the list is the summary of the survival regression call. The third element is the observed concentration data (censored and uncensored). The fourth element is the concentration data predicted by the model. The fifth element provides summary statistics for the predicted concentrations.

Note

The assumed data format is one with columns for water-quality concentration values and a related column for qualification of those values, such as in the case of left-censored values less than a particular value. For example, a water-quality sample was collected and the laboratory analysis indicated that the concentration was less than 0.01 micrograms per liter. The USGS parameter code for simazine is 04035 (U.S. Geological Survey, 2013b). When the data are retrieved through the National Water Information System: Web Interface (http://waterdata.usgs.gov/nwis; U.S. Geological Survey, 2013a), the concentration values are in a column labeled P04035 and the qualification information, or remark codes, are in a column labeled R04035. To use this function, the argument pnames would be the unique identifier for simazine values and qualifications, 04035, and the qwcols argument would be c("R", "P") to indicate that the qualification column starts with an R and the values column starts with a P.
Other users may have data in different format that can be changed to use with this function. For example, a user may have concentration values and qualification codes in one column, such as a column labeled simzaine with the values 0.05, 0.10, <0.01, <0.01, and 0.90. In this case, the less thans and any other qualification codes should be placed in a separate column. The column names for the qualification codes and the concentration values should be the same with the exception of different beginning letters to indicate which column is which. The columns could be named Rsimazine and Psimazine. Then the argument pnames = "simazine" and the argument qwcols = c("R", "P").
Users should exercise caution when their water-quality data have multiple censoring limits and may want to recensor the data to a single censoring level. Censoring and recensoring issues are discussed in the text and Appendix 1 of Ryberg and others (2010).

Author(s)

Aldo V. Vecchia and Karen R. Ryberg

References

Ryberg, K.R., Vecchia, A.V., Martin, J.D., and Gilliom, R.J., 2010, Trends in pesticide concentrations in urban streams in the United States, 1992–2008: U.S. Geological Survey Scientific Investigations Report 2010-5139, 101 p. (Also available at http://pubs.usgs.gov/sir/2010/5139/.)

U.S. Geological Survey, 2013a, National Water Information System: Web Interface, accessed Febaruary 26, 2013, at http://waterdata.usgs.gov.

U.S. Geological Survey, 2013b, Parameter code definition: National Water Information System: Web Interface, accessed Febaruary 26, 2013, at http://nwis.waterdata.usgs.gov/usa/nwis/pmcodes.

Vecchia, A.V., Martin, J.D., and Gilliiom, R.J., 2008, Modeling variability and trends in pesticide concentrations in streams: Journal of the American Water Resources Association, v. 44, no. 5, p. 1308-1324, http://dx.doi.org/10.1111/j.1752-1688.2008.00225.x.

See Also

The functions that fitswavecav calls internally:
prepData and fitMod.

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
data(swData)
modMoRivOmaha<-combineData(qwdat=qwMoRivOmaha, cqwdat=cqwMoRivOmaha)
myfit1 <- fitswavecav(cdat=modMoRivOmaha, cavdat=cqwMoRivOmaha,
tanm="myfit1", pnames=c("04035", "04037", "04041"), yrstart=1995,
yrend=2003, tndbeg=1995, tndend=2003, iwcav=c("flowa30","flowa1"),
dcol="dates", qwcols=c("R","P"))
 ## Not run: 
myfit2 <- fitswavecav(cdat=modMoRivOmaha, cavdat=cqwMoRivOmaha,
tanm="myfit2", pnames=c("04035", "04037", "04041"), yrstart=1995,
yrend=2003, tndbeg=1995, tndend=2003, iwcav=c("seda30","seda1"),
dcol="dates", qwcols=c("R","P"))
myfit3 <- fitswavecav(cdat=modMoRivOmaha, cavdat=cqwMoRivOmaha,
tanm="myfit3", pnames=c("04035", "04037", "04041"), yrstart=1995,
yrend=2003, tndbeg=1995, tndend=2003, iwcav=c("flowa30","flowa1",
"seda30", "seda1"), dcol="dates", qwcols=c("R","P"))

## End(Not run)
# trend model results
myfit1[[1]]
# example regression call
myfit1[[2]][[1]]
# first few lines of observed concentrations
head(myfit1[[3]])
# first few lines of predicted concentrations
head(myfit1[[4]])
# summary statistics for predicted concentrations
head(myfit1[[5]])

Want to suggest features or report bugs for rdrr.io? Use the GitHub issue tracker.