Build an alpha-beta-sigma Matrix for Use with the cplot Function

Description

Creates a 3 \times N (no. of methods) matrix consisting of the estimated alphas, betas, and imprecision sigmas for use with the cplot function.

Usage

alpha.beta.sigma(x)

Arguments

Details

This is primarily a helper function used by the omx function.

Value

A 3 \times N matrix consisting of alphas on the first row, betas on the second row, followed by raw imprecision sigmas.

See Also

cplot, omx.

Examples

## Not run: 
library(OpenMx)
library(merror)
data(pm2.5)
pm <- pm2.5

# OpenMx does not like periods in data column names
names(pm) <- c('ms_conc_1','ws_conc_1','ms_conc_2','ws_conc_2','frm')

# Fit model with FRM sampler as reference
omxfit <- omx(data=pm[,c(5,1:4)],bs.q=c(0.025,0.5,0.975),reps=100)

# Extract the estimates
alpha.beta.sigma(summary(omxfit$fit)$parameters[,c(1,5,6)])

# Make a calibration plot
cplot(pm[,c(5,1:4)],1,2,alpha.beta.sigma=
  alpha.beta.sigma(summary(omxfit$fit)$parameters[,c(1,5,6)]))

# The easier way
cplot(pm[,c(5,1:4)],1,2,alpha.beta.sigma=omxfit$abs)

## End(Not run)

Compute the estimates of betas.

Description

This function is used internally to compute the estimates of betas.

Usage

beta.bar(x)

Arguments

Details

See Jaech, p. 184.

Value

A vector of length N (no. of methods) containing the estimates of beta.

Author(s)

Richard A. Bilonick

References

Jaech, J. L. (1985) Statistical Analysis of Measurement Errors. New York: Wiley.

See Also

cb.pd, ncb.od,lrt

Examples

data(pm2.5)
beta.bar(pm2.5) # estimate betas (accuracy parameter)

Compute accuracy estimates and maximum likelihood estimates of precision for the constant bias measurement error model using paired data.

Description

Compute accuracy estimates and maximum likelihood estimates of precision for the constant bias measurement error model using paired data.

Usage

cb.pd(x, conf.level = 0.95, M = 40)

Arguments

Details

Measurement Error Model:

x[i,k] = alpha[i] + beta[i]*mu[k] + epsilon[i,k]

where x[i,k] is the measurement by the ith method for the kth item, i = 1 to N, k = 1 to n, mu[k] is the true value for the kth item, epsilon[i,k] is the Normally distributed random error with variance sigma[i] squared for the ith method and the kth item, and alpha[i] and beta[i] are the accuracy parameters for the ith method.

The imprecision for the ith method is sigma[i]. If all alphas are zeroes and all betas are ones, there is no bias. If all betas equal 1, then there is a constant bias. Otherwise there is a nonconstant bias.

ME (method of moments estimator) and MLE are the same for N=3 instruments except for a factor of (n-1)/n: MLE = (n-1)/n * ME

Using paired differences forces Constant Bias model (beta[1] = beta[2] = ... = beta[N]). Also, the process variance CANNOT be estimated.

Value

Author(s)

Richard A. Bilonick

References

Jaech, J. L. (1985) Statistical Analysis of Measurement Errors. New York: Wiley.

See Also

ncb.od, lrt

Examples

data(pm2.5)
cb.pd(pm2.5)

Scatter plot of observations for a pair of devices with calibration curve.

Description

Creates a scatter plot for any pair of observations in the data.frame and superimposes the calibration curve.

Usage

cplot(df, i, j, leg.loc="topleft", regress=FALSE, lw=1, t.size=1, alpha.beta.sigma=NULL)

Arguments

Details

By default, cplot displays the corresponding calibration curve for devices i and j based on the parameter estimates for alpha, beta, and sigma computed using ncb.od. You can overide this calibration curve by providing argument alpha.beta.sigma with different estimates. Both naive regression lines (device i regressed on device j, and device j regressed on device i) by setting "regress=TRUE". Note, however, that the calibration curve will fall somewhere between these two regression lines, depending on the the ratio of the imprecision standard deviations (sigmas). (This may not hold if there are missing measurement data values given that ordinary regression requires deleting any item with one or more missing values.)

Value

Produces a scatter plot with the calibration curve and titles that includes the calibration equation and the scale-bias adjusted imprecision standard deviations.

Author(s)

Richard A. Bilonick

References

Jaech, J. L. (1985) Statistical Analysis of Measurement Errors. New York: Wiley.

See Also

merror.pairs

Examples

library(merror)
data(pm2.5)

# Make various calibration plots for pm2.5 measurements
par(mfrow=c(2,2))
cplot(pm2.5,2,1)
cplot(pm2.5,3,1)
cplot(pm2.5,4,1)
# Add the naive regression lines JUST for comparison
cplot(pm2.5,5,1,regress=TRUE,t.size=0.9)

# This is redundant but illustrates using the
# argument alpha.beta.sigma
a <- ncb.od(pm2.5)$sigma.table$alpha.ncb[1:5]
b <- ncb.od(pm2.5)$sigma.table$beta[1:5]
s <- ncb.od(pm2.5)$sigma.table$sigma[1:5]

alpha.beta.sigma <- t(data.frame(a,b,s))

cplot(pm2.5,2,1,alpha.beta.sigma=alpha.beta.sigma)
cplot(pm2.5,2,1,alpha.beta.sigma=alpha.beta.sigma,regress=TRUE)
data(pm2.5)

## Not run: 
# Use omx function to specify the data for alpha.beta.sigma
pm <- pm2.5

# omx uses OpenMx which does not like periods in data column names
names(pm) <- c('ms_conc_1','ws_conc_1','ms_conc_2','ws_conc_2','frm')

# Fit one-factor measurement error model with FRM sampler as reference
omxfit <- omx(data=pm[,c(5,1:4)],bs.q=c(0.025,0.5,0.975),reps=100)

# Make a calibration plot using the results from omx instead of the default ncb.od
cplot(pm[,c(5,1:4)],1,2,alpha.beta.sigma=omxfit$abs)

## End(Not run)

Extracts the estimated measurement errors assuming there is a constant bias and using the original data.

Description

Extracts the estimated measurement errors assuming there is a constant bias and using the original data values.

Usage

errors.cb(x)

Arguments

Details

Errors should have a zero mean and should be Normally distributed with constant variance for a given method.

Value

Author(s)

Richard A. Bilonick

References

Jaech, J. L. (1985) Statistical Analysis of Measurement Errors. New York: Wiley

.

See Also

cb.pd, ncb.od,lrt

Examples

data(pm2.5)
errors.cb(pm2.5)

Extracts the estimated measurement errors assuming there is no bias and using the original data.

Description

Extracts the estimated measurement errors assuming there is no bias and using the original data values.

Usage

errors.nb(x)

Arguments

Details

Errors should have a zero mean and should be Normally distributed with constant variance for a given method.

Value

Author(s)

Richard A. Bilonick

References

Jaech, J. L. (1985) Statistical Analysis of Measurement Errors. New York: Wiley.

See Also

cb.pd, ncb.od,lrt

Examples

data(pm2.5)
errors.nb(pm2.5)

Extracts the estimated measurement errors assuming there is a nonconstant bias and using the original data values.

Description

Extracts the estimated measurement errors assuming there is a nonconstant bias and using the original data values.

Usage

errors.ncb(x)

Arguments

Details

Errors should have a zero mean and should be Normally distributed with constant variance for a given method.

Value

Author(s)

Richard A. Bilonick

References

Jaech, J. L. (1985) Statistical Analysis of Measurement Errors. New York: Wiley.

See Also

cb.pd, ncb.od,lrt

Examples

data(pm2.5)
errors.ncb(pm2.5)

Likelihood ratio test for all betas equalling one.

Description

Likelihood ratio test statistic - H0: all betas = one.

Usage

lrt(x, M = 40)

Arguments

Details

See Jaech, pp. 204-205.

Value

Author(s)

Richard A. Bilonick

References

Jaech, J. L. (1985) Statistical Analysis of Measurement Errors. New York: Wiley.

See Also

ncb.od,cb.pd,pm2.5

Examples

data(pm2.5)

lrt(pm2.5) # compare all 5 samplers (4 personal and 1 frm)

lrt(pm2.5[,1:4]) # compare only the personal samplers

stem(lrt(pm2.5)$beta.bars) # examine the estimated betas

A modified "pairs" plot with all axes haveing the same range.

Description

Creates all pairwise scatter plots.

Usage

merror.pairs(df,labels=names(df))

Arguments

Details

Creates all pairwise scatter plots with the same range for all axes and adds the diagonal line denote the "line of equality" or "no bias".).

Value

Produces a scatter plot with the calibration curve and titles that include the calibration equation and the scale-bias adjusted imprecision standard deviations.

Author(s)

Richard A. Bilonick

See Also

panel.merror

Examples

data(pm2.5)

# All pairwise plots after square root transformation to Normality
merror.pairs(sqrt(pm2.5))

Compute maximum likelihood estimates of precision.

Description

This is an internal function that computes the maximum likelihood estimates of precision for the constant bias model using paired data.

Usage

mle(v, r, ni)

Arguments

Value

An N vector containing the maximum likelihood estimates of precision.

Author(s)

Richard A. Bilonick

References

Jaech, J. L. (1985) Statistical Analysis of Measurement Errors. New York: Wiley.

Compute squared standard errors for imprecision estimates for the constant bias model using paired data.

Description

This is an internal function that computes squared standard errors for imprecision estimates of the constant bias model using paired data.

Usage

mle.se2(v, r, ni)

Arguments

Details

Computes the squared standard errors for the squared precisions. Before calling this function, compute the MLE's

Value

An N+1 symmetric H matrix. See p. 201 of Jaech, 1985, eq. 6.4.2.

Author(s)

Richard A. Bilonick

References

J. L. Jaech, Statistical Analysis of Measurement Errors, Wiley, New York: 1985.

Compute accuracy estimates and maximum likelihood estimates of precision for the nonconstant bias measurement error model using original data.

Description

Compute accuracy estimates and maximum likelihood estimates of precision for the nonconstant bias measurement error model using original data.

Usage

ncb.od(x, beta = beta.bar(x), M = 40, conf.level = 0.95)

Arguments

Details

Measurement Error Model:

x[i,k] = alpha[i] + beta[i]*mu[k] + epsilon[i,k]

where x[i,k] is the measurement by the ith method for the kth item, i = 1 to N, k = 1 to n, mu[k] is the true value for the kth item, epsilon[i,k] is the normally distributed random error with variance sigma[i] squared for the ith method and the kth item, and alpha[i] and beta[i] are the accuracy parameters for the ith method. The product of the betas is constrained to equal one (equivalently, the geometric average of the beta's is constrained to one). When the betas are all equal to one, the average of the alphas equals zero (equivalently, the sum of the alphas is constrained to zero).

The imprecision for the ith method is sigma[i]. If all alphas are zeroes and all betas are ones, there is no bias. If all betas equal 1, then there is a constant bias. If some of the betas differ from one there is a nonconstant bias. Note that the individual betas are not unique - only ratios of the betas are unique. If you divide all the betas by beta_i, then the betas represent the scale bias of the other devices/methods relative to device/method i. Also, when the betas differ from one, the sigmas are not directly comparable because the measurement scales (size of the units) differ. To make the sigmas comparable, divide them by their corresponding beta. This result is shown as bias.adj.sigma.

By using the original data values, the betas can be estimated and also the process variance, that is, the variance of the true values.

Technically, the alphas and betas describe the measurements in terms of the unknown true values (i.e., the unknown true values can be thought of as a latent variable). The "true values" are ALWAYS unknown (unless you have a real, highly accurate reference method/device). The real goal is to calibrate one device/method in terms of another. This is easily accomplished because each measurement is a function of the same unknow true values. By solving the measurement error model (in expectation) for mu and substituting, any two devices/methods i=1 and i=2 can be be related as:

E[x[1,k]] = alpha[1] - alpha[2]*beta[1]/beta[2] + beta[1]/beta[2]*E[x[2,k]]

or equivalently

E[x[2,k]] = alpha[2] - alpha[1]*beta[2]/beta[1] + beta[2]/beta[1]*E[x[1,k]].

Use cplot to display this calibration curve and the corresponding scale-bias adjusted imprecision standard deviations.

The omx function is to be preferred to ncb.od. omx can accomodate missing measurement data values and can provide both likelihood-based confidence intervals and bootstrapped intervals for all parameters and relevant functions of parameters.

Value

Author(s)

Richard A. Bilonick

References

Jaech, J. L. (1985) Statistical Analysis of Measurement Errors. New York: Wiley.

See Also

cb.pd, lrt

Examples

library(merror)
data(pm2.5)
ncb.od(pm2.5)               # nonconstant bias model using original data values
ncb.od(pm2.5,beta=rep(1,5)) # constant bias model using original data values

Compute full information maximum likelihood estimates of accuracy and precision for the nonconstant bias measurement error model using 'OpenMx'.

Description

Compute full information maximum likelihood (FIML) estimates of accuracy (bias) and precision (imprecision) for the nonconstant bias measurement error model. 'OpenMx' functions are used to construct and fit a one latent factor model. Likelihood-based confidence intervals and bootstrapped confidence interval are can be determined for all model parameters and relevant functions of the model parameters.

Usage

omx(data, rvEst=rep(1,ncol(data)), mubarEst=mean(data[,1]), interval=0.95, 
  reps=500,bs.q=c(0.025,0.975), bs=TRUE)

Arguments

Details

Measurement Error Model:

x_{ik}=\alpha_i+\beta_i\mu_k+\epsilon_{ik}

where x_{ik} is the measurement by the ith of N methods for the kth of n items, i = 1 to N\ge 3, k = 1 to n, \mu_k is the true value for the kth item, \epsilon_{ik} is the normally distributed random error with variance \sigma_i^2 for the ith method and the kth item, and \alpha_i and \beta_i are the accuracy parameters for the ith method. The beta for the first column of data) is set to one. The corresponding alpha is set to 0. These constraints or similar are required for model identification.

The imprecision for the ith method is \sigma_i. If all alphas are zeroes and all betas are ones, there is no bias. If all betas equal 1, then there is a constant bias. If some of the betas differ from one there is a nonconstant bias. Note that the individual betas are not unique - only ratios of the betas are unique. If you divide all the betas by \beta_i, then the betas represent the scale bias of the other devices/methods relative to device/method i. Also, when the betas differ from one, the sigmas are not directly comparable because the measurement scales (size of the units) differ. To make the sigmas comparable, divide them by their corresponding beta.

Technically, the alphas and betas describe the measurements in terms of the unknown true values (i.e., the unknown true values can be thought of as a latent variable). The "true values" are ALWAYS unknown (unless you have a real, highly accurate reference method/device). The real goal is to calibrate one device/method in terms of another. This is easily accomplished because each measurement is a linear function of the same unknown true values. For methods 1 and 2, the calibration curve is given by:

E[x_{1k}]=\left(\alpha_1-\alpha_2\beta_1/\beta_2\right)+\left(\beta_1/\beta_2\right)E[x_{2k}]

or equivalently

E[x_{2k}]=\left(\alpha_2-\alpha_1\beta_2/\beta_1\right)+\left(\beta_2/\beta_1\right)E[x_{1k}]

.

Use cplot, with the alpha.beta.sigma argument specified, to display this calibration curve, calibration equation, and the corresponding scale-bias adjusted imprecision standard deviations.

Note that likelihood confidence intervals and bootstrapped confidence intervals can be returned. Wald-type intervals based on the standard errors are alos available by using the confint function on the returned fit object. See examples.

Value

Note

The following names are used to describe the estimates:

1) a1, a2, a3 and so forth denote the alphas (intercept).

2) b1, b2, b3 and so forth denote the betas (scale or slope).

3) ve1, ve2, ve3 and so forth denote the raw (uncorrected for scale) residual random error variances (imprecision variances).

4) se1 denotes the imprecison standard deviation for the reference method.

5) base2, base3 and so forth denote the scale bias-adjusted imprecision standard deviations on the scale of the reference method.

6) mubar is the estimated mean of the true values on the scale of the reference method.

7) sigma2 is the estimated variance of the true values on the scale of the reference method.

8) sigma is the estimated standard deviation of the true values on the scale of the reference method.

Author(s)

Richard A. Bilonick

See Also

ncb.od, alpha.beta.sigma

Examples

## Not run: 
library(OpenMx)
library(merror)

data(pm2.5)

pm <- pm2.5

# OpenMx does not like periods in data column names
names(pm) <- c('ms_conc_1','ws_conc_1','ms_conc_2','ws_conc_2','frm')

# Fit model with FRM sampler as reference
omxfit <- omx(data=pm[,c(5,1:4)],bs.q=c(0.025,0.5,0.975),reps=100)

# Look at results
summary(omxfit$fit)$parameters[,c(1,5,6)] # Parameter estimates and standard errors
round(omxfit$ci[,1:3],3) # Likelihood-based intervals

# Estimated standard errors and quantiles based on bootstrapped samples
round(omxfit$q.boot,3)

# Wald-type intervals
#   - note variances not standard deviations and different ordering
confint(omxfit$fit) 

omxfit$abs # Use with cplot

# Make a calibration plot using the results from omx instead of the default ncb.od
cplot(pm[,c(5,1:4)],1,2,alpha.beta.sigma=omxfit$abs)

## End(Not run)

Draw diagonal line (line of equality) on merror.pairs plots

Description

This functionis used internally by the function merror.pairs.

Usage

panel.merror(x,y, ...)

Arguments

Value

Draws the diagonal line that represents the "line of equality", i.e., the "no bias model".

Author(s)

Richard A. Bilonick

See Also

merror.pairs

PM 2.5 Concentrations from SCAMP Collocated Samplers

Description

Five filter-based samplers for measuring PM 2.5 concentrations were collocated and provided 77 complete sets of concentrations. This data was collected by the Stuebenville Comprehensive Air Monitoring Program (SCAMP) to check the accuracy and precision of the instruments.

Usage

data(pm2.5)

Format

A data frame with 77 sets of PM 2.5 concentrations (micrograms per cubic meter) from the following 5 samplers:

ms.conc.1

- personal sampler 1 - filter MS

ws.conc.1

- personal sampler 1 - filter WS

ms.conc.2

- personal sampler 2 - filter MS

ws.conc.2

- personal sampler 2 - filter WS

frm

- Federal Reference Method sampler

Source

Stuebenville Comprehensive Air Monitoring Program (SCAMP)

Examples

data(pm2.5)
boxplot(pm2.5)
merror.pairs(pm2.5)

# estimates of accuracy and precision
#   for nonconstant bias model using 
#   original data values
ncb.od(pm2.5)

Computes Grubbs' method of moments estimators of precision for the constant bias model using paired differences.

Description

This is an internal function that computes Grubbs' method of moments estimators of precision for the constant bias model using paired differences

Usage

precision.grubbs.cb.pd(x)

Arguments

Details

See Jaech 1985, Chapters 3 & 4, p. 144 in particular.

Value

Estimated squared imprecision estimates (variances).

Author(s)

Richard A. Bilonick

References

Jaech, J. L. (1985) Statistical Analysis of Measurement Errors. New York: Wiley.

See Also

precision.grubbs.ncb.od, ncb.od, cb.pd,lrt

Computes Grubbs' method of moments estimators of precision for the nonconstant bias model using original data values.

Description

This is an internal function that computes Grubbs' method of moments estimators of precision for the nonconstant bias model using original data values.

Usage

precision.grubbs.ncb.od(x, beta.bar.x = beta.bar(x))

Arguments

Details

See Jaech, p. 184.

Value

Grubbs' method of moments estimates of the squared imprecision (variances).

Author(s)

Richard A. Bilonick

References

Jaech, J. L. (1985) Statistical Analysis of Measurement Errors. New York: Wiley.

See Also

precision.grubbs.cb.pd, ncb.od, cb.pd,lrt

Computes iterative approximation to mle precision estimates for nonconstant bias model using original data.

Description

This is an internal function that computes iterative approximation to mle precision estimates for nonconstant bias model using original data.

Usage

precision.mle.ncb.od(x, M = 20, beta.bars = beta.bar(x), jaech.errors = FALSE)

Arguments

Details

Provides iterative approximation to MLE precision estimates for NonConstant Bias model using Original Data. See Jaech, p. 185-186.

Value

Author(s)

Richard A. Bilonick

References

Jaech, J. L. (1985) Statistical Analysis of Measurement Errors. New York: Wiley.

See Also

precision.grubbs.ncb.od,precision.grubbs.cb.pd

Compute process standard deviation

Description

This function computes the process standard deviation and is used internally by the function precision.grubbs.ncb.od.

Usage

process.sd(x)

Arguments

Details

The process standard deviation is the standard deviation of the true values uncontaminated by measurement error. See Jaech, p. 185.

Value

A scalar containing the method of moments estimate of the process standard deviation.

Author(s)

Richard A. Bilonick

References

Jaech, J. L. (1985) Statistical Analysis of Measurement Errors. New York: Wiley.

See Also

precision.grubbs.ncb.od

Examples

data(pm2.5)
process.sd(pm2.5) # estimate of the sd of the "true values using the method of moments")

Compute process variance.

Description

This is an internal function to compute the process variance.

Usage

process.var.mle(sigma2, s, beta.bars, N, n)

Arguments

Details

See Jaech p. 186 equations 6.37 - 6.3.10.

Value

Estimated process variance.

Author(s)

Richard A. Bilonick

References

Jaech, J. L. (1985) Statistical Analysis of Measurement Errors. New York: Wiley.

Compute process variance but with minor error in Jaech Fortran code.

Description

This is an internal function to compute the process variance that replicates the minor error in Jaech's Fortran code. This allows comparing merror estimates to those shown in Jaech 1985.

Usage

process.var.mle.jaech.err(sigma2, s, beta.bars, N, n)

Arguments

Details

See Jaech p. 186 equations 6.37 - 6.3.10. Jaech p. 288 line 2330 has s[i,j] instead of s[j,j]. Jaech p. 288 line 2410 omits "- 1/d2".

Value

Estimated process variance but replicating minor error in Jaech's Fortran code.

Author(s)

Richard A. Bilonick

References

Jaech, J. L. (1985) Statistical Analysis of Measurement Errors. New York: Wiley.

Spectroscopic and Photometric Galaxy Redshift Measurements

Description

The redshift observations were taken from DEEP 2 Galaxy Redshift Survey.

Usage

data(redshift)

Format

Redshift measurements are usually denoted by z.

A data frame with one spectroscopic redshift measurement and six different photometric measurements (by researcher) for 1432 galaxies:

z_spec

Spectroscopic redshift

z_fink

Photometric redshift - S. Finklestein

z_font

Photometric redshift - A. Fontana

z_pfor

Photometric redshift - J. Pforr

z_salv

Photometric redshift - M. Salvator

z_wikl

Photometric redshift - T. Wiklind

z_wuyt

Photometric redshift - S. Wuyts

Details

Because the photometric methods depend on the same color information, a one-factor measurement error model incuding both the spectroscopic and photomentric measurements would not be a viable model because the photometric measurements would tend to be correlated. A two-factor model would be needed but would require at minimum replicated spectroscopic measurments.

Source

<https://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/20140013340.pdf>

References

Newman, Jeffrey A., Michael C. Cooper, Marc Davis, S. M. Faber, Alison L. Coil, Puragra Guhathakurta, David C. Koo et al. "The DEEP2 Galaxy Redshift Survey: Design, observations, data reduction, and redshifts." The Astrophysical Journal Supplement Series 208, no. 1 (2013): 5.

Examples

library(OpenMx)
library(merror)

data(redshift)
merror.pairs(redshift)

# estimates of accuracy and precision
#   parameters for a one-factor
#   measurement error model
head(redshift)
merror.pairs(redshift)

## Not run: 
red <- omx(redshift[,-1],reps=200) # Drop the spectroscopic measurements

summary(red$fit)
red$ci
red$q.boot

cplot(redshift[,-1],1,2,alpha.beta.sigma=red$abs)

## End(Not run)

Computes the ith iteration for computing the squared imprecision estimates.

Description

This is an internal function that computes the ith iteration for computing the squared imprecision estimates.

Usage

sigma_mle(i, s, sigma2, sigma.mu2, beta.bars, N, n)

Arguments

Details

See Jaech p. 185-186 equations 6.3.1 - 6.3.6.

Value

Estimated squared imprecisions (variances) for the ith iteration.

Author(s)

Richard A. Bilonick

References

Jaech, J. L. (1985) Statistical Analysis of Measurement Errors. New York: Wiley.