Nothing
# TODO: Add comment
#
# Author: schueta6
###############################################################################
#' Internal Function Model 2.
#' @param x (numeric) parameter
powfun2simple <- function(x)
{
list(
predictors=list(beta1 = 1),
variables=list(substitute(x)),
term=function(predictors,variables){
paste( predictors[1],"*",variables[1],"^2")
}
)
}
class(powfun2simple) <- "nonlin"
#' Internal Function Model 3.
#' @param x (numeric) parameter
powfun3 <- function(x)
{
cutval <- 0.95
xmax <- max(x)
list(
predictors=list(beta1 = 1, beta2 = 1),
variables=list(substitute(x)),
term=function(predictors,variables){
paste( "exp(",predictors[1],")*(1+(exp(",predictors[2],")-",cutval,")*",variables[1],"/",xmax,")")
}
)
}
class(powfun3) <- "nonlin"
#' Internal Function Model 3.
#' @param x (numeric) parameter
powfun3simple <- function(x)
{
cutval <- 0.95
list(
predictors=list(beta1 = 1, beta2 = 1),
variables=list(substitute(x)),
term=function(predictors,variables){
paste( predictors[1],"+",predictors[2],"*",variables[1],"^2")
}
)
}
class(powfun3simple) <- "nonlin"
#' Internal Function Model 4.
#' @param x (numeric) parameter 1
#' @param potenz (numeric) parameter 2
powfun4 <- function(x, potenz)
{
cutval <- 0.95
xmax <- max(x)
list(
predictors=list(beta1 = 1, beta2 = 1),
variables=list(substitute(x)),
term=function(predictors,variables){
paste( "exp(",predictors[1],")*(1+(exp(",predictors[2],")-",cutval,")*",variables[1],"/",xmax,")^",potenz)
}
)
}
class(powfun4) <- "nonlin"
#' Internal Function Model 4.
#' @param x (numeric) parameter 1
#' @param potenz (numeric) parameter 2
powfun4simple <- function(x, potenz)
{
list(
predictors=list(beta1 = 1, beta2 = 1),
variables=list(substitute(x)),
term=function(predictors,variables){
paste( "(",predictors[1],"+",predictors[2],"*",variables[1],")^",potenz)
}
)
}
class(powfun4simple) <- "nonlin"
#' Internal Function Model 5.
#' @param x (numeric) parameter 1
powfun5 <- function(x)
{
cutval <- 0.95
xmax <- max(x)
list(
predictors=list(beta1 = 1, beta2 = 1),
variables=list(substitute(x)),
term=function(predictors,variables){
paste( "exp(",predictors[1],")*(1+(exp(",predictors[2],")-",cutval,")*",variables[1],"/",xmax,")")
}
)
}
class(powfun5) <- "nonlin"
#' Internal Function Model 5.
#' @param x (numeric) parameter 1
#' @param K (numeric) parameter 2
powfun5simple <- function(x,K)
{
list(
predictors=list(beta1 = 1, beta2 = 1),
variables=list(substitute(x)),
term=function(predictors,variables){
paste( predictors[1],"+",predictors[2],"*",variables[1],"^", K)
}
)
}
class(powfun5simple) <- "nonlin"
#' Internal Function Model 6.
#' @param x (numeric) parameter 1
powfun6 <- function(x)
{
list(
predictors=list(beta1 = 1, beta2 = 1,beta3=1, pot = 1),
variables=list(substitute(x)),
term=function(predictors,variables){
paste( "exp(",predictors[1],")*(1+(atan(",predictors[2],")/pi-0.5)*",variables[1],
"*exp(",predictors[3],")*(1+exp(",predictors[4],")-(",variables[1],"*exp(",predictors[3],
"))^exp(",predictors[4],"))/exp(",predictors[4],"))")
}
)
}
class(powfun6) <- "nonlin"
#' Internal Function Model 6.
#' @param x (numeric) parameter 1
powfun6simple <- function(x)
{
list(
predictors=list(beta1 = 1, beta2 = 1,beta3=1, pot = 1),
variables=list(substitute(x)),
term=function(predictors,variables){
paste( predictors[1],"+",predictors[2],"*",variables[1],"+",predictors[3],"*",variables[1],"^",predictors[4])
}
)
}
class(powfun6simple) <- "nonlin"
#' Internal Function Model 7.
#' @param x (numeric) parameter 1
powfun7 <- function(x)
{
cutval <- 0.95
xmax <- max(x)
list(
predictors=list(beta = 1, beta2=1, pot = 1),
variables=list(substitute(x)),
term=function(predictors,variables){
paste( "exp(",predictors[1],")*(1+(exp(",predictors[2],")-",cutval,")*(",variables[1],"/",xmax,")^((0.1+10*exp(",predictors[3],"))/(1+exp(",predictors[3],"))))")
}
)
}
class(powfun7) <- "nonlin"
#' Internal Function Model 7.
#' @param x (numeric) parameter 1
powfun7simple <- function(x)
{
list(
predictors=list(beta = 1, beta2=1, pot = 1),
variables=list(substitute(x)),
term=function(predictors,variables){
paste( predictors[1],"+",predictors[2],"*",variables[1],"^",predictors[3])
}
)
}
class(powfun7simple) <- "nonlin"
#' Internal Function Model 8.
#' @param x (numeric) parameter 1
#' @param C (numeric) parameter 2
#' @param signJ (numeric) parameter 3
powfun8 <- function(x, C,signJ)
{
cutval <- 0.95
list(
predictors=list(beta1 = 1, beta2 = 1, beta3 = 1),
variables=list(substitute(x/C)),
term=function(predictors,variables){
paste( "exp(",predictors[1],")*(1+(exp(",predictors[2],")-",cutval,")*",variables[1],")^(",signJ,"*(0.1+10*exp(",predictors[3],"))/(1+exp(",predictors[3],")))")
}
)
}
class(powfun8) <- "nonlin"
#' Internal Function Model 8.
#' @param x (numeric) parameter 1
powfun8simple <- function(x)
{
list(
predictors=list(beta1 = 1, beta2 = 1, pot = 1),
variables=list(substitute(x)),
term=function(predictors,variables){
paste( "(",predictors[1],"+",predictors[2],"*",variables[1],")^",predictors[3],sep=" ")
}
)
}
class(powfun8simple) <- "nonlin"
#' Internal Function Model 9.
#' @param x (numeric) parameter 1
powfun9simple <- function(x)
{
list(
predictors=list(beta1 = 1, pot = 1),
variables=list(substitute(x)),
term=function(predictors,variables){
paste( predictors[1],"*",variables[1],"^",predictors[2],sep=" ")
}
)
}
class(powfun9simple) <- "nonlin"
#' Estimate (Im)Precision-Profiles Modeling the Relationship 'Var ~ Mean'.
#'
#' 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 this variance
#' was observed and the respective degrees of freedom. Provided data must contain at least three columns with
#' this information. There are following variance-functions covered: \cr
#' \itemize{
#' \item{constant variance}{ \eqn{sigma^2}}
#' \item{constant CV}{\eqn{ sigma^2=beta_1*u^2 }}
#' \item{mixed constant, proportional variance}{ \eqn{sigma^2=beta_1+beta_2*u^2}}
#' \item{constrained power model, constant exponent}{ \eqn{sigma^2=(beta_1+beta_2*u)^K}}
#' \item{alternative constrained power model}{ \eqn{sigma^2=beta_1+beta_2*u^K}}
#' \item{alternative unconstrained power model for VF's with a minimum}{ \eqn{sigma^2=beta_1+beta_2*u+beta_3*u^J}}
#' \item{alternative unconstrained power model}{ \eqn{sigma^2=beta_1+beta_2*u^J}}
#' \item{unconstrained power model (default model of Sadler)}{ \eqn{sigma^2=(beta_1 + beta_2 * u)^J}}
#' \item{CLSI EP17 similar model}{ \eqn{sigma^2=beta_1 * u^J}}
#' \item{Exact CLSI EP17 model (fit by linear regression on logarithmic scale)}{ \eqn{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. Note that model 10 fits the same predictor function as model 9 using simple linear
#' regression on a logarithmic scale. Therefore, regression statistics like AIC, deviance etc. is not directly comparable
#'to that of models 1 to 9. Despite these possible redundancies, as 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 reparameterized 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.
#'
#' Functions \code{\link{predict.VFP}} and \code{\link{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 \code{\link{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 \code{\link{plot.VFP}} for details. Function \code{\link{predict.VFP}} makes
#' predictions for specified mean-values based on fitted models.
#'
#' Note, that in cases where a valid solution was found in the re-parameterized space but the final fit with 'gnm' in the
#' original parameter-space did not converge no variance-covariance matrix can be estimated. Therefore, no confidence-intervals
#' will be available downstream.
#'
#' @param Data (data.frame, matrix) containing mean-values, estimated variances and degrees of freedom for each sample
#' @param model.no (integer) in 1:10, may be any combination of these integer-values specifying models 1 to 10 which
#' are consequently fitted to the data; defaults to 7
#' @param K (numeric) value defining the constant used in models 4 and 5; defaults to 2
#' @param 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
#' @param 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
#' @param col.mean (character) string specifying a column/variable in 'Data' containing mean-values; defaults to "Mean"
#' @param col.var (character) string specifying a column/variable in 'Data' containing the variances; Defaults to "VC"
#' @param col.df (character) string specifying a column/variable in 'Data' containing the degrees of freedom; defaults to "DF"
#' @param col.sd (character) string (optional) specifying a column/variable in 'Data' containing the SD-values,
#' used for model 10 to derive more precise CV-values compared to derivation from 'col.var' in case
#' 'col.cv' is not specified directly. Note that column "col.var" must nevertheless be set and existing
#' @param 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.
#' Note that column "col.var" must nevertheless be set and existing
#' @param 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 (\code{.Machine$double.eps})
#' @param ... additional parameters passed forward, e.g. 'vc' of function \code{\link{getMat.VCA}} for selecting a
#' specific variance component in case of 'Data' being a list of 'VCA'-objects (see examples)
#'
#' @return (object) of class 'VFP' with elements:\cr
#' \item{Models}{objects of class 'gnm' representing the fitted model}
#' \item{RSS}{residual sum of squares for each fitted model}
#' \item{AIC}{the Akaike information criterion for each fitted model}
#' \item{Deviance}{the deviance for each fitted model}
#' \item{Formulas}{formula(s) as expression for each fitted model}
#' \item{Mu.Func}{function(s) to transform mean value to re-parameterized scale}
#' \item{Mu}{list of transformed mean values}
#' \item{Data}{the original data set provided to fit model(s)}
#' \item{K}{constant as specified by the user}
#' \item{startvals}{start values as specified by the user}
#' \item{errors}{messages of all errors caught}
#' \item{output}{if 'quiet=TRUE' all output that was redirected to a file}
#' \item{warning}{messages of all warnings caught}
#' \item{notes}{all other notes caught during execution}
#'
#' @author Florian Dufey \email{florian.dufey@@roche.com},
#' Andre Schuetzenmeister \email{andre.schuetzenmeister@@roche.com}
#'
#' @seealso \code{\link{plot.VFP}}, \code{\link{predict.VFP}}, \code{\link{predictMean}}
#'
#' @examples
#' \donttest{
#' # 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)
#' #' }
fit.vfp <- function(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, ...) #.Machine$double.eps, ...)
{
Call <- match.call()
stopifnot(model.no %in% 1:10)
stopifnot(K>0)
stopifnot(is.data.frame(Data) || is.matrix(Data) || is(Data, "list") && all(sapply(Data, class) == "VCA"))
if(is(Data, "list") && all(sapply(Data, class) == "VCA"))
Data <- getMat.VCA(Data, ...)
if(is.matrix(Data))
Data <- as.data.frame(Data)
cn <- colnames(Data)
stopifnot(all(c(col.mean, col.var, col.df) %in% cn))
stopifnot(is.numeric(Data[,col.mean]))
stopifnot(is.numeric(Data[,col.var]))
stopifnot(is.numeric(Data[,col.df]))
stopifnot(is.null(col.sd) || is.character(col.sd) && is.numeric(Data[,col.sd]))
stopifnot(is.null(col.cv) || is.character(col.cv) && is.numeric(Data[,col.cv]))
Data[,col.mean] <- as.numeric(Data[,col.mean])
Data[,col.var] <- as.numeric(Data[,col.var])
Data[,col.df] <- as.numeric(Data[,col.df])
ind.NegVC <- which(Data[,col.var] <= 0)
if(length(ind.NegVC) > 0)
{
message("Variance(s) <= 0 detected! This/these will be set to 'minVC'!")
Data[ind.NegVC,col.var] <- minVC
}
if(!is.null(col.sd))
{
Data[,col.sd] <- as.numeric(Data[,col.sd])
colnames(Data)[which(cn == col.sd)] <- "SD"
}
if(!is.null(col.cv))
{
Data[,col.cv] <- as.numeric(Data[,col.cv])
colnames(Data)[which(cn == col.cv)] <- "CV"
}
if(quiet)
{
file.create("./stdout.log")
}
errors <- warnings <- messages <- notes <- NULL
skip5 <- FALSE
if(all(c(3,5) %in% model.no) && K==2)
{
skip5 <- TRUE # models 3 and 5 are identical in case K=2
if(!quiet)
message("Model 5 will not be fitted because it is identical to model 3 with 'K=2'!")
notes <- c(notes, "Model 5 will not be fitted because it is identical to model 3 with 'K=2'!")
}
if(length(model.no) > 1)
startvals <- NULL
if(col.mean != "Mean")
{
if("Mean" %in% cn)
{
Data <- Data[,-which(colnames(Data) == "Mean")]
cn <- colnames(Data)
}
colnames(Data)[which(cn == col.mean)] <- "Mean"
}
if(col.var != "VC")
{
if("VC" %in% cn)
{
Data <- Data[,-which(colnames(Data) == "VC")]
cn <- colnames(Data)
}
colnames(Data)[which(cn == col.var)] <- "VC"
}
if(col.df != "DF")
{
if("DF" %in% cn)
{
Data <- Data[,-which(colnames(Data) == "DF")]
cn <- colnames(Data)
}
colnames(Data)[which(cn == col.df)] <- "DF"
}
ind.NA <- which(apply(Data[c("Mean", "VC", "DF")], 1, function(x) any(is.na(x))))
if(length(ind.NA) > 0)
{
Data <- Data[-ind.NA,,drop=FALSE]
message(paste(length(ind.NA), "observation(s) with at least one NA value were removed!"))
}
Data <- subset(Data,"Mean">0)
Data <- subset(Data,"VC">0)
pweights <- Data[,"DF"]/2 # prior weight parameters
Mean <- Data[,"Mean"]
RSS <- Df.resid <- Deviance <- AIC <- numeric(10)
Formulas <- c( "sigma^{2}==beta[1]", # model 1
"sigma^{2}==beta[1]*u^{2}", # model 2
"sigma^{2}==beta[1]+beta[2]*u^2", # model 3
paste("sigma^{2}==(beta[1]+beta[2]*u)^{",K,"}"), # model 4, replace K by acutal value of K
paste("sigma^{2}==beta[1]+beta[2]*u^{",K,"}"), # model 5, - " -
"sigma^{2}==beta[1]+beta[2]*u+beta[3]*u^{J}", # model 6
"sigma^{2}==beta[1]+beta[2]*u^{J}", # model 7
"sigma^{2}==(beta[1]+beta[2]*u)^{J}", # model 8
"sigma^{2}==beta[1]*u^{J}", # model 9
"cv==beta[1]*u^{J}") # model 10 (true CLSI EP17 model)
names(RSS) <- names(Deviance) <- names(AIC) <- paste0("Model_", 1:10)
res.gnm1 <- res.gnm2 <- res.gnm3 <- res.gnm4 <- res.gnm5 <- res.gnm6 <- res.gnm7 <- res.gnm8 <- res.gnm9 <- res.mod10 <- NULL
mittel <- Mu <- mu.fun <- vector("list", length=9)
eps <- 1.0e-07
# cutval <- 0.95 #determines lower boundary in models, values between 0 and 1 are possible
# find a crude analytical model for the calculation of the start weights
imin=0
for(i in 0:2){
if (i>0) {mindevalt<-mindev}
mindev <- sum((Data[,"DF"]*(1-weighted.mean(Data[,"VC"]/Data[,"Mean"]^i,pweights)*Data[,"Mean"]^i/Data[,"VC"])^2))
if (i>0){
if(mindev<mindevalt){imin=i}
}
}
startVC <- weighted.mean(Data[,"VC"]/Data[,"Mean"]^imin,pweights)*Data[,"Mean"]^imin+0.05*weighted.mean(Data[,"VC"],pweights)
################################################################################################################################
if(1 %in% model.no)
{
cat("\nModel 1 ")
my.form1 <- VC ~ 1 #constant VC
startvals <- log(weighted.mean(Data[,"VC"],pweights)) #exact solution
tmp.res <- conditionHandler(gnm(formula = my.form1, family = Gamma(link = "log"), data = Data, weights = pweights,
start=startvals, trace=TRUE), file="./stdout.log")
if(tmp.res$status != 2)
{
coeffs <- bt.coef(tmp.res$result, model=1)#
my.form1simple <- VC ~ 1
tmp.res <- conditionHandler(gnm(formula = my.form1simple, family = Gamma(link = "identity"), data = Data, weights = pweights,
start=coeffs, trace=TRUE,iterMax=1), file="./stdout.log")
res.gnm1 <- tmp.res$result
RSS[1] <- sum((res.gnm1$y-res.gnm1$fitted.values)^2)
AIC[1] <- res.gnm1$aic
Df.resid[1] <- res.gnm1$df.residual <- sum(Data[,"DF"])-1
Deviance[1] <- res.gnm1$deviance
tmp.msg <- " ... finished."
}
else
{
errors <- c(errors, paste0("Model 1 (Errors original space): ", tmp.res$errors))
tmp.msg <- " ... could not be fitted due to errors."
}
messages <- c(messages, paste0("Model 1 (Messages): ", tmp.res$messages))
warnings <- c(warnings, paste0("Model 1 (Warnings): ", tmp.res$warnings))
startvals <- NULL
cat(tmp.msg)
}
if(2 %in% model.no)
{
cat("\nModel 2 ")
#mu.fun[[2]] <- function(x) 2*log(x)
mu <- Mu[[2]] <- 2*log(Data[,"Mean"])
my.form2 <- VC ~ offset(mu) # constant CV
startvals <- log(weighted.mean(Data[,"VC"]/Data[,"Mean"]^2,pweights)) #exact solution
tmp.res <- conditionHandler(gnm(formula = my.form2, family = Gamma(link = "log"), data = Data, weights = pweights,
start=startvals, trace=TRUE), file="./stdout.log")
if(tmp.res$status != 2)
{
coeffs <- bt.coef(tmp.res$result, model=2)#
mu <- Data[,"Mean"]^2
my.form2simple <- VC ~ powfun2simple(Mean)-1
tmp.res <- conditionHandler(gnm(formula = my.form2simple, family = Gamma(link = "identity"), data = Data, weights = pweights,
start=coeffs, trace=TRUE,iterMax=0), file="./stdout.log")
res.gnm2 <- tmp.res$result
RSS[2] <- sum((res.gnm2$y-res.gnm2$fitted.values)^2)
AIC[2] <- res.gnm2$aic
Df.resid[2] <- res.gnm2$df.residual <- sum(Data[,"DF"])-1
Deviance[2] <- res.gnm2$deviance
tmp.msg <- " ... finished."
}
else
{
errors <- c(errors, paste0("Model 2 (Errors original space): ", tmp.res$errors))
tmp.msg <- " ... could not be fitted due to errors."
}
messages <- c(messages, paste0("Model 2 (Messages): ", tmp.res$messages))
warnings <- c(warnings, paste0("Model 2 (Warnings): ", tmp.res$warnings))
startvals <- NULL
cat(tmp.msg)
}
# models 3 4, and 5 are covered by models and, therefore, do not need to be fitted separately if the more
# general models optimizing the exponent are fitted per default
if (3 %in% model.no) {
cat("\nModel 3 ")
mu <- Mu[[3]] <- Data[, "Mean"]^2
mumax <- max(mu)
startweights <- pweights/startVC^2
if(is.null(startvals))
{
st <- lm(VC ~ mu, data=Data,weights=startweights)
startvals <- t.coef(st$coefficients, Maxi=mumax,model=3)
}
else
{
startvals <- t.coef(startvals,Maxi=mumax,model=3)
}
my.form3 <- VC~ -1 + powfun3(mu) #quadratic model
tmp.res <- conditionHandler(gnm(formula = my.form3, family = Gamma(link = "identity"), data = Data, weights = pweights, start=startvals, trace=TRUE), file="./stdout.log")
if(tmp.res$status != 2)
{
coeffs <- bt.coef(tmp.res$result, model=3)
AIC[3] <- tmp.res$result$aic
tmp.df <- tmp.res$result$df.residual
Deviance[3] <- tmp.res$result$deviance
my.form3simple <- VC ~ powfun3simple(Mean)-1
tmp.res <- conditionHandler(gnm(formula = my.form3simple, family = Gamma(link = "identity"), data = Data, weights = pweights,
start=coeffs, trace=TRUE,iterMax=1000), file="./stdout.log")
if(tmp.res$status != 2)
{
res.gnm3 <- tmp.res$result
RSS[3] <- sum((res.gnm3$y-res.gnm3$fitted.values)^2)
tmp.msg <- " ... finished."
}
else
{
res.gnm3 <- list(coefficients=coeffs, aic=AIC[3], deviance=Deviance[3], df.residual=tmp.df)
tmp.msg <- paste(" ... estimation of covariance-matrix failed, CIs will not be available!", sep="")
}
}
else
{
errors <- c(errors, paste0("Model 3 (Errors transformed space): ", tmp.res$errors))
tmp.msg <- " ... could not be fitted due to errors."
}
messages <- c(messages, paste0("Model 3 (Messages): ", tmp.res$messages))
warnings <- c(warnings, paste0("Model 3 (Warnings): ", tmp.res$warnings))
startvals <- NULL
cat(tmp.msg)
}
if(4 %in% model.no)
{
cat("\nModel 4 ")
mu <- Mu[[4]] <- Data[,"Mean"]
mumax=max(mu)
startweights <- pweights/startVC^(2/K)
VCrootK <- Data[,"VC"]^(1/K)
my.form4 <- VC~-1+powfun4(mu, K) #fixed power model
if (is.null(startvals))
{
st <- lm(VCrootK ~ mu, data=Data,weights=startweights)
startvals <- t.coef(st$coefficients,K=K, Maxi=mumax,model=4)
}
else
{
startvals <- t.coef(startvals, K=K, Maxi=mumax, model=4)
}
tmp.res <- conditionHandler(gnm(formula = my.form4, family = Gamma(link = "identity"), data = Data, weights = pweights,
start=startvals, trace=TRUE), file="./stdout.log")
if(tmp.res$status != 2)
{
coeffs <- bt.coef(tmp.res$result,K=K, model=4)#
AIC[4] <- tmp.res$result$aic
Deviance[4] <- tmp.res$result$deviance
my.form4simple <- VC~-1+powfun4simple(Mean, K) #fixed power model
tmp.res <- conditionHandler(gnm(formula = my.form4simple, family = Gamma(link = "identity"), data = Data, weights = pweights,
start=coeffs, trace=TRUE,iterMax=0), file="./stdout.log")
if(tmp.res$status != 2)
{
res.gnm4 <- tmp.res$result
RSS[4] <- sum((res.gnm4$y-res.gnm4$fitted.values)^2)
Df.resid[4] <- res.gnm4$df.residual <- sum(Data[,"DF"])-2
tmp.msg <- " ... finished."
}
else
{
Df.resid[4] <- sum(Data[,"DF"])-2
res.gnm4 <- list(coefficients=coeffs, aic=AIC[4], deviance=Deviance[4], df.residual=Df.resid[4])
tmp.msg <- paste(" ... estimation of covariance-matrix failed, CIs will not be available!", sep="")
}
}
else
{
res.gnm4 <- NULL
errors <- c(errors, paste0("Model 4 (Errors transformed space): ", tmp.res$errors))
tmp.msg <- " ... could not be fitted due to errors."
}
messages <- c(messages, paste0("Model 4 (Messages): ", tmp.res$messages))
warnings <- c(warnings, paste0("Model 4 (Warnings): ", tmp.res$warnings))
startvals <- NULL # re-set start values
cat(tmp.msg)
}
if(5 %in% model.no){
cat("\nModel 5 ")
if(skip5)
cat(" ... skipped!")
else
{
mu <- Mu[[5]] <- Data[, "Mean"]^K
mumax <- max(mu)
startweights <- pweights/startVC^2
if (is.null(startvals))
{
st <- lm(VC ~ mu, data=Data,weights=startweights)
startvals <- t.coef(st$coefficients,Maxi=mumax, model=5)
}
else
{
startvals <-t.coef(startvals,Maxi=mumax,Model=5)
}
my.form5 <- VC~-1+powfun5(mu) #const model, requires a fixed power on input
tmp.res <- conditionHandler(gnm(formula = my.form5, family = Gamma(link = "identity"), data = Data, weights = pweights,
start=startvals, trace=TRUE), file="./stdout.log")
if(tmp.res$status != 2)
{
coeffs <- bt.coef(tmp.res$result,K=K, model=5)
AIC[5] <- tmp.res$result$aic
Deviance[5] <- tmp.res$result$deviance
my.form5simple <- VC~powfun5simple(Mean,K)-1
tmp.res <- conditionHandler(gnm(formula = my.form5simple, family = Gamma(link = "identity"), data = Data, weights = pweights,
start=coeffs, trace=TRUE,iterMax=0), file="./stdout.log")
if(tmp.res$status != 2)
{
res.gnm5 <- tmp.res$result
RSS[5] <- sum((res.gnm5$y-res.gnm5$fitted.values)^2)
Df.resid[5] <- res.gnm5$df.residual <- sum(Data[,"DF"])-2
tmp.msg <- " ... finished."
}
else
{
Df.resid[5] <- sum(Data[,"DF"])-2
res.gnm5 <- list(coefficients=coeffs, aic=AIC[5], deviance=Deviance[5], df.residual=Df.resid[5])
tmp.msg <- paste(" ... estimation of covariance-matrix failed, CIs will not be available!", sep="")
}
}
else
{
res.gnm5 <- NULL
errors <- c(errors, paste0("Model 5 (Errors transformed space): ", tmp.res$errors))
tmp.msg <- " ... could not be fitted due to errors.\n"
}
messages <- c(messages, paste0("Model 5 (Messages): ", tmp.res$messages))
warnings <- c(warnings, paste0("Model 5 (Warnings): ", tmp.res$warnings))
startvals <- NULL
cat(tmp.msg)
}
}
if(6 %in% model.no)
{
cat("\nModel 6 ")
mu <- Mu[[6]] <- Data[,"Mean"]
my.form6 <- VC~-1+ powfun6(mu) #linear + variable power model
results <- vector("list", 5)
par.ind <- 4
res.ind <- 1
res.gnm6 <- tmp <- NULL
Parms <- vector("list", 5)
SVhist <- matrix(ncol=4, nrow=0)
startweights <- pweights/startVC^2
tmp.msg <- " ... could not be fitted due to errors."
if (is.null(startvals)) # use random start values for gnm
{
tempdev=1.0e100
res.gnm6 <- NULL
for(ldJ in 0:4)
{
J <- 2^ldJ+0.1
muJ <- mu^J
st <- lm(VC ~ mu+muJ, data=Data,weights=startweights)
startvals <-t.coef(c(st$coefficients,J),model=6)
tmp.res <- conditionHandler(gnm(formula = my.form6, family = Gamma(link = "identity"), data = Data, weights = pweights,
start=startvals, trace=TRUE), file="./stdout.log")
if(tmp.res$status != 2)
{
if(!is.null(tmp.res$result) && tmp.res$result$deviance<tempdev)
{
res.gnm6 <- tmp.res$result
tempdev <- res.gnm6$deviance
tmp.msg <- " ... finished."
}
}
else
errors <- c(errors, paste0("Model 6 (Errors transformed space): ", tmp.res$errors))
}
}
else # user-specified start values
{
startvals <- t.coef(startvals,model=6)
tmp.res <- conditionHandler(gnm(formula = my.form6, family = Gamma(link = "identity"), data = Data, weights = pweights,
start=startvals, trace=TRUE), file="./stdout.log")
if(tmp.res$status != 2)
{
res.gnm6 <- tmp.res$result
tmp.msg <- " ... finished.\n"
}
else
{
errors <- c(errors, paste0("Model 6 (Errors transformed space): ", tmp.res$errors))
tmp.msg <- " ... could not be fitted due to errors."
}
}
if(!is.null(res.gnm6))
{
coeffs <- bt.coef(res.gnm6, model=6)#
AIC[6] <- res.gnm6$aic
Deviance[6] <- res.gnm6$deviance
my.form6simple <- VC~-1+powfun6simple(Mean)
tmp.res <- conditionHandler(gnm(formula = my.form6simple, family = Gamma(link = "identity"), data = Data, weights = pweights,
start=coeffs, trace=TRUE, iterMax=0), file="./stdout.log")
if(tmp.res$status != 2)
{
res.gnm6 <- tmp.res$result
RSS[6] <- sum((res.gnm6$y-res.gnm6$fitted.values)^2)
Df.resid[6] <- res.gnm6$df.residual <- sum(Data[,"DF"])-4
}
else
{
Df.resid[6] <- sum(Data[,"DF"])-4
res.gnm6 <- list(coefficients=coeffs, aic=AIC[6], deviance=Deviance[6], df.residual=Df.resid[6])
tmp.msg <- paste(" ... estimation of covariance-matrix failed, CIs will not be available!", sep="")
}
}
messages <- c(messages, paste0("Model 6 (Messages): ", tmp.res$messages))
warnings <- c(warnings, paste0("Model 6 (Warnings): ", tmp.res$warnings))
startvals <- NULL
cat(tmp.msg)
}
if(7 %in% model.no)
{
cat("\nModel 7 ")
mu <- Mu[[7]] <- Data[,"Mean"]
startweights <- pweights/startVC^2
tmp.msg <- " ... could not be fitted due to errors."
if (is.null(startvals) || startvals[3]<0 ) #negative power leads negative sigma^2 for small means
{
tempdev <- 1.0e100
res.gnm7 <- NULL
for(ldJ in -3:3){
J <- 2^ldJ
muJ <- mu^J
st <- lm(VC ~ muJ, data=Data,weights=startweights)
startvals <-t.coef(c(st$coefficients,J),Maxi=max(muJ),model=7)
my.form7 <- VC ~ -1+powfun7(mu)
tmp.res <- conditionHandler(gnm(formula=my.form7, family=Gamma(link = "identity"), data=Data,
weights=pweights, start=startvals, trace=TRUE), file="./stdout.log")
if(tmp.res$status != 2)
{
if(!is.null(tmp.res$result) && tmp.res$result$deviance<tempdev)
{
res.gnm7 <- tmp.res$result
tempdev <- res.gnm7$deviance
tmp.msg <- " ... finished."
}
}
else
errors <- c(errors, paste0("Model 7 (Errors transformed space): ", tmp.res$errors))
}
}
else
{
my.form7 <- VC ~ -1+powfun7(mu)
startvals[3] <- max(0.1,min(10,startvals[3]))#restrict starting values to 0.1<startvals[3] <10
muJmax<-max(mu^startvals[3])
startvals <- t.coef(startvals,Maxi=muJmax,model=7)
tmp.res <- conditionHandler(gnm(formula = my.form7, family = Gamma(link = "identity"), data = Data, weights = pweights,
start=startvals, trace=TRUE), file="./stdout.log")
if(tmp.res$status != 2)
{
res.gnm7 <- tmp.res$result
tmp.msg <- " ... finished."
}
else
{
errors <- c(errors, paste0("Model 7 (Errors transformed space): ", tmp.res$errors))
tmp.msg <- " ... could not be fitted due to errors."
}
}
if(!is.null(res.gnm7))
{
coeffs <- bt.coef(res.gnm7, model=7)#
AIC[7] <- res.gnm7$aic
Deviance[7] <- res.gnm7$deviance
my.form7simple <- VC~-1+powfun7simple(Mean)
tmp.res <- conditionHandler(gnm(formula = my.form7simple, family = Gamma(link = "identity"), data = Data, weights = pweights,
start=coeffs, trace=TRUE,iterMax=0), file="./stdout.log")
if(!is.null(tmp.res))
{
res.gnm7 <- tmp.res$result
RSS[7] <- sum((res.gnm7$y-res.gnm7$fitted.values)^2)
Df.resid[7] <- res.gnm7$df.residual <- sum(Data[,"DF"])-3
}
else
{
Df.resid[7] <- sum(Data[,"DF"])-3
res.gnm7 <- list(coefficients=coeffs, aic=AIC[7], deviance=Deviance[7], df.residual=Df.resid[7])
tmp.msg <- paste(" ... estimation of covariance-matrix failed, CIs will not be available!", sep="")
}
}
messages <- c(messages, paste0("Model 7 (Messages): ", tmp.res$messages))
warnings <- c(warnings, paste0("Model 7 (Warnings): ", tmp.res$warnings))
startvals <- NULL
cat(tmp.msg)
}
if(8 %in% model.no)
{
cat("\nModel 8 ")
mu <- Mu[[8]] <- Data[,"Mean"]
C <- attr(Mu[[8]], "Constant") <- max(mu)
my.form8 <- VC~-1+powfun8(mu, C)
devns <- rep(NA, 1e3)
tmp.msg <- " ... could not be fitted due to errors."
if (is.null(startvals)) # random start values
{
tempdev <- 1.0e100
res.gnm8 <- NULL
for(ldJ in -3:3){
J <- 2^ldJ
signJ <- 1
my.form8 <- VC~-1+powfun8(mu, C,signJ)
startweights <- pweights/startVC^(2/J)
VCrootJ <- Data[,"VC"]^(1/J)
st <- lm(VCrootJ ~ mu, data=Data,weights=startweights)
startvals <- t.coef(c(st$coefficients,signJ*J),Maxi=C,signJ=signJ,model=8)
tmp.res <- conditionHandler(gnm(formula = my.form8, family = Gamma(link = "identity"),
data = Data, weights = pweights ,start=startvals,
trace=TRUE), file="./stdout.log")
if(tmp.res$status != 2)
{
if(!is.null(tmp.res$result))
{
if(tmp.res$result$deviance<tempdev)
{
res.gnm8 <- tmp.res$result
sgnJ <- signJ
tempdev <- res.gnm8$deviance
tmp.msg <- " ... finished."
}
}
}
else
errors <- c(errors, paste0("Model 8 (Errors transformed space): ", tmp.res$errors))
}
for(ldJ in -3:3)
{
J <- -2^ldJ
signJ <- -1
my.form8 <- VC~-1+powfun8(mu, C,signJ)
startweights <- pweights/startVC^(2/J)
VCrootJ <- Data[,"VC"]^(1/J)
st <- lm(VCrootJ ~ mu, data=Data,weights=startweights)
startvals <- t.coef(c(st$coefficients,signJ*J),Maxi=C,signJ=signJ,model=8)
tmp.res <- conditionHandler(gnm(formula = my.form8, family = Gamma(link = "identity"),
data = Data, weights = pweights ,start=startvals,
trace=TRUE), file="./stdout.log")
if(tmp.res$status != 2)
{
if(!is.null(tmp.res$result))
{
if(tmp.res$result$deviance<tempdev)
{
res.gnm8 <- tmp.res$result
sgnJ <- signJ
tempdev <- res.gnm8$deviance
tmp.msg <- " ... finished."
}
}
}
else
errors <- c(errors, paste0("Model 8 (Errors transformed space): ", tmp.res$errors))
}
}
else # user-specified startvalues
{
if(startvals[3]>0){ #restrict range of starting values to 0.1 <|J| <10
startvals[3] <- max(0.1,min(10,startvals[3]))
signJ <- 1
}
else{
signJ <- -1
startvals[3] <- -max(0.1,min(10,-startvals[3]))
}
my.form <- VC~-1+powfun8(Mean, C,signJ)
startvals <-t.coef(startvals,Maxi=C,signJ=signJ,model=8)
tmp.res <- conditionHandler(gnm(formula = my.form8, family = Gamma(link = "identity"),
data = Data, weights = pweights ,start=startvals,
trace=TRUE), file="./stdout.log")
if(tmp.res$status != 2)
{
if(!is.null(tmp.res$result))
{
res.gnm8 <- tmp.res$result
sgnJ <- signJ
tmp.msg <- " ... finished."
}
}
else
{
errors <- c(errors, paste0("Model 8 (Errors transformed space): ", tmp.res$errors))
tmp.msg <- " ... could not be fitted due to errors."
}
}
if(!is.null(res.gnm8))
{
coeffs <- bt.coef(res.gnm8,signJ=sgnJ, model=8)
AIC[8] <- res.gnm8$aic
Deviance[8] <- res.gnm8$deviance
my.form8simple <- VC~-1+powfun8simple(Mean)
tmp.res <- conditionHandler(gnm(formula = my.form8simple, family = Gamma(link = "identity"), data = Data, weights = pweights,
start=coeffs, trace=TRUE,iterMax=0), file="./stdout.log")
if(!is.null(tmp.res))
{
res.gnm8 <- tmp.res$result
RSS[8] <- sum((res.gnm8$y-res.gnm8$fitted.values)^2)
Df.resid[8] <- res.gnm8$df.residual <- sum(Data[,"DF"])-3
}
else
{
Df.resid[8] <- sum(Data[,"DF"])-3
res.gnm8 <- list(coefficients=coeffs, aic=AIC[8], deviance=Deviance[8], df.residual=Df.resid[8])
tmp.msg <- paste(" ... estimation of covariance-matrix failed, CIs will not be available!", sep="")
}
}
messages <- c(messages, paste0("Model 8 (Messages): ", tmp.res$messages))
warnings <- c(warnings, paste0("Model 8 (Warnings): ", tmp.res$warnings))
startvals <- NULL
cat(tmp.msg)
}
if(9 %in% model.no)
{
cat("\nModel 9 ")
mu <- Mu[[9]] <- log(Data[,"Mean"])
logVC <- log(Data[,"VC"])
my.form9 <- formula(paste("VC", "~ mu", sep = " ")) # variable power model without intercept
if (is.null(startvals)) # random start values
{
startweights <- pweights
st <- lm(logVC ~ mu, data=Data,weights=startweights)
startvals <- c(st$coefficients[1],st$coefficients[2])
}
else
{
startvals <-t.coef(startvals,model=9)
}
tmp.res <- conditionHandler(gnm(formula = my.form9, family = Gamma(link = "log"), data = Data, weights = pweights,
start=startvals, trace=TRUE), file="./stdout.log")
if(tmp.res$status != 2)
{
res.gnm9 <- tmp.res$result
tmp.msg <- " ... finished."
}
else
{
errors <- c(errors, paste0("Model 9 (Errors transformed space): ", tmp.res$errors))
tmp.msg <- " ... could not be fitted due to errors."
}
if(!is.null(res.gnm9))
{
coeffs <- bt.coef(res.gnm9, model=9)#
AIC[9] <- res.gnm9$aic
Deviance[9] <- res.gnm9$deviance
my.form9simple <- VC~-1+powfun9simple(Mean)
tmp.res <- conditionHandler(gnm(formula = my.form9simple, family = Gamma(link = "identity"), data = Data, weights = pweights,
start=coeffs, trace=TRUE, iterMax=0), file="./stdout.log")
if(!is.null(tmp.res ))
{
res.gnm9 <- tmp.res$result
RSS[9] <- sum((res.gnm9$y-res.gnm9$fitted.values)^2)
Df.resid[9] <- res.gnm9$df.residual <- sum(Data[,"DF"])-2
}
else
{
Df.resid[9] <- sum(Data[,"DF"])-2
res.gnm9 <- list(coefficients=coeffs, aic=AIC[9], deviance=Deviance[9], df.residual=Df.resid[9])
tmp.msg <- paste(" ... estimation of covariance-matrix failed, CIs will not be available!", sep="")
}
}
messages <- c(messages, paste0("Model 9 (Messages): ", tmp.res$messages))
warnings <- c(warnings, paste0("Model 9 (Warnings): ", tmp.res$warnings))
startvals <- NULL
cat(tmp.msg)
}
if(10 %in% model.no)
{
cat("\nModel 10")
if(!is.null(col.cv)) # minimize error propagation effects by using CV-values directly
{
res.mod10 <- fit.EP17(x=Data$Mean, y=Data$CV, DF=Data$DF, typeY="cv")
}
else
{
if(!is.null(col.sd)) # second best option
{
Data$CV <- 100*Data$SD/Data$Mean
res.mod10 <- fit.EP17(x=Data$Mean, y=Data$SD, DF=Data$DF, typeY="sd")
}
else # worst option, since VC --> SD --> CV
{
Data$CV <- 100*sqrt(Data$VC)/Data$Mean
res.mod10 <- fit.EP17(x=Data$Mean, y=Data$VC, DF=Data$DF, typeY="vc")
}
}
RSS[10] <- res.mod10$RSS
AIC[10] <- res.mod10$AIC
Df.resid[10] <- sum(Data[,"DF"])-2
Deviance[10] <- res.mod10$deviance
cat(" ... finished.")
}
cat("\n")
if(quiet) # read logs and remove temporary files
{
#try(close(outputLog), silent=TRUE)
output <- scan("./stdout.log", sep="\n", what="character", quiet=TRUE)
try(file.remove("./stdout.log"), silent=TRUE)
}
else
output <- NULL
mod <- list(
model1 =res.gnm1,
model2 =res.gnm2,
model3 =res.gnm3,
model4 =res.gnm4,
model5 =res.gnm5,
model6 =res.gnm6,
model7 =res.gnm7,
model8 =res.gnm8,
model9 =res.gnm9,
model10=res.mod10)
ind <- which(!sapply(mod, is.null))
res <- list(
Call=Call,
Models=mod[ind],
RSS=RSS[ind],
AIC=AIC[ind],
GoF.pval=sapply(mod[ind], function(x) pchisq(x$deviance, x$df.residual, lower.tail=FALSE)),
DF.resid=Df.resid[ind],
Deviance=Deviance[ind],
Formulas=Formulas[ind],
Data=Data,
K=K,
startvals=startvals,
errors=errors,
output=output,
warnings=warnings,
notes=notes)
class(res) <- "VFP"
return(res)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.