Correlation and Regression with censored data

Computes three parametric correlation coefficients for one X variable and the corresponding R squared for multiple X variables, and a regression equation for censored data.

cencorreg(y.var, cen.var, x.vars, LOG = TRUE, verbose = 2)

Arguments

y.var	The column of y (response variable) values plus detection limits.
cen.var	The column of indicators, where 1 (or `TRUE`) indicates a detection limit in the `y.var` column, and 0 (or `FALSE`) indicates a detected value is in `y.var`.
x.vars	One or more uncensored explanatory variable(s). For multiple variables it must be a data frame of numeric, character and factor variables. See Details
LOG	Indicator of whether to compute the regression in the original y units, or on their logarithms. The default is to use the logarithms (`LOG = TRUE`). To compute in original units, specify the option `LOG = FALSE` (or `LOG = 0`).
verbose	default `verbose=2`, see details.

Value

When x.vars is one variable, likelihood, rescaled likelihood and McFaddens correlation coefficient (R) are printed. When x.vars is a data.frame of more than one variable, likelihood, rescaled likelihood and McFaddens coefficent of determination (R2) are printed.

Model coefficients (intercept and slopes), Chi-Squared statistic and p-value for the test that all slope coefficients equal zero (overall test), and model AIC and BIC are provided.

A Q-Q plot of model residuals with corresponding Shapiro-Francia W and p-value are plotted for evaluation of model distributional assumptions when verbose=2 (the default).

Details

x.vars: If one x variable only, enter its name. If multiple x variables, enter the name of a data frame of columns of the x variables. Only columns used as X variables in the regression are allowed. Create this by x.frame <- data.frame (Temp, Flow, Time) for 3 variables (temperature, flow and time) used as the X variables in the regression.

AIC and BIC are printed to help evaluate the ‘best’ regression model.

The default is that the Y variable will be log transformed.

verbose option. Default is 2 which provides full output in the console and qqplots in a graphics window. A value of 1 only provides partial results in the console and no qqplots. A value of 0 provides no output; the returning computations will be stored in the specified object.

References

Helsel, D.R., 2011. Statistics for censored environmental data using Minitab and R, 2nd ed. John Wiley & Sons, USA, N.J.

Helsel, D.R., 2005. Nondetects and Data Analysis: Statistics for Censored Environmental Data, 1st ed. John Wiley and Sons, USA, N.J.

Examples


data(Brumbaugh)

# One variable
cencorreg(Brumbaugh$Hg,Brumbaugh$HgCen,Brumbaugh$SedMeHg)
#>  Likelihood R = 0.2765                     AIC = 366.2191 
#> Rescaled Likelihood R = 0.2855            BIC = 373.9126 
#>  McFaddens R = 0.1691 
#>  
#> Call:
#> survreg(formula = "log(Brumbaugh$Hg)", data = "Brumbaugh$SedMeHg", 
#>     dist = "gaussian")
#> 
#> Coefficients:
#>         (Intercept) `Brumbaugh$SedMeHg` 
#>          -1.7594846           0.1238447 
#> 
#> Scale= 0.960393 
#> 
#> Loglik(model)= -179.6   Loglik(intercept only)= -184.9
#> 	Chisq= 10.58 on 1 degrees of freedom, p= 0.00114 
#> n= 133 

# More than one variable for demonstration purposes
cencorreg(Brumbaugh$Hg,Brumbaugh$HgCen,Brumbaugh[,c("SedMeHg","PctWetland")])
#>  Likelihood R2 = 0.1371                     AIC = 359.1937 
#>  Rescaled Likelihood R2 = 0.1461            BIC = 369.7851 
#>  McFaddens R2 = 0.05301555 
#>  
#> Call:
#> survreg(formula = "log(Brumbaugh$Hg)", data = "SedMeHg+PctWetland", 
#>     dist = "gaussian")
#> 
#> Coefficients:
#> (Intercept)     SedMeHg  PctWetland 
#> -1.79140776  0.06634913  0.01402344 
#> 
#> Scale= 0.9264126 
#> 
#> Loglik(model)= -175.1   Loglik(intercept only)= -184.9
#> 	Chisq= 19.61 on 2 degrees of freedom, p= 5.53e-05 
#> n= 133