fit.vfp: Estimate (Im)Precision-Profiles Modeling the Relationship 'Var ~ Mean'.

Description

This function fits one out of ten or any subset of ten non-linear functions at once. This is done on precision data consisting of information about the variance, concentration at which these variance was observed and the respective degrees of freedom. Provided data must contain at least three column with these information. There are following variance-functions covered:

constant variance $sigma^2$
constant CV$ sigma^2=beta_1*u^2 $
mixed constant, proportional variance $sigma^2=beta_1+beta*u^2$
constrained power model, constant exponent $sigma^2=(beta_1+beta_2*u)^C$
alternative constrained power model $sigma^2=beta_1+beta_2*u^C$
alternative unconstrained power model for VF's with a minimum $sigma^2=beta_1+beta_2*u+beta_3*u^J$
alternative unconstrained power model $sigma^2=beta_1+beta_2*u^J$
unconstrained power model (default model of Sadler) $sigma^2=(beta_1 + beta_2 * u)^J$
CLSI EP17 similar model $sigma^2=beta_1 * u^J$
Exact CLSI EP17 model (fit by linear regression on logarithmic scale) $cv=beta_1 * u^J$

Fitting all ten models is somehow redundant if constant 'C' is chosen to be equal to 2, since models 3 and 5 are equivalent and these are constrained versions of model 7 where the exponent is also estimated. The latter also applies to model 4 which is a constrained version of model 8. Nevertheless, since computation time is not critical here for typical precision-profiles (of in-vitro diagnostics precision experiments) we chose to offer batch-processing as well. During computation, all models are internally reparametrized so as to guarantee that the variance function is positive in the range 'u' from 0 to 'u_max'. In models 7 and 8, 'J' is restricted to 0.1<J<10 to avoid the appearance of sharp hooks. Occasionally, these restrictions may lead to a failure of convergence. This is then a sign that the model parameters are on the boundary and that the model fits the data very badly. This should not be taken as reason for concern. It occurs frequently for model 6 when the variance function has no minimum, which is normally the case.

Usage

fit.vfp(Data, model.no = 7, K = 2, startvals = NULL, quiet = T,
  col.mean = "Mean", col.var = "VC", col.df = "DF", col.sd = NULL,
  col.cv = NULL, minVC = NA, ...)

Arguments

Data

(data.frame, matrix) containing mean-values, estimated variances and degrees of freedom for each sample

model.no

(integer) in 1:10, may be any combination of these integer-values specifying models 1 to 9 which are consequently fitted to the data

(numeric) value defining the constant used in models 4 and 5

startvals

(numeric) vector of start values for the optimization algorithm, only used if a single model is specified by the user, otherwise, start values will be sought automatically

quiet

(logical) TRUE = all output messages (warnings, errors) and regular console output is suppressed and can be found in elements "stderr" and "stdout" of the returned object

col.mean

(character) string specifying a column/variable in 'Data' containing mean-values

col.var

(character) string specifying a column/variable in 'Data' containing the variances

col.df

(character) string specifying a column/variable in 'Data' containing the degrees of freedom

col.sd

(character) string (optional) specifying a column/variable in 'Data' containing the SD-values, used for model 10 to derive more precice CV-value compared to derivation from 'col.var' in case 'col.cv' is not specified directly

col.cv

(character) string (optional) specifying a column/variable in 'Data' containing the CV-values, if missing, first 'col.sd' is checked, if missing 'col.var' used to derive per-sample CV-values

minVC

(numeric) value assigned to variances being equal to zero, defaults to NA, which results in removing these observations; could also be set to the smallest possible positive double (.Machine$double.eps)

...

additional parameters passed forward, e.g. 'vc' of function getMat.VCA for selecting a specific variance component in case of 'Data' being a list of 'VCA'-objects (see examples)

Value

(object) of class 'VFP' with elements:

Models

objects of class 'gnm' representing the fitted model

RSS

residual sum of squares for each fitted model

AIC

the Akaike information criterion for each fitted model

Deviance

the deviance for each fitted model

Formulas

formula(s) as expression for each fitted model

Mu.Func

function(s) to transform mean value to re-parameterized scale

list of transformed mean values

Data

the original data set provided to fit model(s)

constant as specified by the user; must be >0

startvals

start values as specified by the user

errors

messages of all errors caught

output

if 'quiet=TRUE' all output that was redirected to a file

warning

messages of all warnings caught

notes

all other notes caught during execution

Details

Functions predict.VFP and predictMean are useful to make inference based on a fitted model. It is possible to derive concentrations at which a predefined variability is reached, which is sometimes referred to as "functional sensitivity" and/or "limit of quantitation" (LoQ). Funtion predictMean returns the fitted value at which a user-defined variance ("vc"), SD or CV is reached with its corresponding 100(1-alpha)% CI derived from the CI of the fitted model. The plotting method for objects of class 'VFP' can automatically add this information to a plot using arguments 'Prediction' and 'Pred.CI' (see plot.VFP for details. Function predict.VFP makes predictions for specified mean-values based on fitted models.

Examples

Run this code

# NOT RUN {
# load VCA-package and data
library(VCA)
data(VCAdata1)
# perform VCA-anaylsis
lst <- anovaVCA(y~(device+lot)/day/run, VCAdata1, by="sample")
# transform list of VCA-objects into required matrix
mat <- getMat.VCA(lst)		# automatically selects "total"
mat

# fit all 9 models batch-wise
res <- fit.vfp(model.no=1:10, Data=mat)

# if 'mat' is not required for later usage, following works
# equally well
res2 <- fit.vfp(lst, 1:10)

# plot best-fitting model
plot(res)
plot(res, type="cv")
plot(res, type="cv", ci.type="lines", ci.col="red",
		Grid=list(col="wheat"), Points=list(pch=2, lwd=2, col="black"))

# now derive concentation at which a specific reproducibility-
# imprecision of 10\% is reached and add this to the plot
pred <- plot(res, type="cv", ci.type="band", 
				ci.col=as.rgb("red", .25), Grid=list(col="orange"), 
				Points=list(pch=2, lwd=2, col="black"),
				Prediction=list(y=10, col="red"), Pred.CI=TRUE)

# (invisibly) returned object contains all relevant information
pred

# same for repeatability
mat.err <- getMat.VCA(lst, "error")
res.err <- fit.vfp(1:10, Data=mat.err)

# without extracting 'mat.err'
res.err2 <- fit.vfp(lst, 1:10, vc="error")

plot(res.err)

#######################################################################
# another example using CA19_9 data from CLSI EP05-A3

data(CA19_9)

# fit reproducibility model to data
fits.CA19_9 <- anovaVCA(result~site/day, CA19_9, by="sample")

# fit within-laboratory-model treating site as fixed effect
fits.ip.CA19_9 <- anovaMM(result~site/(day), CA19_9, by="sample")

# the variability "between-site" is not part of "total"
fits.ip.CA19_9[[1]]
fits.CA19_9[[1]]

# extract repeatability
rep.CA19_9 <- getMat.VCA(fits.CA19_9, "error")

# extract reproducibility
repro.CA19_9 <- getMat.VCA(fits.CA19_9, "total")

# extract intermediate-precision (within-lab)
ip.CA19_9 <- getMat.VCA(fits.ip.CA19_9, "total")

# fit model (a+bX)^C (model 8) to all three matrices

mod8.repro 	<- fit.vfp(repro.CA19_9, 8)
mod8.ip		<- fit.vfp(ip.CA19_9, 8)
mod8.rep		<- fit.vfp(rep.CA19_9, 8)

# plot reproducibility precision profile first
# leave enough space in right margin for a legend
plot(mod8.repro, mar=c(5.1, 7, 4.1, 15), 
		type="cv", ci.type="none", Model=FALSE,
		Line=list(col="blue", lwd=3), 
		Points=list(pch=15, col="blue", cex=1.5),  
		xlim=c(10, 450), ylim=c(0,10),
		Xlabel=list(text="CA19-9, kU/L (LogScale) - 3 Patient Pools, 3 QC Materials",
				cex=1.5), Title=NULL,
		Ylabel=list(text="% CV", cex=1.5),
		Grid=NULL, Crit=NULL, log="x")

# add intermediate precision profile
plot	(mod8.ip, type="cv", add=TRUE, ci.type="none",
		Points=list(pch=16, col="deepskyblue", cex=1.5),
		Line=list(col="deepskyblue", lwd=3), log="x")

# add repeatability precision profile
plot(	mod8.rep, type="cv", add=TRUE, ci.type="none",
		Points=list(pch=17, col="darkorchid3", cex=1.5),
		Line=list(col="darkorchid3", lwd=3), log="x")

# add legend to right margin
legend.rm( x="center", pch=15:17, col=c("blue", "deepskyblue", "darkorchid3"),
		cex=1.5, legend=c("Reproducibility", "Within-Lab Precision", "Repeatability"),
		box.lty=0)
#' 
# }

Run the code above in your browser using DataLab