Omega.GLS: Calculate weighting matrix for GLS regression. (WREG)
In USGS-R/WREG: USGS WREG v. 2.02

Description Usage Arguments Details Value Examples

THe Omega.GLS function calculates the weighting matrix required for generalized least-squares regression, without or without uncertainty in the regional skew.

1 2	Omega.GLS(alpha = 0.01, theta = 0.98, independent, X, Y, recordLengths, LP3, MSEGR = NA, TY = 2, peak = TRUE, distMeth = 2)

`alpha`	A number, required only for “GLS” and “GLSskew”. `alpha` is a parameter used in the estimated cross-correlation between site records. See equation 20 in the WREG v. 1.05 manual. The arbitrary, default value is 0.01. The user should fit a different value as needed.
`theta`	A number, required only for “GLS” and “GLSskew”. `theta` is a parameter used in the estimated cross-correlation between site records. See equation 20 in the WREG v. 1.05 manual. The arbitrary, default value is 0.98. The user should fit a different value as needed.
`independent`	A dataframe containing three variables: `StationID` is the numerical identifier (without a leading zero) of each site, `Lat` is the latitude of the site, in decimal degrees, and `Long` is the longitude of the site, in decimal degrees. The sites must be presented in the same order as `Y`. Required only for “GLS” and “GLSskew”.
`X`	The independent variables in the regression, with any transformations already applied. Each row represents a site and each column represents a particular independe variable. (If a leading constant is used, it should be included here as a leading column of ones.) The rows must be in the same order as the dependent variables in `Y`.
`Y`	The dependent variable of interest, with any transformations already applied.
`recordLengths`	This input is required. #' `recordLengths` should be a matrix whose rows and columns are in the same order as `Y`. Each `(r,c)` element represents the length of concurrent record between sites `r` and `c`. The diagonal elements therefore represent each site's full record length.
`LP3`	A dataframe containing the fitted Log-Pearson Type III standard deviate, standard deviation and skew for each site. The names of this data frame are `S`, `K` and `G`. For “GLSskew”, the regional skew value must also be provided in a variable called `GR`. The order of the rows must be the same as `Y`.
`MSEGR`	A number. The mean squared error of the regional skew. Required only for “GLSskew”.
`TY`	A number. The return period of the event being modeled. Required only for “GLSskew”. The default value is `2`. (See the `Legacy` details below.)
`peak`	A logical. Indicates if the event being modeled is a peak flow event or a low-flow event. `TRUE` indicates a peak flow, while `FALSE` indicates a low-flow event.
`distMeth`	Required for “GLS” and “GLSskew”. A value of `1` indicates that the "Nautical Mile" approximation should be used to calculate inter-site distances. A value of `2` designates the Haversine approximation. See `Dist.WREG`. The default value is `2`. (See the `Legacy` details below.)

This function is largely a subroutine for WREG.GLS.

The weighting matrix is calculated by iteration, as noted in the manual to WREG v. 1.0. As currently implemented the initial estimate of model error variance, GSQ, is taken to range from 0 and 2*var(Y). This interval is broken into 30 equally spaced intervals. The weighting matrix is calculated for each interval endpoint and the deviation from equation 21 in the WREG v. 1.0 manual is recorded. The progam then search for the interval over which the deviatioon changes sign. This interval is then split into 30 finer intervals and the process is repeated. The ine interval with the smallest positive deviation is selected as the best estimate.

This function returns a list with two elements:

`GSQ`	The estimated model error variance.
`Omega`	The estimated weighting matrix. A square matrix.

# Import some example data
peakFQdir <- paste0(
  file.path(system.file("exampleDirectory", package = "WREG"),
    "pfqImport"))
gisFilePath <- file.path(peakFQdir, "pfqSiteInfo.txt")
importedData <- importPeakFQ(pfqPath = peakFQdir, gisFile = gisFilePath)

# Organizing input data
lp3Data <- importedData$LP3f
lp3Data$K <- importedData$LP3k$AEP_0.5
Y <- importedData$Y$AEP_0.5
X <- importedData$X[c("Sand", "OutletElev", "Slope")]

# Compute weighting matrix
weightingResult <- Omega.GLS(alpha = 0.01, theta = 0.98,
  independent = importedData$BasChars, X = X,
  Y = Y, recordLengths = importedData$recLen,
  LP3 = lp3Data, MSEGR = NA, TY = 20, peak = TRUE, distMeth = 2)