Fit seasonal wave and continuous ancillary data for trend analysis
Function to prepare data and fit the seawaveQ model.
1 2 3 4
is the concentration data
is the continuous (daily) ancillary data
is a character identifier that names the trend analysis run. It is used to label output files.
are the parameters (water-quality constituents) to analyze (omit the the starting character, for example for sulfate data indicated by P00945, enter "00945").
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.
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.
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.
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.
is a character vector indicating which continuous ancillary variables to include, if none are used for analysis, use iwcav=c("none").
is the column name for the dates, should be the same for both cdat and cavdat
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).
has not been implemented yet but will provide additional model options.
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:
|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
|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|
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.
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.
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).
Aldo V. Vecchia and Karen R. Ryberg
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.
The functions that
fitswavecav calls internally:
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[] # example regression call myfit1[][] # first few lines of observed concentrations head(myfit1[]) # first few lines of predicted concentrations head(myfit1[]) # summary statistics for predicted concentrations head(myfit1[])
Want to suggest features or report bugs for rdrr.io? Use the GitHub issue tracker.