diff --git a/.Rbuildignore b/.Rbuildignore new file mode 100644 index 0000000..c84a6de --- /dev/null +++ b/.Rbuildignore @@ -0,0 +1,4 @@ +\.lnk$ +\.zip$ + +^\.git \ No newline at end of file diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..b3d2d2c --- /dev/null +++ b/.gitignore @@ -0,0 +1,7 @@ +*Thumbs.db +*.lnk +~* +*~ +.Rhistory +vignettes/* +!vignettes/using_abrem.Rnw \ No newline at end of file diff --git a/ChangeLog b/ChangeLog new file mode 100644 index 0000000..1beb19e --- /dev/null +++ b/ChangeLog @@ -0,0 +1,44 @@ +2014-04-17 Jurgen Symynck + + * upgraded abrem 0.1.14 to 0.1.16: + No changes, trying to counter build problems on r-forge (having to do with linking to RcppArmaddillo) + +2014-04-17 Jurgen Symynck + + * upgraded abrem 0.1.12 to 0.1.14: + General man page updates and reorganization, more error checking, added support for subtracting threshold parameters from the data before plotting, replaced "bernard" with "benard", renamed option "blives" to "unrel", renamed option "conf.n" to "unrel.n", added some debugging datasets and a mixed model (synthetic) dataset. + Many more modifications to the code for debugging and (currently) undocumented features. (choose fit calculation code, choose confidence calculation code, plot superSMITH reports ...) + Removed the two "vignettes" that used to be accessible with browseVignettes(), but not anymore since R 3.1. + + * Abrem.R (Abrem): added support for arguments of class "numeric", + * plot.abrem.R (plot.abrem): "main" title now plots without overlapping the top x-axis labels + * options.abrem.R (options.abrem): added option "mar", added some previously undocumented options, placed options in alphabetical order. + + +2014-03-23 Jurgen Symynck + + * upgraded abrem 0.1.10 to 0.1.12: + Nothing changed, using this version for debugging svn. + + +2013-12-09 Jurgen Symynck + + * upgraded abrem 0.1.8 to 0.1.10: + General man page cleanup and reorganization. Moved the detailed description of some options to the manpage of the most approporiate function. + added support for contour calculation. + removed abrem:::MLEw3p_secant.r; now calling the version in package debias. + removed abrem:::MRRw3pxy; now calling the version in package debias. + + * Abrem.R (Abrem): added fail and susp argument options + * options.abrem.R (options.abrem): added in.legend logical option for + * calculateSingleConf.R (calculateSingleConf): named likelihood ratio bounds to "lrb" instead of "lira" + * calculateSingleFit.R (calculateSingleFit): reduced length of file; more efficient code reusage. + +2013-11-02 Jurgen Symynck + + * upgraded abrem 0.1.7 to 0.1.8: + Added missing information (specifically about dist and method.fit) in several man pages. + General man page cleanup and reorganization. + Change dependencies to debias (>= 0.1.7) and pivotals (>= 0.1.9) + + * calculatesinglefit.R (CalculateSingleFit): replaced RBAw() with RBAbeta() \ No newline at end of file diff --git a/DESCRIPTION b/DESCRIPTION new file mode 100644 index 0000000..71feac9 --- /dev/null +++ b/DESCRIPTION @@ -0,0 +1,14 @@ +Package: abrem +Type: Package +Title: Abernethy Reliability Methods +Version: 0.1.21 +Date: August 23, 2014 +Author: Jurgen Symynck +Maintainer: Jurgen Symynck +Description: Implementation of reliability methods presented in "The New Weibull Handbook", Fifth Edition by Dr. Robert B. Abernethy. This package is dedicated to the control and view of models developed in dependant packages abremPivotals and abremDebias (under construction) implemented using the R object model. +Depends: abremPivotals (>= 0.2.8), debias (>= 0.1.9) +Suggests: MASS, boot +License: GPL-3 +URL: http://r-forge.r-project.org/projects/abernethy, http://www.openreliability.org +BugReports: email the author at +LazyLoad: yes diff --git a/NAMESPACE b/NAMESPACE new file mode 100644 index 0000000..6ff53d1 --- /dev/null +++ b/NAMESPACE @@ -0,0 +1,15 @@ +export( + Abrem, + abrem.fit, + abrem.conf, + params.to.ob, + options.abrem, + plot.abrem, + print.abrem, + contour.abrem) +S3method(plot, abrem) +S3method(print, abrem) +S3method(contour, abrem) +importFrom(abremPivotals) +importFrom(debias) + diff --git a/R/Abrem.R b/R/Abrem.R new file mode 100644 index 0000000..818c04a --- /dev/null +++ b/R/Abrem.R @@ -0,0 +1,154 @@ +# R package 'abrem' +# Abernethy Reliability Methods +# Implementations of lifetime data analysis methods described in +# 'The New Weibull Handbook, Fifth edition' by Dr. Robert B. Abernethy. +# August 2014, Jurgen Symynck +# Copyright 2014, Jurgen Symynck +# +# For more info, visit http://www.openreliability.org/ +# +# For the latest version of this file, check the Subversion repository at +# http://r-forge.r-project.org/projects/abernethy/ +# +# Disclaimer: +# The author is not affiliated with Dr. Abernethy or Wes Fulton - CEO of +# Fulton Findings(TM) and author of the software package SuperSMITH +#------------------------------------------------------------------------------- +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . +# +# +-----------------------------------+ +# | execute this software with R: | +# | http://www.r-project.org/ | +# +-----------------------------------+ + +Abrem <- function(x,...){ + # TODO: add code to convert a Surv() object into an abrem object + #arg <- list(...) + arg <- splitargs(...) + +# arg2 <- arg[!(names(arg) %in% names(options.abrem()))] + # extract the arguments that are NOT abrem options. + opa <- modifyList(options.abrem(), arg$opa) + ret <- list() + class(ret) <- "abrem" + timeorder <- c() +# lowest <- 1 + if(!missing(x)){ + ret$data <- NULL + if(is.vector(x) || is.numeric(x)){ + # TODO: the effects of adding is.numeric(x) are not tested thoroughly yet... + if(opa$verbosity >= 2)message( + 'Abrem: Argument \"x\" is a (numeric) vector of (life-)time observations...') + if(any(is.na(x))) timeorder <- 1:length(x) + else timeorder <- order(x) + # the above is to prevent ordering attempts when NA values are + # present in the lifetime observation vector. + # having NA values implies that the data must be ordered. + ret$data <- data.frame(time=x[timeorder],event=1) + } + if(is.data.frame(x)){ + if(!is.null(x$time) && !is.null(x$event)){ + if(opa$verbosity >= 2)message( + 'Abrem: Argument \"x\" is a dataframe with $time and $event ', + 'columns...') + if(any(is.na(x$time))) timeorder <- 1:length(x$time) + else timeorder <- order(x$time) + ret$data <- as.data.frame(x[timeorder,]) +# ret$data$event <- 1 +# # temporarily set event vector to 1 + }else{ + stop(': Argument \"x\" is missing $time and/or ", + "$event columns...') + } + } + }else{ + ti <- c(arg$rem$time,arg$rem$fail) + if(xor(!is.null(arg$rem$time), !is.null(arg$rem$fail))){ + if(is.vector(ti)){ + if(opa$verbosity >= 2)message( + 'Abrem: Argument \"time\" or \"fail\" is vector of complete (life-)time observations...') + if(any(is.na(ti))) timeorder <- 1:length(arg$rem$time) + else timeorder <- order(ti) + ret$data <- data.frame(time=ti[timeorder],event=1) + }else{stop('Argument \"time\" or fail\" must be vector.')} + } + if(!is.null(arg$rem$susp)){ + if(is.vector(arg$rem$susp)){ + if(opa$verbosity >= 2)message( + 'Abrem: Argument \"susp\" is vector of right-censored (suspended) (life-)time observations...') + timeorder <- order(c(ti,arg$rem$susp)) + ret$data <- data.frame(time=c(ti,arg$rem$susp)[timeorder], + event=c(rep(1,length(ti)),rep(0,length(arg$rem$susp)))[timeorder]) + }else{stop('Argument \"susp\" must be a vector.')} + } + #else{stop("No (life-)time observations were provided.")} + } + + ### setting the event vector correctly ### + if(!is.null(arg$rem$event) && !is.null(ret$data)){ + if(is.vector(arg$rem$event)){ + if(opa$verbosity >= 2)message( + 'Abrem: Argument \"event\" is event vector...') + ret$data$event <- arg$rem$event[timeorder] + } + } + addpppcolumns <- function(ppos){ + if(opa$verbosity >= 2)message(paste0( + 'Abrem: Adding ',ppos,' ranks to (life-)time observations...')) + ret$data <<- cbind(ret$data,ppp=NA) + if(any(is.na(ret$data$time))){ + # experimental code, in combination with support in + # abremPivotals::gePPP for event vector arguments + + ret$data[ret$data$event==1,'ppp'] <<- + abremPivotals::getPPP( + x=ret$data$event, + ppos=ppos)$ppp + }else{ + ret$data[ret$data$event==1,'ppp'] <<- + abremPivotals::getPPP( + x=ret$data$time[ret$data$event==1], + s=ret$data$time[ret$data$event==0], + ppos=ppos)$ppp + } + colnames(ret$data) <<- c( + colnames(ret$data)[-ncol(ret$data)], + paste0("ppp.",ppos)) + # renaming the added column to include the type of ranking + } + + if(any( + c("benard","beta","mean","km","hazen","blom") %in% tolower(opa$ppos))){ + do.call(addpppcolumns,as.list(tolower(opa$ppos))) + } + if(any( + c("kaplan-meier","kaplanmeier","kaplan_meier","kaplan.meier") %in% + tolower(opa$ppos))){ + do.call(addpppcolumns,list("km")) + } + ret$n <- length(ret$data$time) + # TODO: this assumes that any NA time (in any present + # in the time column is there for a good reason: + # accompanied with a censoring indicator (0 or FALSE) + # TODO: check if the above code is still valid! + # this feature must be researched in combination wirh abremPivotals:lslr() + ret$fail <- sum(ret$data$event) + # TODO: Warning; this only works when event can take values of 0 and one! + # This can be solved when events are changed to factors, as they really are + ret$susp <- ret$n-ret$fail + ret$options <- opa + # always store a full copy of the options.abrem structure here + ret + # TODO: check what to do with the automatically added row names that are sometimes out of order +} \ No newline at end of file diff --git a/R/F0inv.R b/R/F0inv.R new file mode 100644 index 0000000..0d6c831 --- /dev/null +++ b/R/F0inv.R @@ -0,0 +1,50 @@ +# R package 'abrem' +# Abernethy Reliability Methods +# Implementations of lifetime data analysis methods described in +# 'The New Weibull Handbook, Fifth edition' by Dr. Robert B. Abernethy. +# August 2014, Jurgen Symynck +# Copyright 2014, Jurgen Symynck +# +# For more info, visit http://www.openreliability.org/ +# +# For the latest version of this file, check the Subversion repository at +# http://r-forge.r-project.org/projects/abernethy/ +# +# Disclaimer: +# The author is not affiliated with Dr. Abernethy or Wes Fulton - CEO of +# Fulton Findings(TM) and author of the software package SuperSMITH +#------------------------------------------------------------------------------- +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . +# +# +-----------------------------------+ +# | execute this software with R: | +# | http://www.r-project.org/ | +# +-----------------------------------+ + +F0 <- function(q) + 1-exp(-exp(q)) + +F0inv <- function(p,log="x"){ + # transformation function to plot its argument + # on the y-axis of the Weibull plot. This transformation function + # lets the Weibull curve appear as a straight line on the weibull paper + # + # This is also the inverse Cumulative Distribution function of the + # standardized Weibull plot with beta=eta=1 + # comparing both implementationss of F0inv() with + # system.time() does not show any significant difference + # log(log(1/(1-p)))} + if(log %in% c("x",""))ret <- log(qweibull(p,1,1)) else ret <- qlnorm(p,0,1) + ret +} \ No newline at end of file diff --git a/R/abrem.conf.R b/R/abrem.conf.R new file mode 100644 index 0000000..c4ff38f --- /dev/null +++ b/R/abrem.conf.R @@ -0,0 +1,59 @@ +# R package 'abrem' +# Abernethy Reliability Methods +# Implementations of lifetime data analysis methods described in +# 'The New Weibull Handbook, Fifth edition' by Dr. Robert B. Abernethy. +# August 2014, Jurgen Symynck +# Copyright 2014, Jurgen Symynck +# +# For more info, visit http://www.openreliability.org/ +# +# For the latest version of this file, check the Subversion repository at +# http://r-forge.r-project.org/projects/abernethy/ +# +# Disclaimer: +# The author is not affiliated with Dr. Abernethy or Wes Fulton - CEO of +# Fulton Findings(TM) and author of the software package SuperSMITH +#------------------------------------------------------------------------------- +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . +# +# +-----------------------------------+ +# | execute this software with R: | +# | http://www.r-project.org/ | +# +-----------------------------------+ + +abrem.conf <- function(x,which="all",...){ + # x is a single Abrem or a list of Abrem objects +# supported_blifeconf <- c("mcpivotals","bbb") + if(missing(x)){ + stop('Argument \"x\" is missing.') + }else{ + if(identical(class(x),"abrem")) x <- list(x) + if(!all(sapply(x,function(x)identical(class(x),"abrem")))){ + stop('Argument \"x\" is not of class \"abrem\" or ", + "a list of \"abrem\" objects.') + } + dr <- findMaxDataRange(x,0) + # for reasonably large confidence bounds + calculateConfsInAbrem <- function(abrem){ + if(!is.null(abrem$fit)){ + abrem$fit <- lapply(abrem$fit,calculateSingleConf, + opadata=abrem$options,datarange=dr,...) + # TODO: add error detection + } + abrem + } + abremlist <- lapply(x,calculateConfsInAbrem) + } + if(length(abremlist)==1) abremlist[[1]] else abremlist +} diff --git a/R/abrem.fit.R b/R/abrem.fit.R new file mode 100644 index 0000000..fd2c59b --- /dev/null +++ b/R/abrem.fit.R @@ -0,0 +1,98 @@ +# R package 'abrem' +# Abernethy Reliability Methods +# Implementations of lifetime data analysis methods described in +# 'The New Weibull Handbook, Fifth edition' by Dr. Robert B. Abernethy. +# August 2014, Jurgen Symynck +# Copyright 2014, Jurgen Symynck +# +# For more info, visit http://www.openreliability.org/ +# +# For the latest version of this file, check the Subversion repository at +# http://r-forge.r-project.org/projects/abernethy/ +# +# Disclaimer: +# The author is not affiliated with Dr. Abernethy or Wes Fulton - CEO of +# Fulton Findings(TM) and author of the software package SuperSMITH +#------------------------------------------------------------------------------- +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . +# +# +-----------------------------------+ +# | execute this software with R: | +# | http://www.r-project.org/ | +# +-----------------------------------+ + +abrem.fit <- function(x,...){ + # x is a single Abrem or a list of Abrem objects + supported_dist <- c( + "weibull","weibull2p","weibull3p", + "lognormal","lognormal2p","lognormal3p") + supported_fit <- c("rr","rr2","mle","mle2","mle3","mle-rba","mle2-rba","mle3-rba") + if(missing(x)){ + stop("Argument \"x\" is missing.") + }else{ + if(identical(class(x),"abrem")) x <- list(x) + if(!all(sapply(x,function(x)identical(class(x),"abrem")))){ + stop("\"x\" argument is not of class \"abrem\" or ", + "a list of \"abrem\" objects.") + } + # from here on, x is a list of one or more abrem objects + if(!is.null(x[[1]]$options)){ + opa <- x[[1]]$options + # use the options of the first abrem object in the list + }else stop("Abrem object has not options set.") + opa <- modifyList(opa, list(...)) + if(is.null(opa$dist)){ + if(opa$verbosity >= 1)message("abrem.fit : ", + ": Target distribution defaults to weibull2p.") + opa$dist <- "weibull2p" + }else{ + if(length(opa$dist)>1) + stop("Too many target distributions supplied.") + if(!any(tolower(opa$dist) %in% supported_dist)){ + if(any(tolower(opa$dist) %in% c("gumbel"))){ + stop("Gumbel is not yet a supported target fit method.") + }else{ + stop(paste0(opa$dist," is not a supported target fit method.")) + } + }else{ + if(is.null(opa$method.fit)){ + if(opa$verbosity >= 1)message("abrem.fit : ", + ': Fit method defaults to \"rr\", \"xony\".') + opa$method.fit <- c("rr","xony") + }else{ + fits <- length(which(opa$method.fit %in% supported_fit)) + if(fits > 1){ + stop("Only one fit method should be supplied.") + }else{ + if(any(c("rr","rr2") %in% tolower(opa$method.fit))){ + if(!any(c("xony","yonx") %in% + tolower(opa$method.fit))){ + if(opa$verbosity >= 1){ + message("abrem.fit : ", + ': Fit method \"rr\" defaults to \"xony\"') + opa$method.fit <- c(opa$method.fit,"xony") + } + } + } + x <- lapply(x,calculateSingleFit,...) + # TODO: only one object or list of abrem objects? + } + } + } + } + } + if(length(x)==1) x[[1]] else x + # return list of abrem objects when argument x was a list + # otherwise return single abrem object +} \ No newline at end of file diff --git a/R/bbb.R b/R/bbb.R new file mode 100644 index 0000000..2ba4971 --- /dev/null +++ b/R/bbb.R @@ -0,0 +1,46 @@ +# R package 'abrem' +# Abernethy Reliability Methods +# Implementations of lifetime data analysis methods described in +# 'The New Weibull Handbook, Fifth edition' by Dr. Robert B. Abernethy. +# August 2014, Jurgen Symynck +# Copyright 2014, Jurgen Symynck +# +# For more info, visit http://www.openreliability.org/ +# +# For the latest version of this file, check the Subversion repository at +# http://r-forge.r-project.org/projects/abernethy/ +# +# Disclaimer: +# The author is not affiliated with Dr. Abernethy or Wes Fulton - CEO of +# Fulton Findings(TM) and author of the software package SuperSMITH +#------------------------------------------------------------------------------- +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . +# +# +-----------------------------------+ +# | execute this software with R: | +# | http://www.r-project.org/ | +# +-----------------------------------+ + +bbb <- function(j,f,CL,beta,eta){ + # function to calculate Beta Binomial Confidence Bounds for B-lives. + # j : rank of failure + # f : number of failures (is NOT the same as sample size!) see + # 'suspended items' fotr more info + # CB : Confidence Bound, Confidence Limit + # beta : Weibull slope or scale parameter + # eta : Weibull shape paramater + # see "The new Weibull handbook, fifth edition" p. 7-3 + # see "The new Weibull handbook, fifth edition" Appendix I + # see also MS. Excel's and GNUmeric's BETAINV() function + eta*(log(1/(1-qbeta(CL,j,f-j+1))))^(1/beta)} diff --git a/R/blifestring.R b/R/blifestring.R new file mode 100644 index 0000000..8348337 --- /dev/null +++ b/R/blifestring.R @@ -0,0 +1,84 @@ +# R package 'abrem' +# Abernethy Reliability Methods +# Implementations of lifetime data analysis methods described in +# 'The New Weibull Handbook, Fifth edition' by Dr. Robert B. Abernethy. +# August 2014, Jurgen Symynck +# Copyright 2014, Jurgen Symynck +# +# For more info, visit http://www.openreliability.org/ +# +# For the latest version of this file, check the Subversion repository at +# http://r-forge.r-project.org/projects/abernethy/ +# +# Disclaimer: +# The author is not affiliated with Dr. Abernethy or Wes Fulton - CEO of +# Fulton Findings(TM) and author of the software package SuperSMITH +#------------------------------------------------------------------------------- +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . +# +# +-----------------------------------+ +# | execute this software with R: | +# | http://www.r-project.org/ | +# +-----------------------------------+ + +Blifestring <- function(B,blicon,signif,...){ + # This functions creates a string for displaying the B-lives in the plot's + # legend. missing input data result in an "NA". For example, the output + # string could look like: + # "B10 = 9.86 | 50.13 | 103.4" + # or + # "B1 = 9.86 | 50.13 | NA" + si <- function(number) + if(!is.null(number))signif(number,signif) + else NA + # shorthand writing of the signif() function + qfun <- function(B,...){ + args <- as.list(unlist(...)) + + ret <- NULL + if(!is.null(args$beta) && !is.null(args$eta)){ + # the fit type was weibull + ret <- qweibull(B,args$beta,args$eta) + if(!is.null(args$t0)){ + # the fit type was weibull3p + ret <- ret+args$t0 + + } + } + if(!is.null(args$mulog) && !is.null(args$sigmalog)){ + # the fit type was lognormal + ret <- qlnorm(B,args$mulog,args$sigmalog) + } + if(!is.null(args$rate)){ + # the fit type was exponential + ret <- qexp(B,args$rate) + } + ret + } + id <- function(x,y)isTRUE(all.equal(x,y)) + c1 <- is.null(blicon$bounds) || is.null(blicon$bounds$lower) + if(!c1) lo <- si(subset(blicon$bounds, + sapply(blicon$bounds$unrel,id,B),lower)) + c2 <- is.null(blicon$bounds) || is.null(blicon$bounds$upper) + if(!c2) up <- si(subset(blicon$bounds, + sapply(blicon$bounds$unrel,id,B),upper)) + ret <- paste(sep = ""," B",signif(100*B)," = ", + ifelse(c1, + "NA",lo), + " | ",si(qfun(B,...)), + " | ",ifelse(c2, + "NA",up)) + ret +# NA +} diff --git a/R/bsll.r b/R/bsll.r new file mode 100644 index 0000000..4bfb5c3 --- /dev/null +++ b/R/bsll.r @@ -0,0 +1,53 @@ +# R package 'abrem' +# Abernethy Reliability Methods +# Implementations of lifetime data analysis methods described in +# 'The New Weibull Handbook, Fifth edition' by Dr. Robert B. Abernethy. +# August 2014, Jurgen Symynck +# Copyright 2014, Jurgen Symynck +# +# For more info, visit http://www.openreliability.org/ +# +# For the latest version of this file, check the Subversion repository at +# http://r-forge.r-project.org/projects/abernethy/ +# +# Disclaimer: +# The author is not affiliated with Dr. Abernethy or Wes Fulton - CEO of +# Fulton Findings(TM) and author of the software package SuperSMITH +#------------------------------------------------------------------------------- +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . +# +# +-----------------------------------+ +# | execute this software with R: | +# | http://www.r-project.org/ | +# +-----------------------------------+ + +bsll <- function(...){ + arg <- list(...) + leline <- list( + legend= NA, + lty= NA, + lwd= NA, + pch= NA, + col= NA) + modifyList(leline,arg) +# leline <- list( +# legend= <- ifelse(is.null(arg$legend),NA,arg$legend) +# title= <- ifelse(is.null(arg$title),NA,arg$title) +# cex= <- ifelse(is.null(arg$cex),NA,arg$cex) +# bg= <- ifelse(is.null(arg$bg),NA,arg$bg) +# lty= <- ifelse(is.null(arg$lty),NA,arg$lty) +# lwd= <- ifelse(is.null(arg$lwd),NA,arg$lwd) +# pch= <- ifelse(is.null(arg$pch),NA,arg$pch) +# col= <- ifelse(is.null(arg$col),NA,arg$col) +} diff --git a/R/buildListOfLegends.R b/R/buildListOfLegends.R new file mode 100644 index 0000000..0bdcdd9 --- /dev/null +++ b/R/buildListOfLegends.R @@ -0,0 +1,56 @@ +# R package 'abrem' +# Abernethy Reliability Methods +# Implementations of lifetime data analysis methods described in +# 'The New Weibull Handbook, Fifth edition' by Dr. Robert B. Abernethy. +# August 2014, Jurgen Symynck +# Copyright 2014, Jurgen Symynck +# +# For more info, visit http://www.openreliability.org/ +# +# For the latest version of this file, check the Subversion repository at +# http://r-forge.r-project.org/projects/abernethy/ +# +# Disclaimer: +# The author is not affiliated with Dr. Abernethy or Wes Fulton - CEO of +# Fulton Findings(TM) and author of the software package SuperSMITH +#------------------------------------------------------------------------------- +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . +# +# +-----------------------------------+ +# | execute this software with R: | +# | http://www.r-project.org/ | +# +-----------------------------------+ + +buildListOfLegends <- function(abrem,v,isplotlegend,...){ + ret <- NULL + if(!is.null(abrem$fit) && any(sapply(abrem$fit,function(fi)!is.null(fi)))){ + # TODO: +# if(!is.null(abrem$fit)){ + # check if any non-NULL list holds only NULL items + # this is needed for dealing with failed fit attempts + # that currently take the form of + # abrem$fit[i] <- list(NULL) +# ret <- unlist(lapply(x$fit,buildSingleFitLegend, +# opadata=x$options,...),FALSE) + ret <- lapply(abrem$fit,buildSingleFitLegend, + opadata=abrem$options,...) + }else{ + if(abrem$options$is.plot.legend && isplotlegend){ + ret <- list(buildSingleDataLegend(abrem,opadata=abrem$options,...)) + if(v >= 1)message( + "buildListOfLegends: This Abrem object contains no fits.") + } + } + ret +} \ No newline at end of file diff --git a/R/buildSingleDataLegend.R b/R/buildSingleDataLegend.R new file mode 100644 index 0000000..6f3c088 --- /dev/null +++ b/R/buildSingleDataLegend.R @@ -0,0 +1,66 @@ +# R package 'abrem' +# Abernethy Reliability Methods +# Implementations of lifetime data analysis methods described in +# 'The New Weibull Handbook, Fifth edition' by Dr. Robert B. Abernethy. +# August 2014, Jurgen Symynck +# Copyright 2014, Jurgen Symynck +# +# For more info, visit http://www.openreliability.org/ +# +# For the latest version of this file, check the Subversion repository at +# http://r-forge.r-project.org/projects/abernethy/ +# +# Disclaimer: +# The author is not affiliated with Dr. Abernethy or Wes Fulton - CEO of +# Fulton Findings(TM) and author of the software package SuperSMITH +#------------------------------------------------------------------------------- +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . +# +# +-----------------------------------+ +# | execute this software with R: | +# | http://www.r-project.org/ | +# +-----------------------------------+ + +buildSingleDataLegend <- function(x,opadata,...){ + arg <- list(...) + si <- function(number)signif(number,opadata$signif) + li <- list() + if(opadata$in.legend){ + li[[10]] <- bsll(legend=paste0("ppos = ",opadata$ppos[1]), + col=opadata$col,pch=opadata$pch,lwd=opadata$lwd.points) + li[[15]] <- bsll(legend=paste0("n (fail | susp.) = ",x$n, + " (",x$fail," | ",x$susp,")")) + } + removeBadLegendEntries <- function(e){ + if(!is.null(e))!is.na(e$legend) else FALSE + } + if(length(li)>0)li <- li[sapply(li,removeBadLegendEntries)] + else li <- "" + # remove list items where the legend text = NA + fu <- function(x,i){if(i %in% names(x))x[[i]]} + fu2 <- function(i,x){lapply(x,fu,i=i)} + items <- c("legend","lty","lwd","pch","col") + le <- lapply(items,fu2,li) + names(le) <- items + if(identical(label <- opadata$label,""))label <- NULL + le$rect <- legend( + "bottomright", + legend=le$legend, + title=label, + cex = opadata$legend.text.size, + plot=FALSE)$rect + le$label <- opadata$label + le$legend.text.size <- opadata$legend.text.size + le +} diff --git a/R/buildSingleFitLegend.R b/R/buildSingleFitLegend.R new file mode 100644 index 0000000..774f926 --- /dev/null +++ b/R/buildSingleFitLegend.R @@ -0,0 +1,132 @@ +# R package 'abrem' +# Abernethy Reliability Methods +# Implementations of lifetime data analysis methods described in +# 'The New Weibull Handbook, Fifth edition' by Dr. Robert B. Abernethy. +# August 2014, Jurgen Symynck +# Copyright 2014, Jurgen Symynck +# +# For more info, visit http://www.openreliability.org/ +# +# For the latest version of this file, check the Subversion repository at +# http://r-forge.r-project.org/projects/abernethy/ +# +# Disclaimer: +# The author is not affiliated with Dr. Abernethy or Wes Fulton - CEO of +# Fulton Findings(TM) and author of the software package SuperSMITH +#------------------------------------------------------------------------------- +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . +# +# +-----------------------------------+ +# | execute this software with R: | +# | http://www.r-project.org/ | +# +-----------------------------------+ + +buildSingleFitLegend <- function(fit,opadata,...){ + arg <- list(...) + if(!is.null(fit$options)){ + opafit <- modifyList(opadata,fit$options) + }else{opafit <- opadata} + opafit <- modifyList(opafit,list(...)) + t0 <- NULL + le <- NULL + + if(opafit$is.plot.legend){ + if(is.logical(opafit$threshold))if(opafit$threshold){ + if(is.logical(opadata$threshold)){if(opadata$threshold) + warning("opafit$threshold and opadata$threshold are logical values but numeric values were expected. Proceeding...") + }else{ + # reuse the t0 value from the data level + t0 <- opadata$threshold + } + } + if(is.numeric(opafit$threshold))t0 <- opafit$threshold + si <- function(number)signif(number,opafit$signif) + # shorter writing form for signif() + li <- list() + if(opadata$in.legend){ + + li[[10]] <- bsll(legend=paste0("ppos = ",opafit$ppos[1],ifelse(is.null(t0),"",paste0(" (t0 = ",si(t0),")"))), + col=opadata$col,pch=opadata$pch,lwd=opadata$lwd.points) + li[[15]] <- bsll(legend=paste0("n (fail | susp.) = ",fit$n, + " (",fit$fail," | ",fit$susp,")")) + } + if(opafit$in.legend){ + li[[20]] <- bsll(legend = paste0(fit$options$dist," (", + paste0(fit$options$method.fit,collapse=", "),")"), + col=opafit$col,lwd=opafit$lwd,lty=opafit$lty) +# li[[30]] <- bsll(legend=ifelse(is.null(fit$rate),NA, +# paste0("rate = ",si(fit$rate)))) + li[[40]] <- bsll(legend=ifelse(is.null(fit$mulog),NA, + paste0("mu(log) = ",si(exp(fit$mulog))," (", + si(fit$mulog),")"))) + li[[50]] <- bsll(legend=ifelse(is.null(fit$sigmalog),NA, + paste0("sigma(log) = ",si(exp(fit$sigma))," (", + si(fit$sigmalog),")"))) + li[[60]] <- bsll(legend=ifelse(is.null(fit$beta),NA, + paste0("beta = ",si(fit$beta)))) + li[[70]] <- bsll(legend=ifelse(is.null(fit$eta),NA, + paste0("eta = ",si(fit$eta)))) + li[[80]] <- bsll(legend=ifelse(is.null(fit$t0),NA, + paste0("t0 = ",si(fit$t0)))) + if(!is.null(fit$gof) && opafit$in.legend.gof){ + if(!is.null(fit$gof$r2)){ + if(!is.null(fit$gof$CCC2)){ + li[[100]] <- bsll(legend=paste0("r^2 | CCC^2 = ", + si(fit$gof$r2)," | ",si(fit$gof$CCC2), + ifelse(fit$gof$r2>=fit$gof$CCC2," (good)"," (BAD)"))) + }else{ + li[[100]] <- bsll(legend=paste0("r^2 = ",si(fit$gof$r2))) + } + } + if(!is.null(fit$gof$loglik)){ + li[[110]] <- bsll(legend=paste0("loglik = ",si(fit$gof$loglik))) + } + li[[120]] <- bsll( + legend=ifelse(is.null(fit$gof$AbPval),NA, + paste0("AbPval = ",si(fit$gof$AbPval)))) + #," (S=",ifelse(is.null(fit$gof$S),"NA",fit$gof$S),")"))) + # AbPval is not created by MC, so S is not used and needed here + } + } + #leconfpos <- length(na.omit(unlist(li))) + 1 + # where displaying confidence info begins + leconf <- legendConf(fit,"blives",opadata=opadata,...) + if(!is.null(leconf))li[[130]] <- bsll(legend="") + li <- c(li,leconf) + removeBadLegendEntries <- function(e){ + if(!is.null(e))!is.na(e$legend) else FALSE + } + if(length(li)>0)li <- li[sapply(li,removeBadLegendEntries)] + else li <- "" + # remove list items where the legend text = NA + fu <- function(x,i){if(i %in% names(x))x[[i]]} + fu2 <- function(i,x){lapply(x,fu,i=i)} + items <- c("legend","lty","lwd","pch","col") + le <- lapply(items,fu2,li) + names(le) <- items + if(identical(label <- opafit$label,""))label <- NULL + le$rect <- legend( + "bottomright", + # "topright", + legend=le$legend, + title=label, + cex = opafit$legend.text.size, + # inset=0.1, + # merge = TRUE, + plot=FALSE)$rect + le$label <- opafit$label + le$legend.text.size <- opafit$legend.text.size + } + le +} diff --git a/R/calculateSingleConf.R b/R/calculateSingleConf.R new file mode 100644 index 0000000..77e84b4 --- /dev/null +++ b/R/calculateSingleConf.R @@ -0,0 +1,425 @@ +# R package 'abrem' +# Abernethy Reliability Methods +# Implementations of lifetime data analysis methods described in +# 'The New Weibull Handbook, Fifth edition' by Dr. Robert B. Abernethy. +# August 2014, Jurgen Symynck +# Copyright 2014, Jurgen Symynck +# +# For more info, visit http://www.openreliability.org/ +# +# For the latest version of this file, check the Subversion repository at +# http://r-forge.r-project.org/projects/abernethy/ +# +# Disclaimer: +# The author is not affiliated with Dr. Abernethy or Wes Fulton - CEO of +# Fulton Findings(TM) and author of the software package SuperSMITH +#------------------------------------------------------------------------------- +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . +# +# +-----------------------------------+ +# | execute this software with R: | +# | http://www.r-project.org/ | +# +-----------------------------------+ + +calculateSingleConf <- function(fit,opadata,datarange,...){ + # fit is a single fit + arg <- list(...) + if(!is.null(list(...)$threshold)) + message("calculateSingleConf: Currently, passing the \'threshold\' argument to abrem.conf is not supported. Proceeding...") + if(missing(fit)){ + stop("Argument \"fit\" is missing.") + }else{ + if(!is.null(fit) && any(sapply(fit,function(fi)!is.null(fi)))){ + if(!is.null(fit$options)) opafit <- modifyList(opadata,fit$options) + opafit$importfile <- NULL + # never use the importfile from abrem.fit; only use it when + # explicitly supplied as a function argument + opaconf <- modifyList(opafit,arg) + if(!is.null(fit$options$dist)){ +# if(opaconf$verbosity >= 1) +# message("calculateSingleConf: Found weibull 2P distribution.") + if("blives" %in% tolower(opaconf$conf.what)){ + if(opaconf$verbosity >= 1) + message("calculateSingleConf: Calculating B-lives confidence bounds.") + mini <- min(c(opaconf$ylim[1]/10,datarange$yrange[1]/10),0.001) + maxi <- max(c((1-(1-opaconf$ylim[2])/10), + (1-(1-datarange$yrange[2])/10),0.999)) + if(opaconf$unrel.n==0){ + # assume that the user supplies enough values in the unrel + # vector to display confidence bounds + # this feature is useful in combination with unit testing + # (testthat) of sumerSMITH iported plotting reports + }else{ + unrel <- c(F0(seq(F0inv(mini),F0inv(maxi), + # TODO: this isn't right... + # unrel <- c(F0(seq(F0inv(1e-3), + # F0inv(0.999), + length.out=opaconf$unrel.n - + length(opaconf$unrel+2))), + opaconf$unrel,0.5,F0(0)) + # effectively ignoring any ylim + # setting per fit. + unrel <- unique(signif(unrel[order(unrel)])) + # signif() needed for eliminating + # identical looking unreliability + # levels that differ only at place far + # from the decimal point + # unrel <- c(F0(seq(par('usr')[3],par('usr')[4], + } + if(is.null(fit$conf)) fit$conf <- list() + atLeastOneBLifeConf <- FALSE + if(is.null(fit$conf$blives)){ + if(opaconf$verbosity >= 2)message( + "calculateSingleConf: Creating the first ", + "B-life confidence calculation in the fit...") + i <- 1 + fit$conf$blives <- list() + }else{ + if(opaconf$verbosity >= 2)message( + "calculateSingleConf: Appending a new ", + "B-life confidence calculation to the fit...") + i <- length(fit$conf$blives)+1 + } + fit$conf$blives[[i]] <- list() + if(!is.null(opaconf$importfile)){ + # ____ __ __ ___ _____ _ _ + # ___ _ _ _ __ ___ _ __/ ___|| \/ |_ _|_ _| | | | + # / __| | | | '_ \ / _ \ '__\___ \| |\/| || | | | | |_| | + # \__ \ |_| | |_) | __/ | ___) | | | || | | | | _ | + # |___/\__,_| .__/ \___|_| |____/|_| |_|___| |_| |_| |_| + # |_| + if(opaconf$verbosity >= 1)message("calculateSingleConf : ", + "importing confidence bounds from superSMITH report file:\n",opaconf$importfile) + try(fi <- file(opaconf$importfile)) + if(!is.null(fi)){ + try(re <- readLines(fi)) + if(!is.null(re)){ + atLeastOneBLifeConf <- TRUE + fit$conf$blives[[i]]$type <- "superSMITH" + fit$conf$blives[[i]]$source <- re + extract <- function(string) + na.omit(as.numeric(unlist(strsplit(gsub(",",".",string)," ")))) + bounds <- data.frame(do.call("rbind",lapply(re[12:length(re)],extract))) + bounds[,1] <- bounds[,1]/100 + names(bounds) <- c("unrel","lower","datum", "upper") + fit$conf$blives[[i]]$bounds <- bounds + } + close(fi) + } + op <- unique(c(names(opafit),names(opaconf))) + if(length(li <- opaconf[sapply(op,function(y){ + !identical(opafit[[y]], opaconf[[y]])})]) > 0){ + fit$conf$blives[[i]]$options <- li + } + return(fit) + }else{ + if("bbb" %in% tolower(opaconf$method.conf.blives)){ + # ____ ____ ____ + # | __ )| __ )| __ ) + # | _ \| _ \| _ \ + # | |_) | |_) | |_) | + # |____/|____/|____/ + + # TODO: bbb for mle? + # TODO: ppp calculation on the fly or taken from x$fit[[i]]$data? + # TODO: what if no ppp are in x$fit[[i]]$data? + # TODO: bbb for lognormal? + + if(opaconf$verbosity >= 1)message( + "calculateSingleConf: Calculating bbb confidence bounds...") + fit$conf$blives[[i]]$type <- "bbb" + fit$conf$blives[[i]]$cl <- opaconf$cl + fit$conf$blives[[i]]$sides <- opaconf$conf.blives.sides + + ### calculate adjusted ranks, just for these BB bounds + sx <- fit$data + if(is.null(fit$data[['ppp',exact=FALSE]])){ + # no ranks available, likely because the fit was mle + if(opaconf$verbosity >= 1) + message('calculateSingleConf: creating ranks for calculating \"bbb\" bounds...') + # ppp <- .Call(ty,fit$data$event, PACKAGE= "pivotals") + ppp <- abremPivotals::getPPP( + x=fit$data$time[fit$data$event==1], + s=fit$data$time[fit$data$event==0], + ppos=opaconf$ppos)$ppp + # TODO: ppos=opaconf$ppos is not compatible with the intended feature that multiple adjuste rank methods can be supplied! + sx$ppp[sx$event] <- ppp + # this assumes fit$data and thus sx is ordered + # TODO: shouldn't I use the ranks that are part of the + # data since they will ALWAYS be part of + # an abrem object? + # The argument can be made that for MLE bbb should also be calculated + # + # TODO: makes sense to have BBB bounds with mle? + } + sx <- sx[order(sx[['ppp',exact=FALSE]]),] + # order data according to rank + sx <- cbind(sx,arank=NA) + sx$rrank <- (fit$n+1-order(sx[['ppp',exact=FALSE]])) + # TODO: does order() completely replace x$ppp? (NA?) + # reverse rank order + # TODO: keep the rrank and arank in fit$data or discard? + parank <- 0 + for (j in 1:fit$n){ + if(!sx$event[j] || is.null(sx$event)){ + sx$arank[j] <- NA + }else{ + sx$arank[j] <- (sx$rrank[j]*parank + fit$n +1)/(sx$rrank[j]+1) + parank <- sx$arank[j] + } + # adjusted_rank = + # (reversed rank * previous adj. rank + n + 1)/(reversed rank + 1) + # see "The new Weibull handbook, fifth edition" p. 2-7, formula 2-5 + } + da <- data.frame( + unrel= sx[['ppp',exact=FALSE]], + lower= bbb(sx$arank,fit$n,(1-opaconf$cl)/2,fit$beta,fit$eta), + upper= bbb(sx$arank,fit$n,1-(1-opaconf$cl)/2,fit$beta,fit$eta)) + lo <- approxfun(F0inv(sx[['ppp',exact=FALSE]]),log(da$lower)) + up <- approxfun(F0inv(sx[['ppp',exact=FALSE]]),log(da$upper)) + bl <- F0inv(unrel) + da <- rbind(da,data.frame(unrel=unrel,lower=exp(lo(bl)),upper=exp(up(bl)))) + da <- da[order(da$unrel),] + da <- da[!duplicated(da$unrel),] + fit$conf$blives[[i]]$bounds <- da + op <- unique(c(names(opafit),names(opaconf))) + # this is needed to add options from opafit into li that + # are NULL in opafit + # TODO: tolower() not needed? + if(length(li <- opaconf[sapply(op,function(y){ + !identical(opafit[[y]], opaconf[[y]])})]) > 0){ + fit$conf$blives[[i]]$options <- li + } + } + if("lrb" %in% tolower(opaconf$method.conf.blives)){ + # _ ____ ____ + # | | | _ \| __ ) + # | | | |_) | _ \ + # | |___| _ <| |_) | + # |_____|_| \_\____/ + + if(opaconf$verbosity >= 1)message( + "calculateSingleConf: Calculating Likelihood Ratio confidence bounds.") + fail <- fit$data$time[fit$data$event==1] + susp <- fit$data$time[fit$data$event==0] + # TODO: this implies only right censoring ! + fit$conf$blives[[i]] <- list() + fit$conf$blives[[i]]$type <- "lrb" + fit$conf$blives[[i]]$cl <- opaconf$cl + fit$conf$blives[[i]]$sides <- opaconf$conf.blives.sides + fit$conf$blives[[i]]$unrel <- opaconf$unrel + is_debias <- ifelse("mle" %in% tolower(fit$options$method.fit),FALSE,TRUE) + ret <- NULL + con <- NULL + try(con <- debias::MLEw2pContour(fail,susp,CL=opaconf$cl,debias=is_debias,show=FALSE)) + if(!is.null(con)){ + fit$conf$blives[[i]]$MLEXContour <- list() + fit$conf$blives[[i]]$MLEXContour[[1]] <- con + + retfit <- abrem.fit(Abrem(fail=fail,susp=susp),method.fit=ifelse(is_debias,"mle-rba","mle")) + fit$conf$blives[[i]]$MLEXContour$MLEpoint <- + data.frame(Eta=retfit$fit[[1]]$eta,Beta=retfit$fit[[1]]$beta) +# MLEpoint <- data.frame(Eta=fit$Eta,Beta=fit$Beta) +# fit$conf$blives[[i+1]] <- list() +# fit$conf$blives[[i+1]]$MLEXContour <- +# list(upper=MLEpoint,lower=MLEpoint,right=MLEpoint,left=MLEpoint) +# # TODO: this hack should produce a point at the MLE location + try(ret <- debias::MLEw2pBounds(fail,susp,Blives=unrel,MLEcontour=con,debias=is_debias,show=FALSE)) + # debias is true in all cases except for regulaer MLE + if(!is.null(ret)){ + atLeastOneBLifeConf <- TRUE + fit$conf$blives[[i]]$bounds <- cbind(unrel,exp(ret[,-1])) + names(fit$conf$blives[[i]]$bounds) <- tolower(names(fit$conf$blives[[i]]$bounds)) + # TODO: ask jacob to provide his dataframe in lowercase + } + } + op <- unique(c(names(opafit),names(opaconf))) + # this is needed to add options from opafit into li that + # are NULL in opafit + # TODO:tolower() not needed? + if(length(li <- opaconf[sapply(op,function(y){ + !identical(opafit[[y]], opaconf[[y]])})]) > 0){ + fit$conf$blives[[i]]$options <- li + } + } + if(any(c("mcpivotals","mcpivotal") %in% tolower(opaconf$method.conf.blives))){ + # _ _ _ + # _ __ ___ ___ _ __ (_)_ _____ | |_ __ _| |___ + # | '_ ` _ \ / __| '_ \| \ \ / / _ \| __/ _` | / __| + # | | | | | | (__| |_) | |\ V / (_) | || (_| | \__ \ + # |_| |_| |_|\___| .__/|_| \_/ \___/ \__\__,_|_|___/ + # |_| + + if(opaconf$verbosity >= 1)message( + "calculateSingleConf: Calculating Monte Carlo Pivotal confidence bounds.") +# i <- length(fit$conf$blives) +# if(length(fit$conf$blives[[i]]) > 0){ +# i <- i+1 +# } + if(any(c("weibull","weibull2p","weibull3p") %in% tolower(fit$options$dist))){ + dst <- "weibull" + dx <- params.to.ob(fit$options$dist,eta=1,beta=1,event=fit$data$event) + } + if(any(c("lognormal","lognormal2p","lognormal3p") %in% tolower(fit$options$dist))){ + dst <- "lognormal" + dx <- params.to.ob(fit$options$dist,mulog=0,sigmalog=1,event=fit$data$event) + } + r1 <- abrem.fit(Abrem(dx[dx$event==1,]),dist=fit$options$dist, + method.fit=fit$options$method.fit) + # TODO: what happens when the above are NULL? + fit$conf$blives[[i]] <- list() + fit$conf$blives[[i]]$type <- "mcpivotals" + fit$conf$blives[[i]]$S <- opaconf$S + fit$conf$blives[[i]]$seed <- opaconf$seed + fit$conf$blives[[i]]$rgen <- opaconf$rgen + fit$conf$blives[[i]]$cl <- opaconf$cl + fit$conf$blives[[i]]$sides <- opaconf$conf.blives.sides + fit$conf$blives[[i]]$unrel <- opaconf$unrel + ret <- NULL +# if(fit$susp != 0){ +# warning( +# "calculateSingleConf: Currently, MC Pivotal bounds for (heavily) censored data\n", +# "are still experimental.") +# } + if(is.null(fit$data[['ppp',exact=FALSE]])){ + message("calculateSingleConf: Currently, only rank regression is supported.") + }else{ +# try(ret <- .Call("pivotalMC", + if(dst=="weibull" && !is.null(r1$fit[[1]]$eta) && !is.null(r1$fit[[1]]$beta)){ + Scale <- r1$fit[[1]]$eta + Shape <- r1$fit[[1]]$beta + } + if(dst=="lognormal" && !is.null(r1$fit[[1]]$mulog) && !is.null(r1$fit[[1]]$sigmalog)){ + Scale <- r1$fit[[1]]$mulog + Shape <- r1$fit[[1]]$sigmalog + } + + try(ret <- abremPivotals::pivotalMC( + x=na.omit(fit$data), + dist=dst, + reg_method=ifelse("xony" %in% tolower(opaconf$method.fit),"xony","yonx"), + R2=0.0,CL=opaconf$cl, + # by passing R2=0.0, only the pivotals are returned in a matrix + unrel=unrel, + P1=Scale, + P2=Shape, + S=opaconf$S, + seed=sample.int(.Machine$integer.max,1), + ProgRpt=FALSE + )) +# ,PACKAGE = "abremPivotals")) + } + if(!is.null(ret)){ + atLeastOneBLifeConf <- TRUE + if(dst=="weibull")fit$conf$blives[[i]]$bounds <- + cbind(unrel,exp(log(fit$eta)+ ret/fit$beta)) + if(dst=="lognormal")fit$conf$blives[[i]]$bounds <- + cbind(unrel,exp(fit$mulog+ ret*fit$sigmalog)) + names(fit$conf$blives[[i]]$bounds) <- c("unrel","lower","datum", "upper") + op <- unique(c(names(opafit),names(opaconf))) + # this is needed to add options from opafit into li that + # are NULL in opafit + # TODO:tolower() not needed? + if(length(li <- opaconf[sapply(op,function(y){ + !identical(opafit[[y]], opaconf[[y]])})]) > 0){ + fit$conf$blives[[i]]$options <- li + } + }else{ + message("calculateSingleConf: Confidence calculation failed.") + fit$conf$blives[[i]] <- NULL + } + } + if(any(c("exp1","exp-1") %in% tolower(opaconf$method.conf.blives))){ + # experimental R based pivotals code + # _ + # _____ ___ __ / | + # / _ \ \/ / '_ \| | + # | __/> <| |_) | | + # \___/_/\_\ .__/|_| + # |_| + if(opaconf$verbosity >= 1)message( + "calculateSingleConf: Calculating EXP-1 confidence bounds.") + dx <- params.to.ob("weibull",beta=1,eta=1, + event=fit$data$event) + r1 <- abrem.fit(Abrem(dx[dx$event==1,]),dist=fit$options$dist, + method.fit=fit$options$method.fit) + # TODO: what happens when the above are NULL? + # TODO: no problems with NA? + fit$conf$blives[[i]] <- list() + fit$conf$blives[[i]]$type <- "exp1" + fit$conf$blives[[i]]$S <- opaconf$S + fit$conf$blives[[i]]$seed <- opaconf$seed + fit$conf$blives[[i]]$rgen <- opaconf$rgen + fit$conf$blives[[i]]$cl <- opaconf$cl + fit$conf$blives[[i]]$sides <- opaconf$conf.blives.sides + fit$conf$blives[[i]]$unrel <- opaconf$unrel + ret <- NULL + daevent <- fit$data$event + MCfun <- function(){ + d2 <- data.frame(time=NA,event=daevent) + d2[daevent == 1,'time'] <- sort( + rweibull(length(daevent[daevent==1]), + r1$fit[[1]]$beta,r1$fit[[1]]$eta)) + if(!is.null(ret))return( + c(u_hat=log(ret[1]),b_hat=1/ret[2])) + else stop() + } + piv <- as.data.frame(t(replicate(opaconf$S,MCfun()))) + wp <- log(qweibull(unrel,1,1)) + #wp <- F0inv(unrel) # identical as the above line + Zp <- function(wp)((piv$u_hat-wp)/piv$b_hat) + # calculate the pivotal quantities for each u_hat and b_hat... + piv <- cbind(piv,sapply(wp,Zp)) + # ... and add them to the dataframe + names(piv) <- c("u_hat","b_hat",signif(unrel)) + + fit$conf$blives[[i]]$bounds <- data.frame(unrel=unrel,row.names=unrel) + Tp <- function(Zp,cl)exp(log(fit$eta)-quantile(Zp,cl)/fit$beta) + fit$conf$blives[[i]]$bounds <- + cbind(fit$conf$blives[[i]]$bounds, + lower =sapply(piv[,c(-1,-2)],Tp,1-(1-opaconf$cl)/2), + datum =sapply(piv[,c(-1,-2)],Tp,0.5), + upper =sapply(piv[,c(-1,-2)],Tp,(1-opaconf$cl)/2)) + op <- unique(c(names(opafit),names(opaconf))) + # this is needed to add options from opafit into li that + # are NULL in opafit + # TODO:tolower() not needed? + if(length(li <- opaconf[sapply(op,function(y){ + !identical(opafit[[y]], opaconf[[y]])})]) > 0){ + fit$conf$blives[[i]]$options <- li + } + } + } + } +#if("weibull3p" %in% tolower(fit$options$dist)){ +# if(is.null(fit$beta) || is.null(fit$eta) || is.null(fit$t0)){ +# stop("Beta, Eta and/or t0 are not available.") +# }else{ +# message("calculateSingleConf: Currently, confidence bounds for Weibull 3P are not supported.") +# } +#} +#if(tolower(fit$options$dist) %in% c("lognormal","lognormal2p")){ +# message("calculateSingleConf: Currently, confidence bounds for lognormal are not supported.") +#} + }else{ + stop("Distribution type was not provided.") + } + }else{ + if(opadata$verbosity >= 1) + # TODO: using opadata since no other location of $verbosity is available. + message("calculateSingleConf: The fit argument is empty or contains no fits.") + } + } + fit +} diff --git a/R/calculateSingleFit.R b/R/calculateSingleFit.R new file mode 100644 index 0000000..1afe3cd --- /dev/null +++ b/R/calculateSingleFit.R @@ -0,0 +1,426 @@ +# R package 'abrem' +# Abernethy Reliability Methods +# Implementations of lifetime data analysis methods described in +# 'The New Weibull Handbook, Fifth edition' by Dr. Robert B. Abernethy. +# August 2014, Jurgen Symynck +# Copyright 2014, Jurgen Symynck +# +# For more info, visit http://www.openreliability.org/ +# +# For the latest version of this file, check the Subversion repository at +# http://r-forge.r-project.org/projects/abernethy/ +# +# Disclaimer: +# The author is not affiliated with Dr. Abernethy or Wes Fulton - CEO of +# Fulton Findings(TM) and author of the software package SuperSMITH +#------------------------------------------------------------------------------- +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . +# +# +-----------------------------------+ +# | execute this software with R: | +# | http://www.r-project.org/ | +# +-----------------------------------+ + +calculateSingleFit <- function(x,...){ + # x is a single Abrem object + # TODO: check the difference between + # x$fit[[i]]$options$dist <<- "weibull2p" + # ... and + # x$fit[[i]]$options$dist <- "weibull2p" + + + ######################### + # auxiliary functions # + ######################### + opadata <- x$options + opafit <- modifyList(opadata,list(...)) + vm <- function(vlevel,mess)if(opafit$verbosity >= vlevel)message(mess) + debug1 <- function()vm(2,paste0( + "calculateSingleFit: Attempting ",opafit$dist," (", + paste0(opafit$method.fit,collapse=", "), + '), ppos = \"',opafit$ppos,'\" fit...')) + debug2 <- function()vm(2,paste0( + "calculateSingleFit: Attempting ",opafit$dist," (", + paste0(opafit$method.fit,collapse=", "),") fit...")) + neededcolumns <- function(pppmethod=NULL){ + # This function selects the appropriate ppp columns for fitting + pppcolumn <- function(colname,pppmethod){ + na <- unlist(strsplit(tolower(colname),".",TRUE)) + identical(na[1],"ppp") && identical(na[2],tolower(pppmethod)) + } + basis <- x$data[,c("time","event")] + if(is.null(pppmethod)){ + basis + }else{ + wh <- which(sapply(names(x$data),pppcolumn,pppmethod,USE.NAMES=FALSE)) + cbind(basis,ppp=x$data[,wh]) + # TODO: explore usefulness of having the ranking methodd embedded in the + # column name at the fit level. + + } + } + goodness_of_fit <- function(npar){ + if(!is.null(x$fit[[i]])){ + if(is.null(x$fit[[i]]$gof)){ + vm(2,"calculateSingleFit: calculating r^2 using cor()...") + x$fit[[i]]$gof <<- list() + x$fit[[i]]$gof$r2 <<- cor(times, ranks, use = "complete.obs")^2 + }else vm(2,"calculateSingleFit: r^2 was already set...") + # if !NULL, then it was already set by the cpp version of the fitting method + + if(npar==2){ + if(identical(x$fit[[i]]$gof$r2,1)){ + vm(2,"calculateSingleFit: r^2 is exactly 1, bypassing AbPval and CCC^2 calculations...") + if(is.null(x$fit[[i]]$gof$AbPval))x$fit[[i]]$gof$AbPval <<- 100 + # TODO: infinity or 100% ? + }else{ + vm(2,"calculateSingleFit: r^2 is lower than 1...") + vm(2,"calculateSingleFit: Calculating AbPval and CCC^2...") + # TODO: prr: see pivotalMC + if(!any(c("mle","mle-rba","mle2","mle2-rba","mle3","mle3-rba") + # TODO: does it make sense to have reports of r^2 and ccc^2 while using MLE? + # Jurgen: yes, for comparing with the rank regression values + # what does Jabob think? + %in% tolower(opafit$method.fit))){ + gof <- abremPivotals::AbPval(x$fit[[i]]$fail,x$fit[[i]]$gof$r2) + if(is.null(x$fit[[i]]$gof$AbPval)) + x$fit[[i]]$gof$AbPval <<- gof[['AbPval']] + if(is.null(x$fit[[i]]$gof$CCC2)) + x$fit[[i]]$gof$CCC2 <<- gof[[2]] + } + } + }else{ + vm(2,paste("calculateSingleFit: bypassing AbPval and CCC^2", + "calculation for three parameter models...")) + } + }else{ + vm(1,"calculateSingleFit: no fit available for goodness-of-fit calculation.") + } + } + + ######################## + # main function body # + ######################## + i <- 1 + atLeastOneFit <- FALSE + if(is.null(x$fit)){ + vm(2,"calculateSingleFit: Creating the first fit in the abrem object...") + i <- 1 + x$fit <- list() + }else{ + vm(2,"calculateSingleFit: Appending a new fit to the existing abrem object...") + i <- length(x$fit)+1 + } + x$fit[[i]] <- list() + op <- unique(c(names(x$options),names(opafit))) + # this is needed to add options from opafit into li that + # are NULL in x$options + # TODO:tolower() needed? + if(length(li <- opafit[sapply(op,function(y){ + !identical(x$options[[y]], opafit[[y]])})]) > 0){ + x$fit[[i]]$options <- li + # the above enlists only options that are different from the abrems + # 'main' options. This excludes options$dist and options$method.fit + } + x$fit[[i]]$n <- x$n + x$fit[[i]]$fail <- x$fail + x$fit[[i]]$susp <- x$susp + + if(!is.null(opafit$importfile)){ + if(opafit$verbosity >= 1)message("calculateSingleFit : ", + "importing fit results from superSMITH report file\n",opafit$importfile) + try(fi <- file(opafit$importfile)) + if(!is.null(fi)){ + try(re <- readLines(fi)) + if(!is.null(re)){ + extract <- function(string) + na.omit(as.numeric(unlist(strsplit(gsub(",",".",string),"[^0123456789.]")))) + he <- data.frame(do.call("rbind",lapply(re[1:10],extract))) + atLeastOneFit <- TRUE + + x$fit[[i]]$options$dist <- "weibull2p" + x$fit[[i]]$options$method.fit <- "superSMITH" + x$fit[[i]]$eta <- he[3,3] + x$fit[[i]]$beta <- he[3,4] + x$fit[[i]]$gof <- list() + x$fit[[i]]$gof$r2 <- he[2,3] + x$fit[[i]]$gof$prr <- he[2,1] + x$fit[[i]]$gof$CCC2 <- he[2,5] + + x$fit[[i]]$n <- he[5,1] + x$fit[[i]]$fail <- he[5,1]-he[5,2] + x$fit[[i]]$susp <- he[5,2] + } + close(fi) + } + return(x) + } + if(any(c("rr","rr2") %in% tolower(opafit$method.fit))){ + # ____ _ _ + # | _ \ __ _ _ __ | | __ _ __ ___ __ _ _ __ ___ ___ ___(_) ___ _ __ + # | |_) / _` | '_ \| |/ / | '__/ _ \/ _` | '__/ _ \/ __/ __| |/ _ \| '_ \ + # | _ < (_| | | | | < | | | __/ (_| | | | __/\__ \__ \ | (_) | | | | + # |_| \_\__,_|_| |_|_|\_\ |_| \___|\__, |_| \___||___/___/_|\___/|_| |_| + # |___/ + debug1() + x$fit[[i]]$data <- neededcolumns(opafit$ppos[1]) + times <- log(x$fit[[i]]$data$time) + npar <- 2 + to_rr2 <- FALSE + # used for redirecting from rr to rr2 in case of three parameter distributions + xy <- ifelse(is_xony <- ("xony" %in% tolower(opafit$method.fit)),"xony","yonx") + if(tolower(opafit$dist) %in% c("weibull","weibull2p","weibull3p")){ + ranks <- log(qweibull(x$fit[[i]]$data[['ppp',exact=FALSE]],1,1)) + x$fit[[i]]$options$dist <- "weibull2p" + dst <- "weibull" + if(tolower(opafit$dist) %in% c("weibull3p")){ + x$fit[[i]]$options$dist <- "weibull3p" + npar <- 3 + } + } + if(tolower(opafit$dist) %in% c("lognormal","lognormal2p","lognormal3p")){ + ranks <- log(qlnorm(x$fit[[i]]$data[['ppp',exact=FALSE]],0,1)) + x$fit[[i]]$options$dist <- "lognormal2p" + dst <- "lognormal" + if(tolower(opafit$dist) %in% c("lognormal3p")){ + x$fit[[i]]$options$dist <- "lognormal3p" + npar <- 3 + } + } + + if("rr" %in% tolower(opafit$method.fit)){ + x$fit[[i]]$options$method.fit <- c("rr",xy) + if(npar==2){ + atLeastOneFit <- TRUE + if(is_xony) x$fit[[i]]$lm <- lm(times ~ ranks,x$fit[[i]]$data) + else x$fit[[i]]$lm <- lm(ranks ~ times,x$fit[[i]]$data) + # TODO: add error checking + B <- coef(x$fit[[i]]$lm)[[2]] + A <- coef(x$fit[[i]]$lm)[[1]] + if(dst=="lognormal"){ + x$fit[[i]]$mulog <- ifelse(is_xony,A,-A/B) + x$fit[[i]]$sigmalog <- ifelse(is_xony,B,1/B) + } + if(dst=="weibull"){ + x$fit[[i]]$eta <- ifelse(is_xony,exp(A),exp(-A/B)) + x$fit[[i]]$beta <- ifelse(is_xony,1/B,B) + } + }else{ + vm(0,paste0('calculateSingleFit: \"rr\" is not defined for three parameter distributions,\n', + 'defaulting to rr2...')) + to_rr2 <- TRUE + } + } + + if("rr2" %in% tolower(opafit$method.fit) || to_rr2){ + x$fit[[i]]$options$method.fit <- c("rr2",xy) + ret <- NULL + me <- NULL + try(ret <- abremPivotals::lslr( + x$fit[[i]]$data, + dist=dst, + npar=npar, + reg_method=xy)) + if(!is.null(ret) && !any(is.nan(ret))){ + atLeastOneFit <- TRUE + if(tolower(opafit$dist) %in% c("weibull","weibull2p","weibull3p")){ + x$fit[[i]]$eta <- ret[['Eta']] + x$fit[[i]]$beta <- ret[['Beta']] + if(tolower(opafit$dist) %in% c("weibull3p")){ + x$fit[[i]]$t0 <- ret[['t0']] + } + } + if(tolower(opafit$dist) %in% c("lognormal","lognormal2p","lognormal3p")){ + x$fit[[i]]$mulog <- ret[['Mulog']] + x$fit[[i]]$sigmalog <- ret[['Sigmalog']] + if(tolower(opafit$dist) %in% c("lognormal3p")){ + x$fit[[i]]$t0 <- ret[['t0']] + } + } + x$fit[[i]]$gof <- list() + x$fit[[i]]$gof$r2 <- ret[['Rsqr']] + if(npar==2) x$fit[[i]]$gof$AbPval <- ret[['AbPval']] + # AbPval is not supported for three parameter distributions + if(!is.null(opafit$threshold)){ + # this overwrites any threshold setting at the data level with a number + # TODO: this is not the way to go when trying to implement support for + # threshold with plot.abrem() + if(is.logical(opafit$threshold) && opafit$threshold) + x$options$threshold <- ret[[3]] + } + }else{ + vm(1,"calculateSingleFit: Fitting failed.") + x$fit[i] <- list(NULL) + # note that is.null(x$fit[[i]]) will exit with an error + # TODO: replace with x$fit[[i]] <- list(NULL) + } + } + goodness_of_fit(npar) + } + + + if(any(c("mle","mle-rba","mle2","mle2-rba","mle3","mle3-rba") %in% tolower(opafit$method.fit))){ + # from email David Silkworth ma 4/11/2013 18:29 + # The MLE you want to call in abrem is MLEw2p_cpp This is the fast one in compiled code. + # + # The other two are simply demonstrations in R. + # - MLEw2p_abrem solves the MLE using the method shown in The Handbook. + # The likelihood function has been differentiated analytically + # so as to separate beta from eta (by others as can be found in + # the literature). Then a relatively simple root finder algorithm + # (Newton method) identifies the beta_hat, followed by + # calculation of eta_hat from the separated function. + # This method is detailed in Appendix C, section C.4 of the Handbook. + # + # - MLEw2p_optim solves the MLE using the optim function the same + # way that surv package would. The optim function is being + # called with a default for the Nelder-Meade "simplex" algorithm. + # + # For the Cpp implementation it was simpler for me to port the + # Newton method as I approached the weibull first (before lognormal). + # Later, in order to do the Cpp implementation of the lognormal MLE, + # I indeed ended up coding a Nelder-Meade simplex algorithm from + # scratch. That algorithm became the basis of the port that is + # called by MLEln2p_cpp. I left the development files in the + # package for eventual tutorial value. + + # "mle" = MLEw2p_abrem + # "mle2" = MLEw2p_cpp , should be default + # "mle3" = MLEw2p_optim + # __ __ _ _____ + # | \/ | | | ____| + # | |\/| | | | _| + # | | | | |___| |___ + # |_| |_|_____|_____| + + debug2() + x$fit[[i]]$data <- neededcolumns(opafit$ppos[1]) + + fa <- x$fit[[i]]$data$time[x$fit[[i]]$data$event==1] + su <- x$fit[[i]]$data$time[x$fit[[i]]$data$event==0] + ret <- NULL + #is_3p <- FALSE + npar <- 2 + if(tolower(opafit$dist) %in% c("weibull","weibull2p","weibull3p")){ + if(tolower(opafit$dist) %in% c("weibull","weibull2p")){ + x$fit[[i]]$options$dist <- "weibull2p" + if(any(c("mle","mle-rba") %in% tolower(opafit$method.fit))){ + x$fit[[i]]$options$method.fit <- "mle" + try(ret <- debias::MLEw2p_abrem(fa,s=su)) + } + if(any(c("mle2","mle2-rba") %in% tolower(opafit$method.fit))){ + x$fit[[i]]$options$method.fit <- "mle2" + try(ret <- debias::MLEw2p_cpp(fa,s=su)) + # TODO: bypass the R code and use the CPP code immediately + } + if(any(c("mle3","mle3-rba") %in% tolower(opafit$method.fit))){ + x$fit[[i]]$options$method.fit <- "mle3" + try(ret <- debias::MLEw2p_optim(fa,s=su)) + # TODO: bypass the R code and use the CPP code immediately + } + } + if(tolower(opafit$dist) %in% c("weibull3p")){ + #is_3p <- TRUE + npar <- 3 + x$fit[[i]]$options$dist <- "weibull3p" + if(any(c("mle2","mle2","mle2-rba","mle3-rba") %in% tolower(opafit$method.fit))){ + vm(0,paste0( + 'calculateSingleFit: \"mle2\" and \"mle3\" are not defined for ", + opafit$dist,", defaulting to \"mle\"...')) + x$fit[[i]]$options$method.fit <- "mle"} + try(ret <- debias::MLEw3p_secant(fa,s=su)) + # secant: purely R code ... + } + if(!is.null(ret)){ + atLeastOneFit <- TRUE + x$fit[[i]]$beta <- ret[[2]] + x$fit[[i]]$eta <- ret[[1]] + if(npar==3){ + x$fit[[i]]$t0 <- ret[[3]] + x$fit[[i]]$gof <- list() + x$fit[[i]]$gof$loglik <- ret[[4]] + }else{ + x$fit[[i]]$gof <- list() + x$fit[[i]]$gof$loglik <- ret[[3]] + } + if(!is.null(opafit$threshold)){ + if(is.logical(opafit$threshold) && opafit$threshold) + x$options$threshold <- ret[[3]] + } + # this overwrites any threshold setting at the data level with a number + # TODO: this is not the way to go when trying to implement support for + # threshold with plot.abrem() + if(any(c("mle-rba","mle2-rba","mle3-rba") %in% tolower(opafit$method.fit))){ + if("mle-rba" %in% tolower(opafit$method.fit)) + x$fit[[i]]$options$method.fit <- "mle-rba" + if("mle2-rba" %in% tolower(opafit$method.fit)) + x$fit[[i]]$options$method.fit <- "mle2-rba" + if("mle3-rba" %in% tolower(opafit$method.fit)) + x$fit[[i]]$options$method.fit <- "mle3-rba" + vm(2,"calculateSingleFit: Applying Abernethy's Bias Reduction ...") + x$fit[[i]]$beta <- ret[[2]]*debias::RBAbeta(length(fa)) + # TODO: set the option: median or mean bias reduction + } + goodness_of_fit(npar) + }else{ + vm(1,"calculateSingleFit: Fitting failed.") + x$fit[i] <- list(NULL) + } + } + + if(tolower(opafit$dist) %in% c("lognormal","lognormal2p")){ + x$fit[[i]]$options$dist <- "lognormal2p" + if(any(c("mle","mle-rba","mle3","mle3-rba") %in% tolower(opafit$method.fit))){ + vm(0,paste0( + 'calculateSingleFit: \"mle\" and \"mle3\" are not defined for ', + opafit$dist,', defaulting to \"mle2\"...'))} + x$fit[[i]]$options$method.fit <- "mle2" + try(ret <- debias::MLEln2p_cpp(fa,s=su)) + if(!is.null(ret)){ + atLeastOneFit <- TRUE + x$fit[[i]]$mulog <- ret[[1]] + x$fit[[i]]$sigmalog <- ret[[2]] + x$fit[[i]]$gof <- list() + x$fit[[i]]$gof$loglik <- ret[[3]] + if(any(c("mle-rba","mle2-rba","mle3-rba") %in% tolower(opafit$method.fit))){ + if("mle-rba" %in% tolower(opafit$method.fit)) + x$fit[[i]]$options$method.fit <- "mle-rba" + if("mle2-rba" %in% tolower(opafit$method.fit)) + x$fit[[i]]$options$method.fit <- "mle2-rba" + if("mle3-rba" %in% tolower(opafit$method.fit)) + x$fit[[i]]$options$method.fit <- "mle3-rba" + vm(2,"calculateSingleFit: Applying Abernethy's Median Bias Reduction ...") + x$fit[[i]]$sigmalog <- ret[[2]]*debias::RBAsigma(length(fa)) + # with RBAsigma, there are no options... + } + goodness_of_fit(npar) + }else{ + vm(1,"calculateSingleFit: Fitting failed.") + x$fit[i] <- list(NULL) + + } + } + } + if(!atLeastOneFit){ + message("*** calculateSingleFit: Nothing has been fitted. ***\n", + '*** Does \"method.fit\" include sensible options? ***') + # x$fit[[i]] <- NULL + } + if(is.numeric(opafit$threshold)) + x$options$threshold <- opafit$threshold + # overwrite any previously set data-level-t0 to the one specified as an argument to abrem.fit() + # Don't know why - here - you MUST use <- in favor of <<- ... + x + # return a single Abrem object +} \ No newline at end of file diff --git a/R/contour.abrem.R b/R/contour.abrem.R new file mode 100644 index 0000000..279c5d1 --- /dev/null +++ b/R/contour.abrem.R @@ -0,0 +1,125 @@ +# R package 'abrem' +# Abernethy Reliability Methods +# Implementations of lifetime data analysis methods described in +# 'The New Weibull Handbook, Fifth edition' by Dr. Robert B. Abernethy. +# August 2014, Jurgen Symynck +# Copyright 2014, Jurgen Symynck +# +# For more info, visit http://www.openreliability.org/ +# +# For the latest version of this file, check the Subversion repository at +# http://r-forge.r-project.org/projects/abernethy/ +# +# Disclaimer: +# The author is not affiliated with Dr. Abernethy or Wes Fulton - CEO of +# Fulton Findings(TM) and author of the software package SuperSMITH +#------------------------------------------------------------------------------- +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . +# +# +-----------------------------------+ +# | execute this software with R: | +# | http://www.r-project.org/ | +# +-----------------------------------+ + +contour.abrem <- function(x,...){ + # +------------------------------+ + # | move abrem objects to list | + # | of abrem objects | + # +------------------------------+ + if(identical(class(x),"abrem")) x <- list(x) + if(!all(sapply(x,function(x)identical(class(x),"abrem")))){ + stop("Argument \"x\" is not of class \"abrem\" or ", + "a list of \"abrem\" objects.") + } + # as of this point, x is always a list of one or more abrem objects + + # +------------------------------------+ + # | create default options arguments | + # +------------------------------------+ + opa <- x[[1]]$options + opa <- modifyList(opa, list(...)) + + + # +--------------------------+ + # | create new plot canvas | + # +--------------------------+ + contourRanges <- findContourRanges(x,opa$verbosity) + if(!is.null(contourRanges)){ + xlimits <- range(contourRanges[,1]) + ylimits <- range(contourRanges[,2]) + opanames <- names(opa) + plotargs <- c(list(x=NA,axes=TRUE), + opa[opanames %in% plot_default_args()]) + plotargs$xlim <- xlimits + plotargs$ylim <- ylimits + plotargs$main <- opa$main.contour + plotargs$sub <- opa$sub.contour + plotargs$log <- "" + plotargs$xlab <- "Eta" + plotargs$ylab <- "Beta" + + do.call("plot.default",plotargs) + if(opa$is.plot.grid){ + abline( + h=pretty(contourRanges[,2],10), + v=pretty(contourRanges[,1],10), + col = opa$col.grid) + # TODO: add userchoice in grid density here + } + }else message("contour.abrem: No contours available in (list of) abrem objects.") +# r <- seq.log(opa$xlim[1]/10,opa$xlim[2]*10,c(1,5)) + + # +------------------+ + # | plot contours | + # +------------------+ + plotContours <- function(abrem){ + if(!is.null(abrem$fit)){ + plotContours2 <- function(fit){ + if(!is.null(fit$options)){ + opafit <- modifyList(abrem$options,fit$options) + }else{opafit <- abrem$options} + is_MLE <- any(c("mle","mle-rba") %in% tolower(fit$options$method.fit)) + if(!is.null(fit$conf$blives)){ + plotContours3 <- function(blicon){ + if(!is.null(blicon$options)){ + opaconf <- modifyList(opafit,blicon$options) + }else{opaconf <- opafit} + opaconf <- modifyList(opaconf,list(...)) + if(!is.null(blicon$MLEXContour)){ + if(!is_MLE) + points(blicon$MLEXContour$MLEpoint,pch=20,col=abrem$options$col) + # if MLE or MLE-RBA was used to calculate the distribution + # parameters, the omit plotting the MLEpoint. In all other cases, + # plot the MLEpoint because the distribution parameters will not match exactly + # the MLE point + if(all(c(!is.null(fit$beta),!is.null(fit$eta)))) + points(x=fit$eta,y=fit$beta,pch=abrem$options$pch,col=abrem$options$col, + lwd=abrem$options$lwd.points,cex=abrem$options$cex.points) + points(blicon$MLEXContour[[1]],type="l",lwd=opaconf$lwd,lty=opaconf$lty,col=opaconf$col) + } + } + #mtrace(plotContours3) + do.call("rbind",lapply(fit$conf$blives,plotContours3)) + # combine the ranges from all MLEXContours + # found in the list of blicons + } + } + do.call("rbind",lapply(abrem$fit,plotContours2)) + # combine the ranges from all MLEXContours + # found in the list of fits + } + } + if(!is.null(contourRanges)) lapply(x,plotContours) + invisible() +} diff --git a/R/findContourRanges.R b/R/findContourRanges.R new file mode 100644 index 0000000..ef6745c --- /dev/null +++ b/R/findContourRanges.R @@ -0,0 +1,72 @@ +# R package 'abrem' +# Abernethy Reliability Methods +# Implementations of lifetime data analysis methods described in +# 'The New Weibull Handbook, Fifth edition' by Dr. Robert B. Abernethy. +# August 2014, Jurgen Symynck +# Copyright 2014, Jurgen Symynck +# +# For more info, visit http://www.openreliability.org/ +# +# For the latest version of this file, check the Subversion repository at +# http://r-forge.r-project.org/projects/abernethy/ +# +# Disclaimer: +# The author is not affiliated with Dr. Abernethy or Wes Fulton - CEO of +# Fulton Findings(TM) and author of the software package SuperSMITH +#------------------------------------------------------------------------------- +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . +# +# +-----------------------------------+ +# | execute this software with R: | +# | http://www.r-project.org/ | +# +-----------------------------------+ + +contourRange <- function(MLEXContour){ + ra <- do.call("rbind",MLEXContour) + data.frame(range(ra[,1]),range(ra[,2])) +} + +findContourRanges <- function(x,v){ + # +---------------------------------------------------+ + # | find absolute maximum and minimum contour range | + # | over the (list of) abrem objects | + # +---------------------------------------------------+ + # x is always a list of abrem object(s) + + findrange1 <- function(abrem){ + if(!is.null(abrem$fit)){ + findrange2 <- function(fit){ + if(!is.null(fit$conf$blives)){ + findrange3 <- function(blicon){ + if(!is.null(blicon$MLEXContour)){ + # a contour is available + #if(deb)mtrace(contourRange) + contourRange(blicon$MLEXContour) + } + } + #if(deb)mtrace(findrange3) + do.call("rbind",lapply(fit$conf$blives,findrange3)) + # combine the ranges from all MLEXContours + # found in the list of blicons + } + } + #if(deb)mtrace(findrange2) + do.call("rbind",lapply(abrem$fit,findrange2)) + # combine the ranges from all MLEXContours + # found in the list of fits + } + } + #if(deb)mtrace(findrange1) + do.call("rbind",lapply(x,findrange1)) +} diff --git a/R/findMaxDataRange.R b/R/findMaxDataRange.R new file mode 100644 index 0000000..7b651e5 --- /dev/null +++ b/R/findMaxDataRange.R @@ -0,0 +1,53 @@ +# R package 'abrem' +# Abernethy Reliability Methods +# Implementations of lifetime data analysis methods described in +# 'The New Weibull Handbook, Fifth edition' by Dr. Robert B. Abernethy. +# August 2014, Jurgen Symynck +# Copyright 2014, Jurgen Symynck +# +# For more info, visit http://www.openreliability.org/ +# +# For the latest version of this file, check the Subversion repository at +# http://r-forge.r-project.org/projects/abernethy/ +# +# Disclaimer: +# The author is not affiliated with Dr. Abernethy or Wes Fulton - CEO of +# Fulton Findings(TM) and author of the software package SuperSMITH +#------------------------------------------------------------------------------- +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . +# +# +-----------------------------------+ +# | execute this software with R: | +# | http://www.r-project.org/ | +# +-----------------------------------+ + +findMaxDataRange <- function(x,v,log=""){ + # +-------------------------------------------+ + # | find absolute maximum and minimum range | + # | over the (list of) abrem objects | + # +-------------------------------------------+ + # x is always a list of abrem object(s) + if(all(sapply(x,function(x)identical(class(x),"abrem")))){ + ret <- do.call("rbind",lapply(x,findRangeInAbrem,v,log)) + }else{ + stop("Argument \"x\" is no list of \"abrem\" objects.") + } + # TODO: the above still needed? because x is always list of abrems? + if(tolower(log) %in% c("x","xy","yx")){ + # if log scale is to be used then omit zero and negative time values + # from the range dataset + ret[ret$xrange <= 0 , 1] <- NA + } + ret +} diff --git a/R/findRangeInAbrem.R b/R/findRangeInAbrem.R new file mode 100644 index 0000000..2955089 --- /dev/null +++ b/R/findRangeInAbrem.R @@ -0,0 +1,51 @@ +# R package 'abrem' +# Abernethy Reliability Methods +# Implementations of lifetime data analysis methods described in +# 'The New Weibull Handbook, Fifth edition' by Dr. Robert B. Abernethy. +# August 2014, Jurgen Symynck +# Copyright 2014, Jurgen Symynck +# +# For more info, visit http://www.openreliability.org/ +# +# For the latest version of this file, check the Subversion repository at +# http://r-forge.r-project.org/projects/abernethy/ +# +# Disclaimer: +# The author is not affiliated with Dr. Abernethy or Wes Fulton - CEO of +# Fulton Findings(TM) and author of the software package SuperSMITH +#------------------------------------------------------------------------------- +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . +# +# +-----------------------------------+ +# | execute this software with R: | +# | http://www.r-project.org/ | +# +-----------------------------------+ + +findRangeInAbrem <- function(abrem,v,log=""){ + if(!is.null(abrem$data)){ + if(!is.null(abrem$data$time)){ + ret <- data.frame(xrange=range(abrem$data$time,na.rm=TRUE)) + }else{ + stop("$data contains no \"$time\" column -> ", + "cannot create plot canvas.") + } + if(length(pppcolumns <- grep("ppp",names(abrem$data))) >= 1){ + ret <- cbind(ret,yrange=range(abrem$data[,pppcolumns],na.rm=TRUE)) + }else{ + stop("$data contains no ppp column(s) -> ", + "cannot create plot canvas.") + } + }else{stop('Argument \"x\" contains no \"$data\" dataframe.')} + ret + } \ No newline at end of file diff --git a/R/getPlotRange.R b/R/getPlotRange.R new file mode 100644 index 0000000..f468ed8 --- /dev/null +++ b/R/getPlotRange.R @@ -0,0 +1,45 @@ +# R package 'abrem' +# Abernethy Reliability Methods +# Implementations of lifetime data analysis methods described in +# 'The New Weibull Handbook, Fifth edition' by Dr. Robert B. Abernethy. +# August 2014, Jurgen Symynck +# Copyright 2014, Jurgen Symynck +# +# For more info, visit http://www.openreliability.org/ +# +# For the latest version of this file, check the Subversion repository at +# http://r-forge.r-project.org/projects/abernethy/ +# +# Disclaimer: +# The author is not affiliated with Dr. Abernethy or Wes Fulton - CEO of +# Fulton Findings(TM) and author of the software package SuperSMITH +#------------------------------------------------------------------------------- +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . +# +# +-----------------------------------+ +# | execute this software with R: | +# | http://www.r-project.org/ | +# +-----------------------------------+ + +getPlotRangeX <- function(log){ + if(log %in% c("x","xy","yx")) 10^par("usr")[1:2] + else par("usr")[1:2] +} + +getPlotRangeY <- function(log){ + if(log %in% c("y","xy","yx")) 10^par("usr")[3:4] + else par("usr")[3:4] +# if(log %in% c("y","xy","yx"))) 10^par("usr")[1:2] +# else par("usr")[1:2] +} diff --git a/R/legendConf.R b/R/legendConf.R new file mode 100644 index 0000000..80a160d --- /dev/null +++ b/R/legendConf.R @@ -0,0 +1,76 @@ +# R package 'abrem' +# Abernethy Reliability Methods +# Implementations of lifetime data analysis methods described in +# 'The New Weibull Handbook, Fifth edition' by Dr. Robert B. Abernethy. +# August 2014, Jurgen Symynck +# Copyright 2014, Jurgen Symynck +# +# For more info, visit http://www.openreliability.org/ +# +# For the latest version of this file, check the Subversion repository at +# http://r-forge.r-project.org/projects/abernethy/ +# +# Disclaimer: +# The author is not affiliated with Dr. Abernethy or Wes Fulton - CEO of +# Fulton Findings(TM) and author of the software package SuperSMITH +#------------------------------------------------------------------------------- +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . +# +# +-----------------------------------+ +# | execute this software with R: | +# | http://www.r-project.org/ | +# +-----------------------------------+ +# +legendConf <- function(fit,conftype,opadata,...){ + if(!is.null(fit)){ + # TODO: what happens when trying to plot superSMITH confidence bounds + # when no fit data are available? + if(!is.null(fit$options)) + opafit <- modifyList(opadata,fit$options) + opafit <- modifyList(opafit,list(...)) + if(identical(tolower(conftype),"blives")){ + if(!is.null(fit$conf$blives)){ + for.each.blicon <- function(blicon){ + if(!is.null(blicon$options)){ + opaconf <- modifyList(opafit,blicon$options) + }else{opaconf <- opafit} + if(opaconf$in.legend){ + # TODO: correct usage of this logical value? + li <- list() + li[[1]] <- bsll(legend=paste0("B-lives, type = ", + ifelse(is.null(blicon$type),"NA", + paste0("\"",blicon$type,"\""))), + col=opaconf$col,lwd=opaconf$lwd,lty=opaconf$lty) + li[[2]] <- bsll(legend=paste0(" CL = ", + ifelse(is.null(blicon$cl),"NA", + paste0(signif(blicon$cl*100,4)," [%]")), + ifelse(is.null(blicon$S),"", + paste0(", S = ",blicon$S)))) + if(opaconf$in.legend.blives){ + params <- unlist(list(beta=fit$beta,eta=fit$eta,t0=fit$t0, + mulog=fit$mulog,sigmalog=fit$sigmalog,rate=fit$rate)) + if(is.null(bl <- blicon$unrel))bl <- opaconf$unrel + fu <- function(bl){ + bsll(legend=Blifestring(bl,blicon,opafit$signif,params)) + } + c(li,lapply(bl,fu)) + }else(li) + }else NULL + } + unlist(lapply(fit$conf$blives,for.each.blicon),FALSE) + # TODO: replace by do.call ? + }else{NULL} + } + }else{NULL} +} \ No newline at end of file diff --git a/R/options.abrem.R b/R/options.abrem.R new file mode 100644 index 0000000..afa6f0f --- /dev/null +++ b/R/options.abrem.R @@ -0,0 +1,140 @@ +# R package 'abrem' +# Abernethy Reliability Methods +# Implementations of lifetime data analysis methods described in +# 'The New Weibull Handbook, Fifth edition' by Dr. Robert B. Abernethy. +# August 2014, Jurgen Symynck +# Copyright 2014, Jurgen Symynck +# +# For more info, visit http://www.openreliability.org/ +# +# For the latest version of this file, check the Subversion repository at +# http://r-forge.r-project.org/projects/abernethy/ +# +# Disclaimer: +# The author is not affiliated with Dr. Abernethy or Wes Fulton - CEO of +# Fulton Findings(TM) and author of the software package SuperSMITH +#------------------------------------------------------------------------------- +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . +# +# +-----------------------------------+ +# | execute this software with R: | +# | http://www.r-project.org/ | +# +-----------------------------------+ + +options.abrem <- function(...){ + # function to handle the many options of the weibull toolkits functions + # the option list should only be manipulated through this function! + + # TODO: WARNING: partial matching is in effect! + # options.abrem()$ylim will return options.abrem()$ylim.default if + # $ylim was set to NULL! + + single <- FALSE + args <- list(...) + + if(!exists(as.character(substitute(options_abrem)))) + # if the globally accessible variable was not defined yet, then + # create it here with default values OR reset to default values + # message ("Resetting Weibulltoolkit options to default values...") + options_abrem <<- list( + dist="weibull", + method.fit=c("rr","xony"), + # TODO: decide between "rr" and "rr2". + # "rr2" is not 100%compatible with params.to.ob + conf.what="blives", + conf.blives.sides="double", + unrel.n=25, + method.conf.blives="mcpivotals", + ppos="benard",#,"beta"),#,"benard_cpp","hazen","mean"), + #"benard", "beta", "mean", "hazen", "Kaplan-Meier", "Blom". + # TODO: make the transition from pp to ppos everywhere in the abrem code! + # TODO: Code has changed from beta to benard as the default! + S=1e4, + pivotals=FALSE, + cl=0.9, + # cl or CL ? + unrel=c(0.1,0.05,0.01), + verbosity=0, + mar=c(5.1,4.1,5.1,2.1), + main="Probability Plot", + main.contour="Contour Plot", + # a default title for Contour plots + sub=NULL, + sub.contour=NULL, + xlim=NULL, + ylim=NULL, + xlab="Time To Failure", + ylab="Unreliability [%]", + log="x", # maybe this should be removed in favor of "canvas" + #canvas="weibull", + coordinate.text.size=0.7, + signif=4, + pch=1, + lwd=2, + lwd.points=2, + cex.points=1, + lty=1, + col="black", + col.grid="gray", + threshold=FALSE, + is.plot.grid=TRUE, + is.plot.fit=TRUE, + is.plot.ppp=TRUE, + is.plot.pppcoordinates=FALSE, + is.plot.legend=TRUE, + # legend.position="bottomright", + legend.text.size=0.7, + label="", + in.legend=TRUE, + in.legend.blives=TRUE, + in.legend.gof=TRUE, + is.plot.cb = TRUE, + persistent=TRUE) + + if (!length(args)) + args <- options_abrem + # return the current option list + else { +# if (all(unlist(lapply(args, is.character)))) + # if all items in the args are characters, then + # treat them as the names of the options. +# args <- as.list(unlist(args)) + # TODO; the above abrem causes bug 5596 + if (length(args) == 1) { + if (is.list(args[[1L]]) | is.null(args[[1L]])) + args <- args[[1L]] + # unlist the first (and only) argument to a string + else if(is.null(names(args))) + # if there is no name to args then + # the arg itself is the name (?) + single <- TRUE + } + } + value <- args + if(options_abrem$persistent){ + options_abrem <<-modifyList(options_abrem, value) + } + if(!is.null(args$persistent)){ + value <- args + if(args$persistent){ + options_abrem <<-modifyList(options_abrem, value) + } + } + # make the options stick between calls of options.abrem() + if(is.null(names(args))) + value <- options_abrem[match(args,names(options_abrem))] + if(single) value <- value[[1L]] + value +} +# TODO :options that are NULL are not shown in the printout diff --git a/R/params.to.ob.R b/R/params.to.ob.R new file mode 100644 index 0000000..caa63c8 --- /dev/null +++ b/R/params.to.ob.R @@ -0,0 +1,111 @@ +# R package 'abrem' +# Abernethy Reliability Methods +# Implementations of lifetime data analysis methods described in +# 'The New Weibull Handbook, Fifth edition' by Dr. Robert B. Abernethy. +# August 2014, Jurgen Symynck +# Copyright 2014, Jurgen Symynck +# +# For more info, visit http://www.openreliability.org/ +# +# For the latest version of this file, check the Subversion repository at +# http://r-forge.r-project.org/projects/abernethy/ +# +# Disclaimer: +# The author is not affiliated with Dr. Abernethy or Wes Fulton - CEO of +# Fulton Findings(TM) and author of the software package SuperSMITH +#------------------------------------------------------------------------------- +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . +# +# +-----------------------------------+ +# | execute this software with R: | +# | http://www.r-project.org/ | +# +-----------------------------------+ + +params.to.ob <- function(dist, ... ){ + # function to generate test data that result in a perfect fit + # (when using RR; X-on-Y, median plot positions) + # beta,eta: slope and shape parameters of 2 parameter Weibull + # event: either an integer with the number of complete observations, or + # a vector with censoring information (e.g.: c(1,1,1,0,0,0,1) + + opa <- options.abrem() + opa <- modifyList(opa, list(...)) + if(!missing(dist)){ + if(!is.null(opa$n)){ + if(opa$n >= 2){ + if(is.null(opa$event)){ + opa$event <- rep(1,opa$n) + }else{ + stop("Either Argument \"n\" or \"event\" ", + "should be supplied, not both.") + } + }else{ + stop("Argument \"n\" (number of failures) must be ", + "at least 2.") + } + } + ### assuming opa$event is present + if(!is.null(opa$event)){ + if(!(length(opa$event) >= 2)){ + stop("Argument \"n\" (number of failures) must be ", + "at least 2.") + } + }else{ + stop("Argument \"event\" should have length of at least 2.") + } + if(all(opa$event==0)){ + stop("Argument \"event\"contains only censored events.") + } + if(any(c("rr","rr2") %in% tolower(opa$method.fit))){ + # TODO: currently, only the "rr2" method is supported. + # Should result in the same numbers though... + # TODO: expand with the new abremPivotal argument "aranks" + if(tolower(dist) %in% c("weibull","weibull2p")){ + if(!is.null(opa$beta) && !is.null(opa$eta)){ + if(!is.null(opa$ppos)){ + ppp <- rep(NA,length(opa$event)) + ppp[opa$event==1] <- abremPivotals::getPPP(x=opa$event,ppos=opa$ppos[1])$ppp + # the above call SHOULD be correct (compated with output from pivotals + # 'mis'using which for getting timelike values from the event vector + # no need to use the real lifetime observations here, + # just getting the ppp + ret <- data.frame(time=qweibull(ppp,opa$beta,opa$eta),event=opa$event) + # a good thing that qweibull deals nicely with NA's! + }else{ + stop("Argument \"ppos\" is missing; no ranking method supplied.") + } + }else{stop("Arguments \"beta\" and/or \"eta\" not supplied.")} + } + if(tolower(dist) %in% c("lognormal","lognormal2p")){ + if(!is.null(opa$mulog) && !is.null(opa$sigmalog)){ + if(!is.null(opa$ppos)){ + ppp <- rep(NA,length(opa$event)) + ppp[opa$event==1] <- abremPivotals::getPPP(x=which(opa$event==1),#s=which(opa$event==0), + ppos=opa$ppos[1])$ppp + ret <- data.frame(time=qlnorm(ppp,opa$mulog,opa$sigmalog),event=opa$event) + }else{ + stop("Argument \"ppos\" is missing; no ranking method supplied.") + } + }else{stop("Arguments \"mulog\" and/or \"sigmalog\" not supplied.")} + } + }else{ + stop("Currently, only rank regression is supported.") + ret <- NULL + } + ret + }else{ + stop("Argument \"dist\" is missing.") + ret <- NULL + } +} diff --git a/R/plot.abrem.R b/R/plot.abrem.R new file mode 100644 index 0000000..b2536a3 --- /dev/null +++ b/R/plot.abrem.R @@ -0,0 +1,163 @@ +# R package 'abrem' +# Abernethy Reliability Methods +# Implementations of lifetime data analysis methods described in +# 'The New Weibull Handbook, Fifth edition' by Dr. Robert B. Abernethy. +# August 2014, Jurgen Symynck +# Copyright 2014, Jurgen Symynck +# +# For more info, visit http://www.openreliability.org/ +# +# For the latest version of this file, check the Subversion repository at +# http://r-forge.r-project.org/projects/abernethy/ +# +# Disclaimer: +# The author is not affiliated with Dr. Abernethy or Wes Fulton - CEO of +# Fulton Findings(TM) and author of the software package SuperSMITH +#------------------------------------------------------------------------------- +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . +# +# +-----------------------------------+ +# | execute this software with R: | +# | http://www.r-project.org/ | +# +-----------------------------------+ + +plot.abrem <- function(x,...){ + # +------------------------------+ + # | move abrem objects to list | + # | of abrem objects | + # +------------------------------+ + if(identical(class(x),"abrem")) x <- list(x) + if(!all(sapply(x,function(x)identical(class(x),"abrem")))){ + stop("Argument \"x\" is not of class \"abrem\" or ", + "a list of \"abrem\" objects.") + } + # as of this point, x is always a list of one or more abrem objects + + # +------------------------------------+ + # | create default options arguments | + # +------------------------------------+ + opa <- x[[1]]$options + opa <- modifyList(opa, list(...)) + + # +--------------------------------------+ + # | dealing with threshold parameters | + # +--------------------------------------+ + + if(!is.null(list(...)$threshold)) + message("Currently, passing the \'threshold\' argument to plot.abrem is not supported. Proceeding...") + + # +--------------------------+ + # | create new plot canvas | + # +--------------------------+ + ra <- findMaxDataRange(x,opa$verbosity,opa$log) + # NA values can be part of ra, when log scales are to be used + # and there are negative failure times + xlimits <- range(ra$xrange,na.rm=TRUE) + ylimits <- range(ra$yrange,na.rm=TRUE) + if(is.null(opa$xlim)){ + opa$xlim <- c(10^(log10(xlimits[1])-0.5), + 10^(log10(xlimits[2])+1)) + } + if(is.null(opa$ylim)){ + if(ylimits[1] < 0.01) opa$ylim <- c(signif(ylimits[1],1),0.99) + else opa$ylim <- c(0.01,0.99) + # do not care about the upper limit + } + opanames <- names(opa) + plotargs <- c(list(x=NA,axes=FALSE), + opa[opanames %in% plot_default_args()]) + if(!is.null(plotargs$ylim)){ + plotargs$ylim <- F0inv(plotargs$ylim,opa$log) + } + plotargs$main <- NULL + # do not plot "main" just yet... + if(!is.null(opa$mar))par(mar=opa$mar) + if(!is.null(opa$mai))par(mai=opa$mai) + do.call(plot.default,plotargs) + if(opa$is.plot.grid){ + abline( + h=F0inv(seq.wb(opa$ylim[1]/10,1-(1-opa$ylim[2])/10),opa$log), + v=seq.log(opa$xlim[1]/10,opa$xlim[2]*10,seq(0,10,1)), + col = opa$col.grid) + } + r <- seq.log(opa$xlim[1]/10,opa$xlim[2]*10,c(1,5)) + #lin <- 0.0 + for(t in c(1,3)){ + axis(t,at=seq.log(opa$xlim[1]/10,opa$xlim[2]*10,seq(0,10,0.2)), + labels=NA,tcl=-0.25)#,line=0.0 + # plot top and bottom axis tickmarks + axis(t,at=r,labels=r,tcl=-0.75)#,line=0.0 + # plot top and bottom axis labels + } + r <- c(seq.wb(opa$ylim[1]/10,1-(1-opa$ylim[2])/10,c(1,2,5)),0.9) + for(t in c(2,4)){ + # TODO: rewrite as do.call() or apply() + axis(t,at=F0inv(seq.wb(opa$ylim[1]/10,1-(1-opa$ylim[2])/10), + opa$log),labels=NA,tcl=-0.25)#,line=0.0 + # plot left and right axis tickmarks + axis(t,at=F0inv(r,opa$log), + labels=r*100,tcl=-0.75)#,line=0.0 + # plot left and right axis labels + } + abline(h=0,lty = 3,col = opa$col.grid) + # plot the 63.2 [%] unreliability line + title(main=opa$main,line=3) + + # +--------------------------+ + # | plot confidence bounds | + # +--------------------------+ + lapply(x,plotConfsInAbrem,v=opa$verbosity,...) + + # +-------------+ + # | plot fits | + # +-------------+ + lapply(x,plotFitsInAbrem,v=opa$verbosity,xl=opa$xlim,yl=opa$ylim,...) + + # +-----------------------------------+ + # | plot probability plot positions | + # +-----------------------------------+ + lapply(x,plotSingleDataSet,opa$is.plot.ppp,...) + + # +----------------+ + # | plot legends | + # +----------------+ + lolegends <- NULL + lolegends <- unlist(lapply(x,buildListOfLegends, + v=opa$verbosity, + isplotlegend=opa$is.plot.legend,...),FALSE) + # TODO: likely, unlist is NOT the best way to go here, investigate + lolegends <- lolegends[sapply(lolegends,function(lol)!is.null(lol))] + # omit any list entries that contain only NULL + if(!is.null(lolegends) && any(sapply(lolegends,function(lol)!is.null(lol)))){ + lx <- rep(lolegends[[1]]$rect$left,length(lolegends)) + ly <- lolegends[[1]]$rect$top + + c(0,cumsum(sapply(lolegends,function(le)le$rect$h)[-1])) + if(opa$log %in% c("x","xy","yx")) lx <- 10^lx + if(opa$log %in% c("y","xy","yx")) ly <- 10^ly + # TODO: F0(ly): looks very suspicious that this works -> investigate! + for(i in 1:length(lolegends)){ + plotSingleLegend(lolegends[[i]],lx[i],ly[i]) + # TODO: replace with lapply + } + }else{ + if(opa$verbosity >= 1)message( + "plot.abrem: There is no legend to plot.") + } +# if(opa$log == "x") legend("top",legend=NA,title="Weibull",bg="white") +# if(opa$log == "xy") legend("top",legend=NA,title="Lognormal",bg="white") +# if(opa$log %in% c("","y")) legend("top",legend=NA,title="xxx",bg="white") + invisible() + # TODO: return the abrem object with updated graphical options + # TODO: check if this makes sense when supplying a list +} diff --git a/R/plotConfsInAbrem.R b/R/plotConfsInAbrem.R new file mode 100644 index 0000000..1beaf2e --- /dev/null +++ b/R/plotConfsInAbrem.R @@ -0,0 +1,44 @@ +# R package 'abrem' +# Abernethy Reliability Methods +# Implementations of lifetime data analysis methods described in +# 'The New Weibull Handbook, Fifth edition' by Dr. Robert B. Abernethy. +# August 2014, Jurgen Symynck +# Copyright 2014, Jurgen Symynck +# +# For more info, visit http://www.openreliability.org/ +# +# For the latest version of this file, check the Subversion repository at +# http://r-forge.r-project.org/projects/abernethy/ +# +# Disclaimer: +# The author is not affiliated with Dr. Abernethy or Wes Fulton - CEO of +# Fulton Findings(TM) and author of the software package SuperSMITH +#------------------------------------------------------------------------------- +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . +# +# +-----------------------------------+ +# | execute this software with R: | +# | http://www.r-project.org/ | +# +-----------------------------------+ + +plotConfsInAbrem <- function(abrem,v,...){ + #opadata <- modifyList(x$options, arg) + if(!is.null(abrem$fit)){ + ret <- lapply(abrem$fit,plotConfsInFit,opadata=abrem$options,...) + }else{ + if(v >= 1)message( + "plotConfsInAbrem: This Abrem object contains no fits ", + "or confidence calculations.") + } +} \ No newline at end of file diff --git a/R/plotConfsInFit.R b/R/plotConfsInFit.R new file mode 100644 index 0000000..c27462b --- /dev/null +++ b/R/plotConfsInFit.R @@ -0,0 +1,45 @@ +# R package 'abrem' +# Abernethy Reliability Methods +# Implementations of lifetime data analysis methods described in +# 'The New Weibull Handbook, Fifth edition' by Dr. Robert B. Abernethy. +# August 2014, Jurgen Symynck +# Copyright 2014, Jurgen Symynck +# +# For more info, visit http://www.openreliability.org/ +# +# For the latest version of this file, check the Subversion repository at +# http://r-forge.r-project.org/projects/abernethy/ +# +# Disclaimer: +# The author is not affiliated with Dr. Abernethy or Wes Fulton - CEO of +# Fulton Findings(TM) and author of the software package SuperSMITH +#------------------------------------------------------------------------------- +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . +# +# +-----------------------------------+ +# | execute this software with R: | +# | http://www.r-project.org/ | +# +-----------------------------------+ + +plotConfsInFit <- function(fit,opadata,...){ + arg <- list(...) + if(!is.null(fit$conf$blives)){ + if(!is.null(fit$options)){ + opafit <- modifyList(opadata,fit$options) + }else{opafit <- opadata} + lapply(fit$conf$blives,plotSingleConfBound,opafit=opafit,opadatathreshold=opadata$threshold,...) + } +# else{if(arg$v >= 1)message(match.call()[[1]], +# ": This fit contains no confidence calculations for B-lives.")} +} \ No newline at end of file diff --git a/R/plotFitsInAbrem.R b/R/plotFitsInAbrem.R new file mode 100644 index 0000000..0c742c4 --- /dev/null +++ b/R/plotFitsInAbrem.R @@ -0,0 +1,44 @@ +# R package 'abrem' +# Abernethy Reliability Methods +# Implementations of lifetime data analysis methods described in +# 'The New Weibull Handbook, Fifth edition' by Dr. Robert B. Abernethy. +# August 2014, Jurgen Symynck +# Copyright 2014, Jurgen Symynck +# +# For more info, visit http://www.openreliability.org/ +# +# For the latest version of this file, check the Subversion repository at +# http://r-forge.r-project.org/projects/abernethy/ +# +# Disclaimer: +# The author is not affiliated with Dr. Abernethy or Wes Fulton - CEO of +# Fulton Findings(TM) and author of the software package SuperSMITH +#------------------------------------------------------------------------------- +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . +# +# +-----------------------------------+ +# | execute this software with R: | +# | http://www.r-project.org/ | +# +-----------------------------------+ + +plotFitsInAbrem <- function(abrem,v,xl,yl,...){ + opadata <- modifyList(abrem$options,list(xl,yl)) + if(!is.null(abrem$fit)){ + ret <- lapply(abrem$fit,plotSingleFit, + v=v,opadata=opadata,...) + }else{ + if(v >= 1)message( + "plotFitsInAbrem: This Abrem object contains no fits.") + } +} \ No newline at end of file diff --git a/R/plotSingleConfBound.R b/R/plotSingleConfBound.R new file mode 100644 index 0000000..054a36b --- /dev/null +++ b/R/plotSingleConfBound.R @@ -0,0 +1,60 @@ +# R package 'abrem' +# Abernethy Reliability Methods +# Implementations of lifetime data analysis methods described in +# 'The New Weibull Handbook, Fifth edition' by Dr. Robert B. Abernethy. +# August 2014, Jurgen Symynck +# Copyright 2014, Jurgen Symynck +# +# For more info, visit http://www.openreliability.org/ +# +# For the latest version of this file, check the Subversion repository at +# http://r-forge.r-project.org/projects/abernethy/ +# +# Disclaimer: +# The author is not affiliated with Dr. Abernethy or Wes Fulton - CEO of +# Fulton Findings(TM) and author of the software package SuperSMITH +#------------------------------------------------------------------------------- +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . +# +# +-----------------------------------+ +# | execute this software with R: | +# | http://www.r-project.org/ | +# +-----------------------------------+ + +plotSingleConfBound <- function(blc,opafit,opadatathreshold,...){ + if(!is.null(blc$options)){ + opaconf <- modifyList(opafit,blc$options) + }else{opaconf <- opafit} + opaconf <- modifyList(opaconf,list(...)) + if(opaconf$is.plot.cb){ + t0 <- 0 + if(is.logical(opadatathreshold))if(opafit$threshold) + warning ("opadata$threshold is a logical value but numeric value was expected. Proceeding...") + if(is.numeric(opadatathreshold))t0 <- opadatathreshold + # efffectively ignore any threshold argument set at the conf level + + if(!is.null(blc$bounds$datum)) + lines(y=F0inv(blc$bounds$unrel,opaconf$log), + x=blc$bounds$datum-t0, + col=opaconf$col,lwd=1,lty=2) + if(!is.null(blc$bounds$lower)) + lines(y=F0inv(blc$bounds$unrel,opaconf$log), + x=blc$bounds$lower-t0,col=opaconf$col, + lwd=opaconf$lwd,lty=opaconf$lty) + if(!is.null(blc$bounds$upper)) + lines(y=F0inv(blc$bounds$unrel,opaconf$log), + x=blc$bounds$upper-t0,col=opaconf$col, + lwd=opaconf$lwd,lty=opaconf$lty) + } +} diff --git a/R/plotSingleDataSet.R b/R/plotSingleDataSet.R new file mode 100644 index 0000000..07fd2b6 --- /dev/null +++ b/R/plotSingleDataSet.R @@ -0,0 +1,60 @@ +# R package 'abrem' +# Abernethy Reliability Methods +# Implementations of lifetime data analysis methods described in +# 'The New Weibull Handbook, Fifth edition' by Dr. Robert B. Abernethy. +# August 2014, Jurgen Symynck +# Copyright 2014, Jurgen Symynck +# +# For more info, visit http://www.openreliability.org/ +# +# For the latest version of this file, check the Subversion repository at +# http://r-forge.r-project.org/projects/abernethy/ +# +# Disclaimer: +# The author is not affiliated with Dr. Abernethy or Wes Fulton - CEO of +# Fulton Findings(TM) and author of the software package SuperSMITH +#------------------------------------------------------------------------------- +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . +# +# +-----------------------------------+ +# | execute this software with R: | +# | http://www.r-project.org/ | +# +-----------------------------------+ + +plotSingleDataSet <- function(x,isplotppp,...){ + if(isplotppp){ + # TODO: possibly, this does not allow much flexibility in plotting. + opadata <- modifyList(x$options,list(...)) + if(!is.null(x$data) && + !is.null(ti <- x$data$time) && + !is.null(ra <- x$data[,paste0("ppp.",tolower(opadata$ppos[1]))])){ + # TODO: add support for plotting all ppp columns, not just the first one + # this will need added support for vectors of col, lwd, lty, pch, ... +# if(opadata$log %in% c("x","xy","yx")){ +# replace +# } + t0 <- 0 + if(is.logical(opadata$threshold))if(opadata$threshold) + warning ("opadata$threshold is a logical value but numeric value was expected. Proceeding...") + if(is.numeric(opadata$threshold))t0 <- opadata$threshold + points(ti-t0,F0inv(ra,opadata$log),pch = opadata$pch, + col = opadata$col,lwd = opadata$lwd.points,cex=opadata$cex.points) + # option "log" should only be set and read from either + # the arguments of plot.abrem + # Other instances should be ignored + }else{ + stop("This Abrem object contains no probability plot positions.") + } + } +} \ No newline at end of file diff --git a/R/plotSingleFit.R b/R/plotSingleFit.R new file mode 100644 index 0000000..6893cad --- /dev/null +++ b/R/plotSingleFit.R @@ -0,0 +1,95 @@ +# R package 'abrem' +# Abernethy Reliability Methods +# Implementations of lifetime data analysis methods described in +# 'The New Weibull Handbook, Fifth edition' by Dr. Robert B. Abernethy. +# August 2014, Jurgen Symynck +# Copyright 2014, Jurgen Symynck +# +# For more info, visit http://www.openreliability.org/ +# +# For the latest version of this file, check the Subversion repository at +# http://r-forge.r-project.org/projects/abernethy/ +# +# Disclaimer: +# The author is not affiliated with Dr. Abernethy or Wes Fulton - CEO of +# Fulton Findings(TM) and author of the software package SuperSMITH +#------------------------------------------------------------------------------- +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . +# +# +-----------------------------------+ +# | execute this software with R: | +# | http://www.r-project.org/ | +# +-----------------------------------+ +# +plotSingleFit <- function(fit,v,opadata,...){ + opafit <- opadata + if(!is.null(fit$options)){ + opafit <- modifyList(opadata,fit$options) + } + opafit <- modifyList(opafit,list(...)) + if(opafit$is.plot.fit){ + t0 <- 0 + if(is.logical(opafit$threshold))if(opafit$threshold){ + if(is.logical(opadata$threshold)){if(opadata$threshold) + warning("opafit$threshold and opadata$threshold are logical values but numeric values were expected. Proceeding...") + }else{ + # reuse the t0 value from the data level + t0 <- opadata$threshold + } + } + if(is.numeric(opafit$threshold))t0 <- opafit$threshold + if(v >= 1)message("plotSingleFit: Adding fit...") + fitt0 <- 0 + cret <- NULL + if(!is.null(fit$t0)) fitt0 <- fit$t0 + if(!is.null(fit$beta) && !is.null(fit$eta)) + cret <- curve(F0inv(pweibull(x-fitt0+t0,fit$beta,fit$eta),opafit$log), + add=TRUE,n=1001, + # n=1001 is needed for displaying the extreme + # curvature towards -Inf with low Beta values + # like 0.1 + col=opafit$col,lwd=opafit$lwd,lty=opafit$lty, + xlim=getPlotRangeX(opafit$log), + log=opafit$log) + if(!is.null(fit$mulog) && !is.null(fit$sigmalog)) + cret <- curve(F0inv(plnorm(x-fitt0+t0,fit$mulog,fit$sigmalog),opafit$log), + add=TRUE,n=1001, + col=opafit$col,lwd=opafit$lwd,lty=opafit$lty, + xlim=getPlotRangeX(opafit$log), + log=opafit$log) + if(!is.null(cret)){ + cret$y[is.infinite(cret$y)] <- NA + # works for weibull canvas + cret$y[cret$y==0] <- NA + # replacing zero's is needed for lognormal canvas. + imin <- which.min(cret$y) + lines(rep(cret$x[imin],2), + y=c(cret$y[imin],getPlotRangeY(opafit$log)[1]), + col=opafit$col,lwd=opafit$lwd,lty=opafit$lty) + # plot vertical line towards -Inf + }else{ + if(opafit$verbosity >= 1)message("plotSingleFit: No distribution parameters available for plotting the fit...") + } +# if(!is.null(fit$rate)){ +# ### exponential ### +# if(opafit$verbosity >= 1)message( +# "plotSingleFit: Adding Exponential fit ...") +# curve(F0inv(pexp(x+t0,fit$rate),opafit$log),add=TRUE, +# col=opafit$col,lwd=opafit$lwd,lty=opafit$lty, +# xlim=getPlotRangeX(opafit$log), +# log=opafit$log) +# } + } + invisible() +} \ No newline at end of file diff --git a/R/plotSingleLegend.R b/R/plotSingleLegend.R new file mode 100644 index 0000000..fb06d70 --- /dev/null +++ b/R/plotSingleLegend.R @@ -0,0 +1,56 @@ +# R package 'abrem' +# Abernethy Reliability Methods +# Implementations of lifetime data analysis methods described in +# 'The New Weibull Handbook, Fifth edition' by Dr. Robert B. Abernethy. +# August 2014, Jurgen Symynck +# Copyright 2014, Jurgen Symynck +# +# For more info, visit http://www.openreliability.org/ +# +# For the latest version of this file, check the Subversion repository at +# http://r-forge.r-project.org/projects/abernethy/ +# +# Disclaimer: +# The author is not affiliated with Dr. Abernethy or Wes Fulton - CEO of +# Fulton Findings(TM) and author of the software package SuperSMITH +#------------------------------------------------------------------------------- +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . +# +# +-----------------------------------+ +# | execute this software with R: | +# | http://www.r-project.org/ | +# +-----------------------------------+ + +plotSingleLegend <- function(le,x,y){ + if(identical(label <- le$label,""))label <- NULL + if(is.null(le$legend))le$legend <- "" + legend( + x=x, + y=y, + legend=le$legend, + title=label, +# title.col=le$lcol, + cex = le$legend.text.size, + bg = "white", + lty = unlist(le$lty), + lwd = unlist(le$lwd), + pch = unlist(le$pch), + col = unlist(le$col), +# inset=0.1, + text.col = "black", + xpd=TRUE +# merge = TRUE + ) + # TODO: Warning: unlist coerces numeric colors to character! +} \ No newline at end of file diff --git a/R/plot_default_args.R b/R/plot_default_args.R new file mode 100644 index 0000000..60b58f4 --- /dev/null +++ b/R/plot_default_args.R @@ -0,0 +1,64 @@ +# R package 'abrem' +# Abernethy Reliability Methods +# Implementations of lifetime data analysis methods described in +# 'The New Weibull Handbook, Fifth edition' by Dr. Robert B. Abernethy. +# August 2014, Jurgen Symynck +# Copyright 2014, Jurgen Symynck +# +# For more info, visit http://www.openreliability.org/ +# +# For the latest version of this file, check the Subversion repository at +# http://r-forge.r-project.org/projects/abernethy/ +# +# Disclaimer: +# The author is not affiliated with Dr. Abernethy or Wes Fulton - CEO of +# Fulton Findings(TM) and author of the software package SuperSMITH +#------------------------------------------------------------------------------- +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . +# +# +-----------------------------------+ +# | execute this software with R: | +# | http://www.r-project.org/ | +# +-----------------------------------+ + + +plot_default_args <- function(){ + paronly <- c("ask","fig", "fin","lheight","mai", "mar", "mex", "mfcol", + "mfrow", "mfg","new","oma", "omd", "omi","pin", "plt", "ps", "pty", + "usr","xlog", "ylog","ylbias") + # parameters that can only be set using par() + # see $par() for the origin of this list + parreadonly <- c("xlog", "ylog", "adj", "ann", "ask", + "bg", "bty", "cex", "cex.axis", "cex.lab", + "cex.main", "cex.sub", "col", "col.axis", "col.lab", + "col.main", "col.sub", "crt", "err", "family", + "fg", "fig", "fin", "font", "font.axis", + "font.lab", "font.main", "font.sub", "lab", "las", + "lend", "lheight", "ljoin", "lmitre", "lty", + "lwd", "mai", "mar", "mex", "mfcol", + "mfg", "mfrow", "mgp", "mkh", "new", + "oma", "omd", "omi", "pch", "pin", + "plt", "ps", "pty", "smo", "srt", + "tck", "tcl", "usr", "xaxp", "xaxs", + "xaxt", "xpd", "yaxp", "yaxs", "yaxt", + "ylbias") + # par() parameter that can be set + # par(no.readonly=TRUE) + parplot <- unique(sort(c(parreadonly[!(parreadonly %in% paronly)], + "type","xlim","ylim","log","main","sub","xlab","ylab", + "ann","axes","frame.plot","panel.first","panel.last","asp"))) + # all valid (?) graphical parameters that can be supplied + # to plot.default + parplot +} \ No newline at end of file diff --git a/R/print.abrem.R b/R/print.abrem.R new file mode 100644 index 0000000..fa4dd3d --- /dev/null +++ b/R/print.abrem.R @@ -0,0 +1,59 @@ +# R package 'abrem' +# Abernethy Reliability Methods +# Implementations of lifetime data analysis methods described in +# 'The New Weibull Handbook, Fifth edition' by Dr. Robert B. Abernethy. +# August 2014, Jurgen Symynck +# Copyright 2014, Jurgen Symynck +# +# For more info, visit http://www.openreliability.org/ +# +# For the latest version of this file, check the Subversion repository at +# http://r-forge.r-project.org/projects/abernethy/ +# +# Disclaimer: +# The author is not affiliated with Dr. Abernethy or Wes Fulton - CEO of +# Fulton Findings(TM) and author of the software package SuperSMITH +#------------------------------------------------------------------------------- +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . +# +# +-----------------------------------+ +# | execute this software with R: | +# | http://www.r-project.org/ | +# +-----------------------------------+ + +print.abrem <- function(x,...){ + if(!is.null(x$data)){ + cat("\n$data :\n") + print(x$data) + } + if(!is.null(x$n)) cat("\n$n :",x$n) + if(!is.null(x$fail)) cat("\n$fail :",x$fail) + if(!is.null(x$susp)) cat("\n$susp :",x$susp,"\n") + # TODO: check what to do with mentionin this data at the x$fit[[i]] level + + + if(!is.null(x$options)){ + cat("\n$options\n ") + if((le <- length(x$options)) >=1){ + cat(paste0(" named list of ",le," items.\n")) + }else{ + cat(" empty list.\n") + } + } + if(!is.null(x$fit)){ + cat("\n$fit(s)\n ") + print(x$fit) + } + invisible(x) +} diff --git a/R/seq.R b/R/seq.R new file mode 100644 index 0000000..56d3474 --- /dev/null +++ b/R/seq.R @@ -0,0 +1,43 @@ +# R package 'abrem' +# Abernethy Reliability Methods +# Implementations of lifetime data analysis methods described in +# 'The New Weibull Handbook, Fifth edition' by Dr. Robert B. Abernethy. +# August 2014, Jurgen Symynck +# Copyright 2014, Jurgen Symynck +# +# For more info, visit http://www.openreliability.org/ +# +# For the latest version of this file, check the Subversion repository at +# http://r-forge.r-project.org/projects/abernethy/ +# +# Disclaimer: +# The author is not affiliated with Dr. Abernethy or Wes Fulton - CEO of +# Fulton Findings(TM) and author of the software package SuperSMITH +#------------------------------------------------------------------------------- +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . +# +# +-----------------------------------+ +# | execute this software with R: | +# | http://www.r-project.org/ | +# +-----------------------------------+ + +seq.wb <- function(from,to,base=seq(1,9,1)){ + # define gridline positions for y axis + r <- c(seq.log(from,.9,base),rev(1-seq.log(1-to,0.1,base))[-1]) + r[r >= from & r<=to]} + +seq.log <- function(from,to,base=c(1,2,5)){ + r <- NULL + for(i in floor(log10(from)):floor(log10(to)))r <- c(r,base*10^i) + r[r >= from & r<=to]} diff --git a/R/splitargs.R b/R/splitargs.R new file mode 100644 index 0000000..32059af --- /dev/null +++ b/R/splitargs.R @@ -0,0 +1,45 @@ +# R package 'abrem' +# Abernethy Reliability Methods +# Implementations of lifetime data analysis methods described in +# 'The New Weibull Handbook, Fifth edition' by Dr. Robert B. Abernethy. +# August 2014, Jurgen Symynck +# Copyright 2014, Jurgen Symynck +# +# For more info, visit http://www.openreliability.org/ +# +# For the latest version of this file, check the Subversion repository at +# http://r-forge.r-project.org/projects/abernethy/ +# +# Disclaimer: +# The author is not affiliated with Dr. Abernethy or Wes Fulton - CEO of +# Fulton Findings(TM) and author of the software package SuperSMITH +#------------------------------------------------------------------------------- +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . +# +# +-----------------------------------+ +# | execute this software with R: | +# | http://www.r-project.org/ | +# +-----------------------------------+ + +splitargs <- function(...){ + arg <- list(...) + argnames <- names(arg) + parplot <- plot_default_args() + ret <- list() + opanames <- names(options.abrem()) + ret$opa <- arg[tolower(argnames) %in% tolower(opanames)] + # ret$opa can be an emply list, which is compatible with modifyList() + ret$rem <- arg[!(tolower(argnames) %in% tolower(opanames))] + ret +} diff --git a/data/abrem13.rda b/data/abrem13.rda new file mode 100644 index 0000000..4daee94 Binary files /dev/null and b/data/abrem13.rda differ diff --git a/data/abrem13c1.rda b/data/abrem13c1.rda new file mode 100644 index 0000000..6ac7c2b Binary files /dev/null and b/data/abrem13c1.rda differ diff --git a/data/abrem13c2.rda b/data/abrem13c2.rda new file mode 100644 index 0000000..f1e2724 Binary files /dev/null and b/data/abrem13c2.rda differ diff --git a/data/abrem13c3.rda b/data/abrem13c3.rda new file mode 100644 index 0000000..0da18bc Binary files /dev/null and b/data/abrem13c3.rda differ diff --git a/data/abrem13c4.rda b/data/abrem13c4.rda new file mode 100644 index 0000000..8af176e Binary files /dev/null and b/data/abrem13c4.rda differ diff --git a/data/abrem2.rda b/data/abrem2.rda new file mode 100644 index 0000000..012caa6 Binary files /dev/null and b/data/abrem2.rda differ diff --git a/data/abrem3.rda b/data/abrem3.rda new file mode 100644 index 0000000..80629c3 Binary files /dev/null and b/data/abrem3.rda differ diff --git a/data/abrem_mix1.rda b/data/abrem_mix1.rda new file mode 100644 index 0000000..e133a81 Binary files /dev/null and b/data/abrem_mix1.rda differ diff --git a/demo/00Index b/demo/00Index new file mode 100644 index 0000000..cad62e8 --- /dev/null +++ b/demo/00Index @@ -0,0 +1 @@ +abrem calculating fits from (life-)time observation data. diff --git a/demo/abrem.r b/demo/abrem.r new file mode 100644 index 0000000..b7ddb64 --- /dev/null +++ b/demo/abrem.r @@ -0,0 +1,191 @@ +objectlist <- ls() # for removing unused demo objects later +abrem.defaults <- options.abrem() +# +----------------------------------------------------------------------------- +options.abrem(sub=paste0( + "R package abrem ",packageVersion("abrem"), + " (https://r-forge.r-project.org/projects/abernethy)"), + method.fit=c("rr","xony"),method.conf.blives="mcpivotals") +# +----------------------------------------------------------------------------- +lto <- runif(8,100,1000) +da <- list(Abrem(lto,ppos="Benard",col="black",pch=1), + Abrem(lto,ppos="beta",col="blue",pch=2), + Abrem(lto,ppos="mean",col="green",pch=3), + Abrem(lto,ppos="KM",col="yellow3",pch=4), + Abrem(lto,ppos="Hazen",col="orange",pch=5), + Abrem(lto,ppos="Blom",col="red",pch=6)) +da <- abrem.fit(da) +plot.abrem(da,main='Comparing different ranking methods.') +#readline(prompt = "Hit to see next plot.\n") +# +----------------------------------------------------------------------------- + +da <- Abrem(runif(8,100,1000)) +da <- abrem.fit(da) +da <- abrem.conf(da) +da <- abrem.conf(da,method.conf.blives="bbb",lty=2) +plot(da) +#readline(prompt = "Hit to see next plot.\n") +# +----------------------------------------------------------------------------- +if(require(boot)){ + data(aircondit,aircondit7,package="boot") + da1 <- abrem.fit( + Abrem(time=aircondit$hours,label='\"aircondit\" dataset (package \"boot\")')) + da2 <- abrem.fit( + Abrem(time=aircondit7$hours,label='\"aircondit7\" dataset',col="red",pch=4)) + plot.abrem(list(da1,da2),xlim=c(0.01,1e5), + main='\"aircondit\" dataset (package \"boot\")', + xlab="Time To Failure (hours)") +}else message("Recommended package \"boot\" is not available, skipping this demo...") +# +----------------------------------------------------------------------------- +if(require(boot)){ + data(hirose,package="boot") + names(hirose) <- c("volt","time","event") + da <- list() + hc <- rev(rainbow(4,end=4/6)) + # create temperature-like colors + da[[1]] <- abrem.fit(Abrem(subset(hirose,volt==5), + label="Voltage= 5 [V]",col=hc[1],pch=0)) + da[[2]] <- abrem.fit(Abrem(subset(hirose,volt==7), + label="Voltage= 7 [V]",col=hc[2],pch=1)) + da[[3]] <- abrem.fit(Abrem(subset(hirose,volt==10), + label="Voltage= 10 [V]",col=hc[3],pch=2)) + da[[4]] <- abrem.fit(Abrem(subset(hirose,volt==15), + label="Voltage= 15 [V]",col=hc[4],pch=5)) + + def <- options.abrem() + options.abrem(legend.text.size=0.5,xlim=NULL) + + #lapply(da,plot.abrem) + plot.abrem(da,xlim=c(1,1e6),main='\"Hirose\" dataset (package \"boot\")') +# legend("topright",c("5 [V]","7 [V]","10 [V]","15 [V]"), +# title="Test voltages:",bg="white",inset=0.05, +# cex=0.7,col=hc,lty=1,lwd=2) + invisible(options.abrem(def)) +}else message("Recommended package \"boot\" is not available, skipping this demo...") + +# +----------------------------------------------------------------------------- +if(require(MASS)){ + def <- options.abrem() + data(motors,package="MASS") + names(motors) <- c("temp","time","event") + mo <- list() +# mo[[4]] <- NULL #abrem.fit(Abrem(subset(motors,temp==150), +# #label="Temp = 150")) + # note that subset(motors,temp==150) has no failures, + # so it cannot be analysed by this package yet. + hc <- rev(rainbow(3,end=4/6)) + # create temperature-like colors + mo[[1]] <- abrem.fit(Abrem(subset(motors,temp==170), + label="Temp = 170",col=hc[1],pch=0)) + mo[[2]] <- abrem.fit(Abrem(subset(motors,temp==190), + label="Temp = 190",col=hc[2],pch=1)) + mo[[3]] <- abrem.fit(Abrem(subset(motors,temp==220), + label="Temp = 220",col=hc[3],pch=2)) + + mo <- lapply(mo,abrem.conf) + +# mo[[2]] <- abrem.conf(abrem.fit(Abrem(time=mo[[2]]$time,event=mo[[2]]$status), +# col=hc[1],pch=0)) +# mo[[3]] <- abrem.conf(abrem.fit(Abrem(time=mo[[3]]$time,event=mo[[3]]$status), +# col=hc[2],pch=1)) +# mo[[4]] <- abrem.conf(abrem.fit(Abrem(time=mo[[4]]$time,event=mo[[4]]$status), +# col=hc[3],pch=2)) + plot.abrem(mo,xlim=c(5,1e6),main='\"Motors\" dataset (package \"MASS\")') + options.abrem(def) +}else message("Recommended package \"MASS\" is not available, skipping this demo...") + +# +----------------------------------------------------------------------------- +da <- Abrem(rweibull(8,1.5,1000)+600) + # store any threshold value for later use while plotting +da <- abrem.fit(da,dist="weibull2p",col="blue") +da <- abrem.fit(da,dist="weibull3p",col="steelblue") +da <- abrem.fit(da,dist="lognormal2p",col="red") +da <- abrem.fit(da,dist="lognormal3p",col="orange") + +plot(da,main="Comparison between weibull and lognormal, two and three parameter") +legend("topleft",legend="Weibull plotting canvas",bg="white") +plot(da,main="Comparison between weibull and lognormal, two and three parameter",log="xy") +legend("topleft",legend="Lognormal plotting canvas",bg="white") +# +----------------------------------------------------------------------------- +da <- abrem.fit(Abrem(rweibull(5,3,1000)),dist="weibull") +message("This might take some time...") +for(i in 1:20){ + da <- abrem.conf(da,col="blue",S=100,lwd=1,in.legend=i==TRUE) + # just allow the first confidence bound in the legend, drop the rest +} +for(i in 1:20){ + da <- abrem.conf(da,col="red",S=1e4,lwd=1,in.legend=i==TRUE) + # just allow the first confidence bound in the legend, drop the rest +} + +plot(da,main='Variability in \"mcpivotals\" B-life confidence for S=100 and 1e4') +# +---------------------------------------------------------------------------- +data(abrem_mix1) +earlyda <-abrem_mix1[1:10] +midda <-abrem_mix1[11:131] +endda <-abrem_mix1[132:200] + +for(th in c(FALSE,TRUE)){ + da <-Abrem(abrem_mix1,col="gray",pch=1, + label=" abrem_mix1 (Complete, unaltered dataset) ") + da21 <-Abrem(fail=earlyda,susp=c(midda,endda),col="black",pch=19) + da22 <-Abrem(fail=midda,susp=c(earlyda,endda),col="blue",pch=3) + da23 <-Abrem(fail=endda,susp=c(earlyda,midda),col="green3",pch=4) + + ### without threshold parameter corrected plotting ### + da21 <- abrem.fit(da21, + label='\"Early\" data',dist="weibull2p",threshold=th) + da22 <- abrem.fit(da22, + label='\"Mid-\" data',dist="weibull3p",threshold=th) + da23 <- abrem.fit(da23, + label='\"End\" data',dist="weibull3p",threshold=th) + plot.abrem(list(da,da21,da22,da23),xlim=c(0.5,1e5), + main=paste0("Diaphragm life data of acid gas compressor, ", + ifelse(th,"with","without"), + " threshold parameter corrected plotting")) +} + +# +---------------------------------------------------------------------------- +data(abrem_mix1) +earlyda <-abrem_mix1[1:10] +midda <-abrem_mix1[11:131] +endda <-abrem_mix1[132:200] +da <-Abrem(abrem_mix1,col="gray",pch=1, + label=" abrem_mix1 (Complete, unaltered dataset) ") +da21 <-Abrem(fail=endda,susp=c(earlyda,midda),col="green2",pch=19) +da22 <-Abrem(fail=endda,susp=c(earlyda,midda),col="green3",pch=3) +da23 <-Abrem(fail=endda,susp=c(earlyda,midda),col="green4",pch=4) +da21 <- abrem.fit(da21, + label='\"End\" data, threshold=FALSE',dist="weibull3p",threshold=FALSE) +da22 <- abrem.fit(da22, + label='\"End\" data, threshold=TRUE',dist="weibull3p",threshold=TRUE) +da23 <- abrem.fit(da23, + label='\"End\" data, threshold=5000',dist="weibull3p",threshold=1800) +plot.abrem(list(da,da21,da22,da23),xlim=c(0.5,1e5), + main="Diaphragm life data of acid gas compressor, threshold parameter usage") + +## +---------------------------------------------------------------------------- +fail1 <- c(535,269,1066,788,288) +fail2 <- c(1043,843,1303,1847,2590) +da1 <- abrem.fit(Abrem(time=fail1,pch=0,col="blue")) +da2 <- abrem.fit(Abrem(time=fail2,pch=4,col="red")) + +da1 <- abrem.conf(da1,method.conf.blives="lrb",cl=0.5, + in.legend.blives=F,is.plot.cb=F,lty=5) +da1 <- abrem.conf(da1,method.conf.blives="lrb",cl=0.9) +da1 <- abrem.conf(da1,method.conf.blives="lrb",cl=0.95, + in.legend.blives=F,is.plot.cb=F,lty=5) +da2 <- abrem.conf(da2,method.conf.blives="lrb",cl=0.5, + in.legend.blives=F,is.plot.cb=F,lty=5) +da2 <- abrem.conf(da2,method.conf.blives="lrb",cl=0.9) +da2 <- abrem.conf(da2,method.conf.blives="lrb",cl=0.95, + in.legend.blives=F,is.plot.cb=F,lty=5) +par(mfrow=c(1,2)) +dataset <- list(da1,da2) +plot.abrem(dataset,main="Weibull 2p") +contour.abrem(dataset) +# +---------------------------------------------------------------------------- + +invisible(options.abrem(abrem.defaults)) + # resetting the all abrem options. +rm(list = ls()[!(ls() %in% objectlist)]) + # removing demo objects diff --git a/man/Abrem.Rd b/man/Abrem.Rd new file mode 100644 index 0000000..c8bb8ec --- /dev/null +++ b/man/Abrem.Rd @@ -0,0 +1,150 @@ +\name{Abrem} +\alias{Abrem} +\alias{ppos} +\alias{ppp} +%- Also NEED an '\alias' for EACH other topic documented here. +\title{ + Create an \code{abrem} Object for Lifetime and Reliability Analysis +} +\description{ + This function creates an object of class \code{"abrem"} for further processing + by the other functions of \pkg{abrem}. +} +\usage{Abrem(x,\dots)} +\arguments{ + \item{x}{ + Either a dataframe containing at least \code{$time} and \code{$event} + columns, or a vector of class \code{"numeric"} or \code{"integer"} with + (life-)time observations. + + See section "Details" for other data passing arguments. + } + \item{\dots}{ + Graphical options for plotting the \code{abrem} object; see section "Details". + %% Extend this with "type" for type of censoring + } +} +\details{ + There are several methods to passing arguments for building an \code{abrem} + object. + \itemize{ + \item When a single unnamed vector of class \code{"numeric"} + or \code{"integer"} is supplied, it is treated as a vector + of (life-)time observations. + \item If argument \code{time} or \code{fail} is provided, it is treated as + a vector of (life-)time observations. Take care \emph{not} to supply both + \code{time} and \code{fail} in the same function call. + \item If argument \code{event} is provided, it is treated as + a vector of event indicators with possible values of + \code{0} and \code{1}. See section "Value" for more details on + event vectors. + \item If argument \code{susp} is provided, it is treated as + a vector of right-censored (life-)time observations (also called + suspended observations or suspensions). + In that case, argument \code{time} or \code{fail} is mandatory and is + treated as a vector of failure times. + \item If argument \code{x} is of class \code{"data.frame"}, + then it should at least contain \code{$time} and \code{$event} + columns. Additional columns in the dataframe will be reused in the + \code{abrem} object, allowing for extra information like + serial numbers to be included (see section "Examples"). + } + + \code{Abrem} \emph{always} generates probability plotting positions for graphically + displaying the (life-)time observations and for later usage + by \code{\link{abrem.fit}} when using rank regression. + The type of plotting positions to be calculated (also quite confusingly + called "ranks") is controlled by argument \code{ppos}, a vector of + class \code{"character"}. Currently, all ranking calculations are done by + function \code{\link[abremPivotals]{getPPP}} from + package \pkg{abremPivotals}. Supported ranking methods include: + \itemize{ + \item \code{"beta"}: Rankin using the incomplete beta function, also called "exact median ranks". + \item \code{"benard"}: A good approximation to exact median ranks, + and currently the default method. + \item \code{"mean"}, also known as Herd-Johnson. + \item \code{"km"}: Kaplan-Meier ranking with modification for final complete failure. + \item \code{"hazen"} or \emph{modified} Kaplan-Meier. + \item \code{"blom"}. + } + + All methods can be passed in the \code{ppos} argument vector + at the same time but currently, only the first element in the \code{ppos} + vector will be used for further fitting by \code{\link{abrem.fit}}. + + Note that is is currently allowed to have \code{NA} values in the + \code{time} argument. In that case, the vector is expected to be + ordered and no ordering will be applied by \code{Abrem}. This feature is + useful in combination with the output of \code{\link{params.to.ob}}. + + Additionally, one can supply any options available from \code{options.abrem}, + such as \code{col} or \code{is.plot.legend}. Some of these options + will be used when plotting the (life-)time observations using \code{plot.abrem}. + Subsequent calls to \code{abrem.fit} and \code{abrem.conf} will inherit these options. + + % \code{\link[abremout:plot.abrem]{plot.abrem}}. +} +\value{ + A named list of class \code{"abrem"}. The first list + item (\code{$data}) is a dataframe with at least three columns: + \describe{ + \item{\code{$time}}{ + An ordered vector with (life-)time observations. + } + \item{\code{$event}}{ + A vector of class \code{"numeric"} with right-censoring indicators. + Values of \code{1} mean "dead" or "failed" while \code{0} + mean "alive" or "right-censored"/"suspended" observations. + This censoring indicator scheme is modeled after the + \code{Surv} function of the \pkg{survival} + package. + %% add support for factors + } + \item{\code{$ppp. \dots }}{ + Depending on the argument \code{\link{ppos}} (either passed as + an argument or taken from \code{options.abrem}), a vector + of class \code{"numeric"} with the probability plot positions' y-coordinates. + The exact column name depends on the selected ranking method. + + Defaults to \code{$ppp.benard}. + } + } +} +\author{Jurgen Symynck \email{jusy@openreliability.org}} +\examples{ +## These code lines all generate the same object ## +Abrem(c(500,1200,900,1300,510)) +Abrem(time=c(500,1200,900,1300,510)) +Abrem(time=c(500,1200,900,1300,510),event=c(1,1,1,1,1)) +Abrem(fail=c(500,1200,900,1300,510)) +Abrem(fail=c(500,1200,900,1300,510),susp=c()) +da1 <- data.frame( + serial=c("S12","S16","S17","S3","S5"), + time=c(500,1200,900,1300,510), + event=c(1,1,1,1,1)) +Abrem(da1,label="complete dataset",pch=1) +da1 <- Abrem(da1,label="complete dataset",pch=3,col="orange3") + +## Generate a similar dataset, but with suspensions ## +Abrem(time=c(500,1200,900,1300,510),event=c(1,1,1,0,0)) +Abrem(data.frame(time=c(500,1200,900,1300,510),event=c(1,1,1,0,0))) +Abrem(fail=c(500,1200,900),susp=c(1300,510)) +Abrem(time=c(500,1200,900),susp=c(1300,510)) +da2 <- Abrem(fail=c(500,1200,900,1300,510), + event=c(1,1,1,0,0),label="censored dataset",pch=1,col="blue") + +## plot datasets ## +plot.abrem(list(da1,da2)) + +## different ranking methods ## +## note that ppos is implemented case insensitive ## +lto <- runif(8,100,1000) +da3 <- list(Abrem(lto,ppos="Benard",col="black",pch=1), + Abrem(lto,ppos="beta",col="blue",pch=2), + Abrem(lto,ppos="mean",col="green",pch=3), + Abrem(lto,ppos="KM",col="yellow3",pch=4), + Abrem(lto,ppos="Hazen",col="orange",pch=5), + Abrem(lto,ppos="Blom",col="red",pch=6)) +da3 <- abrem.fit(da3) +plot.abrem(da3,main='Comparing different ranking methods.') +} diff --git a/man/abrem-package.Rd b/man/abrem-package.Rd new file mode 100644 index 0000000..9a99996 --- /dev/null +++ b/man/abrem-package.Rd @@ -0,0 +1,72 @@ +\name{abrem-package} +\alias{abrem-package} +\alias{abrem} +\docType{package} +\title{ + Abernethy Reliability Methods +} +\description{ + Implementation of reliability methods presented in \dQuote{The New Weibull Handbook}, Fifth Edition by Dr. Robert B. Abernethy. This package is dedicated to the control and view of models developed in dependant packages abremPivotals and abremDebias (under construction) implemented using the R object model. It is based on the now outdated + and unmaintained \pkg{weibulltoolkit} package, still available from: + + http://sourceforge.net/projects/weibulltoolkit/. + % TODO: replace with Jacob's version from DESCRIPTION +} +\author{ +Jurgen Symynck + +Maintainer: Jurgen Symynck +} +\note{} +\references{ + \itemize{ + \item + \url{http://www.openreliability.org/} + Homepage of the Abernethy Reliability Methods project. + \item + Robert B. Abernethy, + \emph{The New Weibull Handbook, Fifth Edition} + (Robert B. Abernethy, 2004) + \item + Jerald F. Lawless, + \emph{Statistical Models and Methods for Lifetime Data, + 2nd edition} + (Wiley-Interscience, Hoboken N.J., 2003) + \item William Q. Meeker and Luis A. Escobar, + \emph{Statistical Methods for Reliability Data} + (Wiley-Interscience, New York, 1998) + + \item Chi-Chao Lui, + \emph{A Comparison Between The Weibull And Lognormal Models Used + To Analyse Reliability Data} + (dissertation from University of Nottingham, 1997) +% \item +% Efron, B. and Tibshirani, R., +% \emph{An Introduction to the Bootstrap} +% (Chapman & Hall, 1993) + \item Jurgen Symynck, Filip De Bal, + \emph{Weibull analysis using R, in a nutshell} + (New Technologies and Products in Machine Manufacturing Technology, + Stefan cel Mare University of Suceava, 2010). + +% Available on +%\url{http://mechanics.kahosl.be/fatimat/index.php/downloads-and-information/40/171}. + \item Jurgen Symynck, Filip De Bal, + \emph{Monte Carlo pivotal confidence bounds for Weibull analysis, + with implementations in R} (New Technologies and Products in + Machine Manufacturing Technology, Stefan cel Mare University of + Suceava, 2011). + + Available on +\url{http://www.fim.usv.ro/conf_1/tehnomusjournal/pagini/journal2011/files/7.pdf} +%and \url{http://mechanics.kahosl.be/fatimat/index.php/downloads-and-information/40/198}. + + \item + \emph{SuperSMITH(TM) Weibull}, + \url{http://www.Weibullnews.com/contents.htm} + \item \url{http://r-forge.r-project.org/projects/abernethy/} + \item \url{http://sourceforge.net/projects/weibulltoolkit/} + } +} + +\keyword{ package } diff --git a/man/abrem.conf.Rd b/man/abrem.conf.Rd new file mode 100644 index 0000000..b0a4bc7 --- /dev/null +++ b/man/abrem.conf.Rd @@ -0,0 +1,174 @@ +\name{abrem.conf} +\alias{abrem.conf} +\alias{contours} +\alias{MLE contours} +\alias{cl} +%% cl or CL +\alias{unrel.n} +\alias{conf.what} +\alias{method.conf.blives} +\alias{conf.blives.sides} +\alias{S} +\alias{in.legend} + +\title{Add Confidence to \code{abrem} Objects} +\description{ + This function adds confidence calculations to + various entities in \code{abrem} objects. +} +\usage{abrem.conf(x,which="all",\dots)} +\arguments{ + \item{x}{Object of class \code{"abrem"}.} + \item{which}{Calculate which fit in the \code{abrem} object will be processed.} + \item{\dots}{Options for calculating confidence, and for plotting the results.} +} +\details{ + This function adds confidence calculations to various entities in + \code{abrem} objects and adds them to the object alongside any pre-existing + confidence calculations. + + Additional options for calculating B-life confidence are passed with: + + \describe{ + \item{\code{cl}}{ + Confidence level: A single number from the interval \code{[0,[1} + specifying the confidence level for various confidence calculations. + + Defaults to \code{0.9}. + } + \item{\code{conf.blives.sides}}{ + Either \code{"lower"}, \code{"upper"} or \code{"double"}, + specifying the type of bound(s) to be calculated. + + Defaults to \code{c("double")}, the other options are currently + not implemented. + } + \item{\code{unrel.n}}{ + An integer controlling the amount of unreliability levels for + which B-life confidence bounds are calculated and ultimately plotted. + + Higher numbers will result in smoother confidence bounds. In any + case, confidence intervals will be calculated for: + \itemize{ + \item the B-lives at unreliability levels specified with option + \code{\link{unrel}} + \item the B-life at \code{50 [\%]} unreliability + \item the B-life at the calculcate characteristic life + or logmean (depending on the fitted distribution) + } + + Note: When plotting fits and confidence bounds that are adjusted with + a threshold (see option \code{"threshold"}), it is often the case that + the bounds appear to be cut of on the left. This can be countered by + dramatically increasing \code{unrel.n}, resulting in confidence + bounds that extend to the edge of the plotting area. + + Defaults to \code{25}. + } + \item{\code{conf.what}}{ + A vector of class \code{"character"} describing for which entities + that confidence should be calculated. + + Defaults to \code{c("blives")}, the only type currently supported. + } + \item{\code{unrel}}{ + An unordered numeric vector with unreliability levels for which + B-life confidence will be calculated. + + Defaults to \code{c(0.1,0.05,0.01)}. + } + + \item{\code{method.conf.blives}}{ + A vector of class \code{"character"} describing the technique to be + used for calculating confidence for B-lives. Possible values are + \code{"bbb"} (Beta Binomial confidence bounds), + \code{"lrb"} (Likelihood Ratio confidence bounds) and + \code{"mcpivotal"} or \code{"mcpivotals"} (Monte Carlo Pivotal + confidence bounds). + + Monte Carlo Pivotal confidence bounds use a large number of + simulations to calculate the confidence bounds. See option + \code{"S"} for more info. + + Defaults to \code{c("mcpivotals")}. + } + \item{\code{S}}{ + An integer describing the number of Monte Carlo simulations on + which the Monte Carlo pivotal confidence bounds and calculation + of the "prr" goodness-of-fit indicator are based. + + High values are needed for good confidence bounds at the lower + end of the fitted model, especially for data with heavy censoring. + + Note that \code{S >= 100} and that \code{S} must be divisible by 10. + + Defaults to \code{10000}. + } + \item{\code{in.legend}}{ + Logical value controlling the inclusion of confidence calculation + results in the legend. + + If \code{in.legend=FALSE} is passed , + the resulting confidence calculations will be omitted from the legend. + + Defaults to \code{TRUE}. + } + } + + Additionally, one can pass any options available from \code{options.abrem}, + such as \code{col} or \code{is.plot.legend}. The graphical options + will be used when plotting the (life-)time observations using \code{plot.abrem}. +} +\value{ + The function returns its argument \code{x}, extended with the confidence + calculations and any optional graphical and calculation arguments + as passed to the function. +} +\author{Jurgen Symynck \email{jusy@openreliability.org}} +\note{ + \itemize{ + \item Currently, only \code{which = "all"} is supported, meaning that a + call to \code{abrem.conf} attempts calculation of confidence for all + fits in the \code{abrem} object. + \item Currently, only \code{conf.what = "blives"} and + \code{conf.blives.sides = "double"} are supported. + \item Currently, the Monte Carlo + pivotal confidence bounds are identical to superSMITH's + MC pivotal bounds for complete, uncensored data. For heavily censored + datasets with few failures however, the bounds appear more optimistic than + superSMITH's bounds. Research on this issue is ongoing. + } +} +\seealso{ + \code{\link{Abrem}}, + \code{\link{abrem.fit}}, + \code{\link{contour.abrem}} +} +\examples{ +## full dataset ## +da1 <- Abrem(runif(10,100,1e4),label="Complete data") +da1 <- abrem.fit(da1) +da1 <- abrem.conf(da1,method.conf.blives="mcpivotals",col="red") +da1 <- abrem.conf(da1,method.conf.blives="bbb",col="orange") +da1 <- abrem.conf(da1,method.conf.blives="lrb",col="yellow3") +print(da1$fit[[1]]$conf$blives[[1]]) +plot(da1,main="Comparison between MC Pivotal bounds and BB Bounds") + +## censored dataset: generates a warning for MC Pivotal confidence bounds ## +da2 <- runif(8,100,1e4) +da2 <- Abrem(fail=da2,susp=rep(max(da2),2),label="Type II censored data") + # generate a 'type 2' censored dataset +da2 <- abrem.fit(da2) +da2 <- abrem.conf(da2,method.conf.blives="mcpivotals",col="blue1") +da2 <- abrem.conf(da2,method.conf.blives="bbb",col="steelblue") +da2 <- abrem.conf(da2,method.conf.blives="lrb",col="cyan3") +plot(da2,main="Comparison between different bound types.") + +## show variability in Monte Carlo Pivotal bounds with low S ## +da3 <- Abrem(rweibull(5,3,1000)) +da3 <- abrem.fit(da3) +for(i in 1:20) da3 <- abrem.conf(da3,S=300,lwd=1,col="red") + # just keep adding bounds to the abrem object... +plot(da3,is.plot.legend=FALSE, + main="Variability in MC Pivotal Conf. Bounds for S=300") +} diff --git a/man/abrem.fit.Rd b/man/abrem.fit.Rd new file mode 100644 index 0000000..a9db0b5 --- /dev/null +++ b/man/abrem.fit.Rd @@ -0,0 +1,213 @@ +\name{abrem.fit} +\alias{abrem.fit} +\alias{dist} +\alias{method.fit} +\alias{prr} +\alias{P-value} +\alias{AbPval} +\alias{threshold} +\alias{in.legend} + +\title{ + Fit Distributions to \code{Abrem} Objects +} +\description{ + This function fits probability distributions to \code{Abrem} objects. +} +\usage{abrem.fit(x,\dots)} +\arguments{ + \item{x}{ + Object of class \code{"abrem"}. + } + \item{\dots}{ + Options for fitting the (life-)time observations, + and for plotting the results. + } +} +\details{ + This function calculates fits for the (life-)time observations in the + \code{abrem} object and adds them to the object alongside any + pre-existing fits. + + Fitting options are passed with the \code{dist} and \code{method.fit} + arguments: + + \describe{ + \item{\code{dist}}{ + A character string with the target distribution for fitting. + Possible values are \code{"weibull"} or \code{"weibull2p"}, + \code{"weibull3p"} (three-parameter Weibull), \code{"lognormal"} + or \code{"lognormal2p"} and \code{"lognormal3p"} (three-parameter lognormal). + + Defaults to \code{"weibull"}. + } + \item{\code{in.legend}}{ + Logical value controlling the inclusion of various elements in + the legend. + + If \code{in.legend=FALSE} is passed, + the resulting fit calculations will be omitted from the legend, + leaving only observation summary data. + + Defaults to \code{TRUE}. + } + \item{\code{method.fit}}{ + A vector of class \code{"character"} with fitting options. + + Defaults to \code{c("rr","xony")}. + + \itemize{ + \item \code{"rr"}: Rank Regression (RR). The parameters are calculated + using least square linear regression on the probability plot positions that + were calculated from the (life-)time observations during the creation of the + \code{Abrem} object (see option \code{\link{ppos}} and + function \code{\link{Abrem}}). The \code{"rr"} argument invokes + the \code{\link[stats]{lm}} function of the \pkg{stats} package, + and its return object can be retrieved from the \code{abrem} + object under \code{x$fit[[i]]$lm} where \code{x} is the return + object of \code{abrem.fit} and \code{i} the number of the fit + in the abrem object. + + By passing \code{"rr2"}, one can explicitly use the fast C++ compiled + linear regression code from \pkg{abremPivotals} that is also + used for the Monte Carlo pivotal confidence bound calculation. + + If Rank Regression is used then it + is \emph{mandatory} to additionally specify either X-on-Y or + Y-on-X regression. + + \item \code{"xony"},\code{"yonx"}: Differentiate between X-on-Y and Y-on-X + regression, respectively. For rank regression in lifetime analysis, + it is best practice to use the X values ((life-)time observations) + as the response variables whose horizontal distance to the fit line + must be minimized, and the Y values (unreliabilities) as the + explanatory variable. + + \item \code{"mle"}: Maximum Likelihood Estimation (MLE), using + many functions of the \pkg{debias} package. + + \item \code{"mle-rba"}: Maximum Likelihood Estimation with Reduced Bias + based on the median bias of the distributions MLE. + This options uses functions \code{\link[debias]{RBAbeta}} and + \code{\link[debias]{RBAsigma}}of the + \pkg{debias} package. + } + + Similar to Rank Regression, by passing \code{"mle2"} or + \code{"mle2-rba"} one can explicitly use fast C++ compiled code + from the \pkg{debias} package. Support for these extra + implementations (in the case of MLE) might be removed in + future versions. + } + + Additionally, one can pass any options available from \code{options.abrem}, + such as \code{col} or \code{is.plot.legend}. The graphical options + will be used when plotting the (life-)time observations using \code{plot.abrem}. + Subsequent calls to \code{abrem.conf} will inherit these options. + + The "prr" goodness-of-fit indicator is also calculated here when appropriate. + See the "Examples" section on how to retrieve it, see the "References" + section for additional information on the prr, pve and AbPval values. + % \code{\link[abremout:plot.abrem]{plot.abrem}}. + + Passing the \code{threshold} parameter here will result in plotting the + fit (and its associated probability plotting positions) with a + threshold value subtracted. + + If three-parameter models are used and + \code{threshold = TRUE}, the calculated third parameter (\code{t0}) of + the \emph{last three-parameter fit} in the abrem object will be used for + plotting the graphical entities. If a numeric value is passed, then + only the value from the \emph{last} call to \code{abrem.fit} is used. + % what if different t0 values are given with different fits? this should + % result in only one fit used. + % test the statement that the t0 of the last three-parameter fit + % is used + + Currently, there is no graceful error recovery after attempting to fit + lifetime data including negative time observations, for example + \code{abrem.fit(Abrem(-5:10)).} + + } +} +\value{ + The function returns its argument \code{x}, extended with the + calculated fit and the optional graphical and calculation arguments as + provided to the function. +} +\author{Jurgen Symynck \email{jusy@openreliability.org}} +\note{ + The \code{$time} column of the \code{abrem} object can have \code{NA} + values when accompanied with the 'suspended' indicator in the \code{$event} + column (= \code{0}). While this poses no problem for \code{abrem.fit} when + using Rank Regression by means of \code{"rr"}, it results in + different fits when using the Maximum Likelihood Estimantion method of + \code{\link[survival]{survreg}}. Also, \code{"rr2"} currently doesn't + support \code{NA} values fully. +} +%\section{To Do}{ +% \itemize{ +% \item Research the effect of \code{NA} values in the (life-)time +% observations of the \code{abrem} argument when not using rank regression. +% } +%} +\seealso{ + \code{\link{Abrem}}, + \code{\link[abremPivotals]{AbPval}} +} +\references{ + \itemize{ + \item + \emph{Improved Goodness of Fit: + P-value of the Correlation Coefficient}, Wes Fulton + + \url{http://www.barringer1.com/jun05prb.htm} + \item + \url{http://www.openreliability.org/r_and_d.html} + + Discussion on prr, pve and AbPval by Jacob Ormerod. + } +} +\examples{ +da1 <- Abrem(runif(5,100,2000)) +da1 <- abrem.fit(da1,dist="weibull",method.fit=c("rr","xony"),pch=3) + +AbPval <- da1$fit[[1]]$gof$AbPval +message("Abernethy's P-value goodness-of-fit of first fit:") +message(paste0(AbPval,ifelse(AbPval >= 10," -> Good fit."," -> BAD fit!"))) + +da1 <- abrem.fit(da1,dist="weibull",method.fit="mle",col="red1") +da1 <- abrem.fit(da1,dist="weibull",method.fit="mle-rba",col="orange3",lty=2) +da1 <- abrem.fit(da1,dist="lognormal",method.fit=c("rr","xony"),col="steelblue3",pch=8) +#da1 <- abrem.conf(da1) +me <- "Comparison between RR, MLE, MLE-RBA and lognormal2p" +plot(da1,main=me);message(me) +message(paste0( + " RR : beta=",signif(da1$fit[[1]]$beta), + ", eta=",signif(da1$fit[[1]]$eta))) +message(paste0( + " MLE: beta=",signif(da1$fit[[2]]$beta), + ", eta=",signif(da1$fit[[2]]$eta))) +message(paste0( + " MLE: beta=",signif(da1$fit[[3]]$beta), + ", eta=",signif(da1$fit[[3]]$eta))) + +### threshold parameter usage demo ### +data(abrem_mix1) +earlyda <-abrem_mix1[1:10] +midda <-abrem_mix1[11:131] +endda <-abrem_mix1[132:200] +da <-Abrem(abrem_mix1,col="gray",pch=1, + label="Complete, unaltered dataset") +da21 <-Abrem(fail=endda,susp=c(earlyda,midda),col="black",pch=19) +da22 <-Abrem(fail=endda,susp=c(earlyda,midda),col="blue",pch=3) +da23 <-Abrem(fail=endda,susp=c(earlyda,midda),col="green3",pch=4) +da21 <- abrem.fit(da21, + label="threshold=FALSE",dist="weibull3p",threshold=FALSE) +da22 <- abrem.fit(da22, + label="threshold=TRUE",dist="weibull3p",threshold=TRUE) +da23 <- abrem.fit(da23, + label="threshold=5000",dist="weibull3p",threshold=1800) +plot.abrem(list(da,da21,da22,da23),xlim=c(10,1e5), + main="Threshold parameter usage") +} diff --git a/man/contour.abrem.Rd b/man/contour.abrem.Rd new file mode 100644 index 0000000..a9da580 --- /dev/null +++ b/man/contour.abrem.Rd @@ -0,0 +1,98 @@ +\name{contour.abrem} +\alias{contour.abrem} +%- Also NEED an '\alias' for EACH other topic documented here. +\title{ + MLE Contour Plotting +} +\description{ + This function adds the \code{.abrem} method to \code{\link[graphics]{contour}} + from the \pkg{graphics} package. + + Currently, the function plots all MLE contours, found in an \code{abrem} + object or a list of \code{abrem} objects. +} +\usage{\method{contour}{abrem}(x, \dots)} +\arguments{ + \item{x}{ + Object of class \code{"abrem"} or a list of \code{abrem} objects. + } + \item{\dots}{ + Options for plotting the contours in the \code{abrem} object; see section "Details". + } +} +\details{ + The \code{\dots} argument can be any graphical parameter that can be + supplied to \code{\link[graphics]{plot.default}}, and any option that can be + set by the function \code{\link{options.abrem}}. The options set + in this way are applied to all graphical elements of the plot, overriding + any previously supplied options. + + One can pass a list of \code{abrem} objects to \code{\link{contour.abrem}}; in + that case it is mandatory to use the full method name: \code{contour.abrem(\dots)} + and not \code{contour(\dots)}. + + Currently, MLE contours are added to the abrem object by calculating + Likelihood Ratio confidence bounds for B-lives (see \code{\link{abrem.conf}}); + this function calls the \code{MLEXXXContour} function from package \pkg{debias}, where + \code{XXX} is an abbreviation for a supported distribution. Currently, only + two-parameter Weibull is supported, so MLE contours can only + be calculated for this distribution. The latter function's output is stored + in the abrem object and its output can be used to plot the contours. + + When the distribution fitting type is either MLE or MLE-RBA, + the MLE point (by definition always calculated with MLE or MLE-RBA) + of the contours coincides exactly with the calculated distribution + parameters. In all other cases, there will be a spatial difference + between the two. In that case, the correct MLE point is added to the plot as + a dot. + + In the future, calculating contours using a function like \code{abrem.contour()} + for adding them to the \code{abrem} object as well displaying them will be added + to package \pkg{abrem}. In the mean time, check out the contour + functions from package \pkg{debias} for customized contour calculations. + + As this function is still in development, no legends are currently plotted. + However, clever usage ot \code{plot.abrem} can provide all necessary + information for interpreting the contour plot (see section "Examples"). + +} +\value{ + Currently, the function returns no value. +} + +\author{Jurgen Symynck \email{jusy@openreliability.org}} + +\examples{ +## some standard options ## +defaults <- options.abrem() +options.abrem(method.conf.blives="lrb",is.plot.cb=FALSE) + +## simple example ## +da <- abrem.fit(Abrem(c(round(rweibull(5,3,1000))))) +da <- abrem.conf(da) +par(mfrow=c(1,2)) +plot(da) +contour(da) + +## multiple datasets and contours ## +## while prevention excessive plotting of CB ## +fail1 <- c(round(rweibull(5,1,1000))) +fail2 <- c(round(rweibull(8,2,3000))) + +da1 <- abrem.fit(Abrem(time=fail1,pch=0,col="red",label="First test batch")) +da2 <- abrem.fit(Abrem(time=fail2,pch=4,col="blue",label="Second test batch")) + +da1 <- abrem.conf(da1,cl=0.5 ,in.legend=FALSE) +da1 <- abrem.conf(da1,cl=0.9 ,is.plot.cb=TRUE) +da1 <- abrem.conf(da1,cl=0.95,in.legend=FALSE) +da2 <- abrem.conf(da2,cl=0.5 ,in.legend=FALSE) +da2 <- abrem.conf(da2,cl=0.9 ,is.plot.cb=TRUE) +da2 <- abrem.conf(da2,cl=0.95,in.legend=FALSE) + # prevent plotting of some confidence bounds to not overload the plot, + # also prevent inclusion in the legend. +set <- list(da1,da2) +par(mfrow=c(1,2)) +plot.abrem(set,xlim=c(5,5e5)) +contour.abrem(set) +invisible(options.abrem(defaults)) +} diff --git a/man/datasets.Rd b/man/datasets.Rd new file mode 100644 index 0000000..c1047b9 --- /dev/null +++ b/man/datasets.Rd @@ -0,0 +1,29 @@ +\name{datasets} +\alias{abrem_mix1} +\docType{data} +\title{ + Synthetic Datasets. +} +\description{ + A collection of synthetic datasets. Currently, only one dataset is provided +} + +\usage{ + abrem_mix1 +} +\format{ +% Either a numeric vector, or an ordered dataframe with columns \code{$time} +% ((life-)time observation s) and \code{$event} (censoring indicators). + + \describe{ + \item{\code{abrem_mix1}}{ + "Diaphragm life data of an acid gas compressor". A numeric vector + with complete, mixed data (200 failures) with multiple failure modes, + modeled on a real (proprietary) dataset. + } + } + + All datasets have comments; see \code{\link[base]{comment}} for more info. + +} +\keyword{datasets} \ No newline at end of file diff --git a/man/datasets_for_debugging.Rd b/man/datasets_for_debugging.Rd new file mode 100644 index 0000000..013c81d --- /dev/null +++ b/man/datasets_for_debugging.Rd @@ -0,0 +1,70 @@ +\name{debugging datasets} +\alias{abrem2} +\alias{abrem3} +\alias{abrem13} +\alias{abrem13c1} +\alias{abrem13c2} +\alias{abrem13c3} +\alias{abrem13c4} +\docType{data} +\title{ + Synthetic Datasets, for Demo's and Debugging \pkg{abrem} . +} +\description{ + A collection of synthetic datasets for debugging the \pkg{abrem} package, + created using the \code{\link{params.to.ob}} function. When these datasets + are fitted with a two-parameter Weibull distribution using X-on-Y rank + regression on exact 'median' probability plotting positions + (\code{ppos = "beta"}), the resulting parameters \code{beta} and + \code{eta} are (almost) exactly 3 and 1000, respectively. +} +\note{ + The observations in this dataset are rounded output values from + \code{params.to.ob}. Even higher precision can be achived by using the output of + \code{params.to.ob} directly, as evidenced by comparing the goodness-of-fit + indicators from both fits. +} +\usage{ + abrem2 + abrem3 + abrem13 + abrem13c1 + abrem13c2 + abrem13c3 + abrem13c4 +} +\format{ + An ordered dataframe with columns \code{$time} + ((life-)time observations) and \code{$event} (censoring indicators). + + All datasets have comments; see \code{\link[base]{comment}} for more info. + + \describe{ + \item{\code{abrem2}}{ + Two failures, no censoring. + } + \item{\code{abrem3}}{ + Three failures, no censoring. + } + \item{\code{abrem13}}{ + Thirteen failures, no censoring. + } + \item{\code{abrem13c1}}{ + Type II censored dataset with 13 observations: two failures + and 11 suspensions (at the same time as the latest failure). + } + \item{\code{abrem13c2}}{ + Type II censored dataset with 13 observations: three failures + and 10 suspensions (at the same time as the latest failure). + } + \item{\code{abrem13c3}}{ + Type II censored dataset with 13 observations: ten failures + and three suspensions (at the same time as the latest failure). + } + \item{\code{abrem13c4}}{ + Type II censored dataset with 13 observations: 12 failures + and one suspension (at the same time as the latest failure). + } + } +} +\keyword{datasets} diff --git a/man/options.abrem.Rd b/man/options.abrem.Rd new file mode 100644 index 0000000..503c475 --- /dev/null +++ b/man/options.abrem.Rd @@ -0,0 +1,431 @@ +\name{options.abrem} +\alias{options.abrem} +\alias{cex.points} +\alias{cl} +\alias{col} +\alias{col.grid} +\alias{conf.blives.sides} +\alias{conf.what} +\alias{coordinate.text.size} +\alias{dist} +\alias{in.legend} +\alias{in.legend.blives} +\alias{in.legend.gof} +\alias{is.plot.cb} +\alias{is.plot.fit} +\alias{is.plot.grid} +\alias{is.plot.legend} +\alias{is.plot.ppp} +\alias{is.plot.pppcoordinates} +\alias{label} +\alias{legend.text.size} +\alias{log} +\alias{lty} +\alias{lwd} +\alias{lwd.points} +\alias{main} +\alias{main.contour} +\alias{mar} +\alias{method.conf.blives} +\alias{method.fit} +\alias{pch} +\alias{persistent} +\alias{pivotals} +\alias{ppos} +\alias{S} +\alias{signif} +\alias{sub} +\alias{threshold} +\alias{unrel} +\alias{unrel.n} +\alias{verbosity} +\alias{xlab} +\alias{xlim} +\alias{ylab} +\alias{ylim} + +%- Also NEED an '\alias' for EACH other topic documented here. +\title{ + \pkg{abrem} Options +} +\description{ + This function handles the various calculation, printing and plotting options + of the \pkg{abrem} package. +} +\usage{options.abrem(\dots)} +\arguments{ + \item{\dots}{ + Options for calculating, printing and plotting with the \pkg{abrem} + package. + } +} + +\details{ + This function borrows its internal structure from + the \code{\link[graphics]{par}} function of package \pkg{graphics}. + It can be used in different ways: + + \code{options.abrem()} + + Returns the currently used options and there values. + + + \code{options.abrem()$dist \cr + options.abrem("dist")} + + Returns the current value of an option. + + \code{options.abrem(cl=0.95) \cr + options.abrem(list(cl=0.95,S=5e4))} + + Sets the specified options. + + Currently, there is no way to reset the options to the default values using + this function. One might, before changing any options, store the option + list in a temporary variable like + + \code{abrem.defaults <- options.abrem()} + + for restoring it later by running \code{options.abrem(abrem.defaults)}. + + The function creates a globally accessible list named \code{options_abrem}, + holding the options. One should always use the \code{options.abrem} + function to access the option list, do not access this list directly. + +} + +\section{Abrem options}{ + + \describe{ + \item{\code{cex.points}}{ + A number describing the relative size of the datapoint glyphs. + + Defaults to \code{1}. + } + \item{\code{cl}}{ + Confidence level: A single number from the interval \code{[0,[1} + specifying the confidence level for various confidence calculations. + + Defaults to \code{0.9}. + } + \item{\code{col}}{ + An integer or character string describing the color of a + graphical entity. + See function \code{\link[grDevices]{colors}} for available colors. + + Defaults to \code{"black"}. + } + \item{\code{col.grid}}{ + An integer or character string describing the color of the grid + + Defaults to \code{"gray"}. + } + \item{\code{conf.blives.sides}}{ + Either \code{"lower"}, \code{"upper"} or \code{"double"}, + specifying the type of bound(s) to be calculated. + + Defaults to \code{c("double")}. + } + \item{\code{conf.what}}{ + A vector of class \code{"character"} describing for which entities + that confidence should be calculated. + + Defaults to \code{c("blives")}. + } + \item{\code{coordinate.text.size}}{ + A number determining the relative coordinate text label size. + + Defaults to \code{0.7}. + } + \item{\code{dist}}{ + A character string with the target distribution for fitting. + See \code{\link{abrem.fit}} for in-depth discussion of the settings. + + Defaults to \code{"weibull"}. + } + \item{\code{in.legend}}{ + Logical value controlling the inclusion of various elements in + the legend. See \code{\link{abrem.fit}} and \code{\link{abrem.conf}} + for more details. + + Defaults to \code{TRUE}. + } + \item{\code{in.legend.blives}}{ + Logical value controlling the inclusion of B-life confidence bounds + in the legend. + + Defaults to \code{TRUE}. + } + \item{\code{in.legend.gof}}{ + Logical value controlling the inclusion of goodness-of-fit + indicators in the legend. Note that this does \emph{not} prohibit + calculation of sometimes time consuming \code{prr} values, it only + prevents inclusion in the legend. + %% TODO: check if this is still valid (prr <-> pve <-> AbPval) + + Defaults to \code{TRUE}. + } + \item{\code{is.plot.cb}}{ + Logical value controlling the plotting of various types of + confidence bounds (if any present in the \code{abrem} object). + + Defaults to \code{TRUE}. + } + \item{\code{is.plot.fit}}{ + Logical value controlling the plotting of the fitted line. + + Defaults to \code{TRUE}. + } + \item{\code{is.plot.grid}}{ + Logical value controlling the plotting of the grid. + + Defaults to \code{TRUE}. + } + \item{\code{is.plot.legend}}{ + Logical value controlling the plotting of the legend. + + Defaults to \code{TRUE}. + } + \item{\code{is.plot.ppp}}{ + Logical value controlling the plotting of the + (life-)time observations' probability plot positions. + + Defaults to \code{TRUE}. + } +%% \item{\code{is.plot.pppvalues}}{ +%% Logical value controlling the plotting of the text labels +%% (the coordinates on the plot) next to the plot positions. +%% +%% Defaults to \code{FALSE}. +%% } + \item{\code{label}}{ + A character string with the title of the legend box, can be used to + label a dataset or fit. + + Defaults to \code{NULL}. + } + \item{\code{legend.text.size}}{ + A number determining the relative legend text size. + + Defaults to \code{0.7}. + } + \item{\code{log}}{ + A character string describing the type of plotting canvas to use. + Possible values are \code{"x"} (weibull canvas), \code{"xy"} + (lognormal canvas), \code{""} (weibull Y axis, linear X axis) + and \code{"y"} (lognormal Y axis, linear X axis). + + + Defaults to \code{"x"}, Weibull canvas. + For contour plots, \code{log} is always \code{""} and cannot be set. + } + \item{\code{lty}}{ + An integer describing the line type of both the fitted + line and the confidence bounds, if any. + See function \code{\link[graphics]{par}} for available line types. + + Defaults to \code{1}. + } + \item{\code{lwd}}{ + An integer describing the line width of both the fitted + line and the confidence bounds, if any. + + Defaults to \code{2}. + } + \item{\code{lwd.points}}{ + An integer describing the thickness of the datapoint glyphs. + + Defaults to \code{2}. + } + \item{\code{main}, \code{main.contour}}{ + A character string; the main title of the plot or contour plot. + + Defaults to \code{"Probability Plot"} or \code{"Contour Plot"} . + } + \item{\code{mar}}{ + A numerical vector of the form \code{c(bottom, left, top, right)} + which gives the number of lines of margin to be specified on + the four sides of the plot. + See function \code{\link[graphics]{par}} for more info. + + Defaults to \code{c(5.1,4.1,5.1,2.1)}. + } + \item{\code{method.conf.blives}}{ + A vector of class \code{"character"} describing the techniques used + for calculating confidence for B-lives. + See \code{\link{abrem.conf}} for in-depth discussion of the settings. + + Defaults to \code{c("mcpivotals")}. + } + \item{\code{method.fit}}{ + A vector of class \code{"character"} with fitting options. + See \code{\link{abrem.fit}} for in-depth discussion of the settings. + + Defaults to \code{c("rr","xony")}. + } + \item{\code{pch}}{ + An integer or single character describing the plotting symbol, used for + plotting the datapoints. For more info, see + \code{\link[graphics]{points}}. + + Defaults to \code{1}. + } + \item{\code{persistent}}{ + Experimental parameter, do not use. + + Defaults to \code{TRUE}. + } + \item{\code{pivotals}}{ + Currently unused. Defaults to \code{FALSE}. + } + \item{\code{ppos}}{ + Short for "probability plotting positions", it is a vector of + class \code{"character"} describing the type of plotting positions to be + calculated. + See \code{\link{Abrem}} for in-depth discussion of the settings. + + Defaults to \code{c("benard")}. + } + \item{\code{S}}{ + An integer describing the number of Monte Carlo simulations on + which the Monte Carlo pivotal confidence bounds and calculation + of the "prr" goodness-of-fit indicator are based (See + \code{\link[abremPivotals]{pivotalMC}} of package \pkg{abremPivotals} for more + details). + % TODO: rewrite, since there is no mention anymore of prr in abremPivotals + % Actually, there is access to prr somewhere... + + Defaults to \code{10000}. + } + \item{\code{signif}}{ + An integer describing the significant digits of various numbers + that are displayed in the legend. + + Defaults to \code{4}. + } + \item{\code{sub}, \code{sub.contour}}{ + A character string; the subtitle of the plot or contour plot. + + Defaults to \code{""}. + } + \item{\code{threshold}}{ + A logical value specifying whether a threshold time value should be + substracted from the graphical elements before plotting. This is + particulary useful when fitting three-parameter models like + \code{weibull3p}, where the third model parameter is a threshold value. + + If \code{TRUE}, the software will subtract the threshold value + (that is to be calculated by \code{abrem.fit}) from its + associated graphical elements like plot poistions, fits and confidence bounds. + + If a numeric value is passed, the software will subtract this value + from the appropriate graphical elements. + + Note: currently, \code{"threshold"} should only be used in combination with + \code{Abrem} and \code{abrem.fit} and not with + \code{plot.abrem}. See \code{abrem.fit} for more info. + + + Defaults to \code{FALSE} + } + \item{\code{unrel}}{ + An unordered numeric vector with unreliability levels for which + B-life confidence will be calculated. + + Defaults to \code{c(0.1,0.05,0.01)}. + } + \item{\code{unrel.n}}{ + An integer controlling the amount of unreliability levels for + which confidence bounds are calculated and ultimately plotted. + See \code{\link{abrem.conf}} for more details. + + Defaults to \code{25} + } + \item{\code{verbosity}}{ + An integer from \code{c(1,2,3)} setting the level of verbosity of + all \pkg{abrem} functions. Use this for debugging. + + Defaults to \code{1}. + } + \item{\code{xlab}, \code{ylab}}{ + Character strings with the labels for the X and Y axis. + + Default to \code{"Time To Failure"} and + \code{"Unreliability [\%]"} respectively. + For contour plots, Default to \code{"Eta"} and \code{"Beta"} respectively. + } + \item{\code{xlim}}{ + A numeric vector with two values determining the plotting range + of the horizontal axis of the plot. + \code{plot.abrem} calculates horizontal limits automatically + from its data argument. + + Defaults to \code{NULL}. + } + \item{\code{ylim}}{ + A numeric vector with two values determining the plotting range + of the vertical axis of the plot. Allowed values come + from the interval \code{]0,1[}. \code{plot.abrem} calculates + vertical limits automatically from its data argument; if vertical + plotting positions do not exceed the standard range of + \code{c(0.01,0.99)}, then the latter is used. + Otherwise, the standard range is expanded to include the + extreme plot positions. + + The standard vertical range is currently hardcoded to \code{c(0.01,0.99)}. + + Defaults to \code{NULL}. + } + } +} + +\value{ + Executing \code{options.abrem} without arguments returns a named list + containing the currently active global options of the \pkg{abrem} package. + + When arguments were supplied, these are returned in a named list. +} +\author{Jurgen Symynck \email{jusy@openreliability.org}} +\note{ + Typical usage of \pkg{abrem} involves calling a sequence of functions like: + + \code{da <- Abrem(c(10,11,27))}\cr + \code{da <- abrem.fit(da)}\cr + \code{da <- abrem.conf(da)}\cr + \code{plot(da)} + + Do not call \code{options.abrem} in between these functions because some + options are locked and cannot be altered further in this chain. This is an + implication of the way the \code{abrem} object is structured. + + The correct time to specify an option is when it is needed for the first + time. For example, when the color setting option \code{col = "red"} is + passed as an argument of function \code{Abrem}, it will be used for + datapoints, fits and confidence bounds. If supplied to \code{\link{abrem.conf}}, + only the confidence bounds will have the specified color, hereby overriding + any previously inherited color settings from \code{\link{abrem.fit}} + or \code{\link{Abrem}}. +} +%\section{To Do}{ +% \itemize{ +% \item Add the possibility to reset to the default values. +% \item Implement returning of the pivotals dataframe in the +% \code{abrem} object. +% \item Add detection of unsupported arguments. +% } +%} +\examples{ +## backup options ## +abrem.defaults <- options.abrem() + +## setting new options ## +options.abrem(S=5e5,cl=0.99) +%%options.abrem(sub="Testing options.abrem()") + +## listing options ## +options.abrem() +options.abrem()$main + +## restore options ## +options.abrem(abrem.defaults) +} diff --git a/man/params.to.ob.Rd b/man/params.to.ob.Rd new file mode 100644 index 0000000..8e1d675 --- /dev/null +++ b/man/params.to.ob.Rd @@ -0,0 +1,85 @@ +\name{params.to.ob} +\alias{params.to.ob} +%- Also NEED an '\alias' for EACH other topic documented here. +\title{ + Create (un-)Censored Testdata from Distribution Parameters} +\description{ + Create a set of (life-)time observations (possibly with censoring), + perfectly matching a given distribution. +} +\usage{params.to.ob(dist, \ldots)} +%- maybe also 'usage' for other objects documented here. +\arguments{ + \item{dist}{ + The target distribution for creating the (life-)time observations. + } + \item{\ldots}{ + Named arguments for calculating the dataset such as slope, shape, + number and event vector. See section "Details". + } +} +\details{ + This function can be used for testing purposes. Internally, it is used for + the experimental calculation of Monte Carlo Pivotal confidence bounds for + right censored (life-)time observations. + + \describe{ + \item{generating Weibull datasets}{ + \code{params.to.ob("weibull",beta=3,eta=1000,n=5)} + } + \item{generating lognormal datasets}{ + \code{params.to.ob("lognormal",mulog=log(1000),sigmalog=log(2),n=5)} + } + \item{censoring, event vector}{ + \code{params.to.ob("weibull",beta=3,eta=1000,event=c(1,1,1,0,0))} + } + } + When \code{\link{abrem.fit}} is called on an \code{abrem} object based on + these (life-)time observations, the same fit parameters will be found as + those used to generate the dataset. +} +\value{ + A dataframe with two columns: + \describe{ + \item{\code{$time}}{ + An ordered vector with (life-)time observations. + } + \item{\code{$event}}{ + A vector of class \code{"numeric"} with right-censoring indicators. + See \code{\link{Abrem}} for more details on the indicators. + } + +%% If argument \code{x} was a dataframe with additional columns, these will also be returned. + } +} +\author{Jurgen Symynck \email{jusy@openreliability.org}} +\note{ + Currently, only distributions fitted with + \code{method.fit = c("rr","xony")} + are supported. +} +\examples{ +## generate three synthetic datasets ## +d1 <- params.to.ob("weibull",beta=3,eta=1000,n=10) +d1 <- abrem.fit(Abrem(d1),lwd=1) +print(d1$data) +message(paste0( + " beta=",d1$fit[[1]]$beta, + ", eta=",d1$fit[[1]]$eta)) + +d2 <- params.to.ob("weibull",beta=3,eta=1000, + event=c(1,1,0,0,0,1,0,1,0,0)) +d2 <- abrem.fit(Abrem(d2,pch=3,col="red",cex.points=1.5),lty=3,lwd=3) +print(d2$data) +message(paste0( + " beta=",d2$fit[[1]]$beta, + ", eta=",d2$fit[[1]]$eta)) + +d3 <- params.to.ob("lognormal",mulog=log(1000),sigmalog=log(2),n=10) +d3 <- abrem.fit(Abrem(d3,pch=0),dist="lognormal") +print(d3$data) +message(paste0( + " mulog=",d3$fit[[1]]$mulog, + ", sigmalog=",d3$fit[[1]]$sigmalog)) +plot.abrem(list(d1,d2,d3),main="Demo of params.to.ob()") +} \ No newline at end of file diff --git a/man/plot.abrem.Rd b/man/plot.abrem.Rd new file mode 100644 index 0000000..aac3579 --- /dev/null +++ b/man/plot.abrem.Rd @@ -0,0 +1,87 @@ +\name{plot.abrem} +\alias{plot.abrem} +%- Also NEED an '\alias' for EACH other topic documented here. +\title{ + \code{Abrem} Object Plotting +} +\description{ + This function adds the \code{.abrem} method to \code{\link[graphics]{plot}} + from the \pkg{graphics} package. + + Currently, the function plots the (life-)time observations, fits (if any) + and confidence bounds for B-lives (if any) of an + \code{abrem} object or a list of \code{abrem} objects on Weibull or + lognormal probability paper. + + For each fit in the (list of) \code{abrem} object(s), legends are added to + the plot, displaying the fit parameters and (if available) goodness-of-fit + indicators and confidence information. +} +\usage{\method{plot}{abrem}(x, \dots)} +\arguments{ + \item{x}{ + Object of class \code{"abrem"} or a list of \code{abrem} objects. + } + \item{\dots}{ + Options for plotting the \code{abrem} object; see section "Details". + } +} +\details{ + The \code{\dots} argument can be any graphical parameter that can be + supplied to \code{\link[graphics]{plot.default}}, and any option that can be + set by the function \code{\link{options.abrem}}. The options set + in this way are applied to all graphical elements of the plot, overriding + any previously supplied options. + + One can pass a list of \code{abrem} objects to \code{\link{plot.abrem}}; in + that case it is mandatory to use the full method name: \code{plot.abrem(\dots)} + and not \code{plot(\dots)}. + + The calculated Weibull or lognormal distribution fits are plotted + on Weibull probability paper by default, but by passing the argument + \code{\link{log} = "xy"} to the function, lognormal paper is used + (see \code{\link{options.abrem}}). + + When a \emph{list} of abrem objects is passed, the plot window is generated + with the options of the first \code{abrem} object in the list. + % TODO: list these options + +% A thin dashed line appears close to the fitted line (if any) when confidence +% are present and plotted. +} +\value{ + Currently, the function returns no value. +} +%\section{To Do/Known issues/Bugs}{ +% \itemize{ +% \item The legend positions are slightly to the right, overlapping the +% right Y axis (especially for small horizontal ranges). +% \item Have the function return the updated \code{abrem} arguments. +% } +%} +\author{Jurgen Symynck \email{jusy@openreliability.org}} + +\seealso{ + \code{\link{options.abrem}} +} +\examples{ +da1 <- Abrem(rweibull(5,3,1000)) +da1 <- abrem.fit(da1,dist="weibull", method.fit=c("rr","xony"),pch=3) +da1 <- abrem.fit(da1,dist="weibull", method.fit="mle",col="red1") +da1 <- abrem.fit(da1,dist="weibull", method.fit="mle-rba",col="orange",lty=2) +da1 <- abrem.fit(da1,dist="lognormal",method.fit=c("rr","xony"), + col="steelblue3",pch=8) +da1 <- abrem.conf(da1) +plot(da1,col="black") # plot all in black +plot(da1,main="Weibull distributed observations") # plot with specified colors + +### usage of lists ### +options.abrem(blives=0.1) # make the legend boxes a bit shorter... +da2 <- abrem.conf(abrem.fit(Abrem(runif(5,10,100),col="red"))) +da3 <- abrem.conf(abrem.fit(Abrem(rweibull(5,2,1000),col="green4",pch=3))) +da4 <- abrem.conf(abrem.fit(Abrem(rlnorm(5,log(500),log(2)),col="blue3",pch=8), + dist="lognormal")) + +plot.abrem(list(da2,da3,da4),xlim=c(1,1e6), + main="Uniformly distributed observations") +} diff --git a/man/print.abrem.Rd b/man/print.abrem.Rd new file mode 100644 index 0000000..7791fee --- /dev/null +++ b/man/print.abrem.Rd @@ -0,0 +1,26 @@ +\name{print.abrem} +\alias{print.abrem} +%- Also NEED an '\alias' for EACH other topic documented here. +\title{ + \code{Abrem} Object Printing +} +\description{ + This function prints the \code{abrem} object while - for clarity - omitting + the long \code{$options} named list. +} +\usage{\method{print}{abrem}(x, \dots)} +\arguments{ + \item{x}{ + Object of class \code{"abrem"}. + } + \item{\dots}{ + Further arguments passed to or from other methods. + } +} + +\author{Jurgen Symynck \email{jusy@openreliability.org}} + +\examples{ +da <- abrem.fit(Abrem(rweibull(13,3,1000))) +print(da) +} diff --git a/vignettes/using_abrem.Rnw b/vignettes/using_abrem.Rnw new file mode 100644 index 0000000..0854c0d --- /dev/null +++ b/vignettes/using_abrem.Rnw @@ -0,0 +1,495 @@ +%% LyX 2.0.7 created this file. For more info, see http://www.lyx.org/. +%% Do not edit unless you really know what you are doing. +\documentclass[american,english,noae]{article} +\usepackage{lmodern} +\renewcommand{\sfdefault}{lmss} +\renewcommand{\ttdefault}{lmtt} +\usepackage[T1]{fontenc} +\usepackage[latin9]{inputenc} + +\makeatletter +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Textclass specific LaTeX commands. +<>= + if(exists(".orig.enc")) options(encoding = .orig.enc) +@ + +\makeatother + +\usepackage{babel} +\begin{document} + + + + +<>= +invisible(library(abrem)) +invisible(options.abrem(sub=paste0("(plot created with abrem ",packageVersion("abrem"),")"))) +#invisible(options.abrem(mar=c(4,4,4,1))) +@ + + + + +\title{Using abrem} + + +\author{David Silkworth, Jurgen Symynck} +\maketitle +\begin{abstract} +Functions of the abrem package provide an application layer for entering +life data, selecting various options for the treatment of that data, +generation of a graphical display, and presenting meaningful statistical +output in a legend. This application layer is executed in the R console, +based on scripts that are likely prepared in a separate editor. + +Package abrem makes use of the R object model to prepare objects for +plotted output. Various options are handled to extend the object with +parametric fitting, preparation of confidence interval bounds, and +preparation of the graphical output. +\end{abstract} + +\section{Using of abrem for default two--parameter Weibull plotting} + +For a rapid Weibull two-parameter fit and display, an object is constructed +using the \texttt{Abrem()} function on the data of interest. Then +a fit element is added using function \texttt{abrem.fit()}. This object +can then be plotted with the \texttt{.abrem} method of the standard +R function \texttt{plot()}. To demonstrate this simple use, the initial +problem sets from Chapter 2 of \textquotedblleft{}The New Weibull +Handbook, Fifth edition\textquotedblright{} \cite{abernethy2008new}are +shown below. For each example just copy the code script and paste +it into the R console to observe the graphical output in full size. + +\begin{center} +<>= +library(abrem) +Prob2.1Data <- c(150,85,250,240,135,200,190) + # data from Problem 2.1 from the Handbook +Prob2.1Object <- Abrem(Prob2.1Data) +Prob2.1Object <- abrem.fit(Prob2.1Object) +plot(Prob2.1Object, main="Problem 2-1") +@ +\par\end{center} + +This is a first, simple, plotting example. In the script above the +first line merely assures that the abrem package and all its dependencies +have been loaded in this R session. + +The second line places the example data into a numeric vector. + +The third line uses function \texttt{Abrem()} to construct an abrem +object from the data. The Abrem function has been designed for some +flexibility of data input. In this case, a single numeric vector has +been passed as the first argument, so it is taken to be a set of complete +failure times. This initial object could have been plotted at this +point and one would observe the points, only, plotted on the default +canvas for a Weibull analysis. By default the fitting distribution +is \texttt{``weibull2p\textquotedblright{}}, and the plotting positions +have been determined using the beta binomial function for most technically +accurate median rank positions. These positions are minimally different +than the Benard approximation method used in the \textquotedblleft{}Handbook\textquotedblright{} +text. + +The fourth line of the script adds a fit of the data into the abrem +object by applying the\texttt{ abrem.fit()} function on the object, +extending it. By default this fit is by rank regression using the +preferred X--on--Y regression convention. + +Finally, the fifth line plots the object based on the exposed elements +of the abrem object, which has been provided as a first argument to +the \texttt{plot()} function. In this case the main title of the chart +has been specified. Upon execution of this line an R Graphics Device +window is opened with the resulting display. + + +\subsection*{Now consider example Problem 2--2 from the text} + +\begin{center} +<>= +library(abrem) +serials <- c(831, 832, 833, 834, 835, 836, 837, 838) +times <- c( 9, 6,14.6, 1.1, 20, 7, 65, 8) +events <- c( 1, 1, 0, 1, 1, 0, 1, 0) +Prob2.2Data <- data.frame(serial=serials,time=times,event=events) +Prob2.2Object <- Abrem(Prob2.2Data, col="red") +Prob2.2Object <- abrem.fit(Prob2.2Object) +plot(Prob2.2Object, main="Problem 2-2") +@ +\par\end{center} + +In the script above the data is entered into three vectors in the +exact order as presented in a table in the text. These vectors are +used to create a dataframe to present the data with specifically named +fields of \texttt{``time''} and \texttt{``event''}. Note that +the event vector is made up of 1's for failures and 0's for suspensions. +The \texttt{Abrem()} function will accept this format as a first argument. +It will parse the dataframe for the two specifically named fields +as shown singular and all lower case. Additional information that +may be stored in the dataframe is not used on the abrem object that +will be created by function \texttt{Abrem()}. In this example some +color has been added to the object. A default rank regression fit +is then added to the object as in the previous example. Finally the +plot is executed providing a new main title for the chart. + + +\subsection*{Placing more than one fitted object on one chart.} + +\begin{center} +<>= +library(abrem) +Prob2.1Data <- c(150,85,250,240,135,200,190) +Prob2.1Object <- Abrem(Prob2.1Data, col="blue", + label="Problem 2-1") +Prob2.1Object <- abrem.fit(Prob2.1Object) +P2.2Failures <- c(9.0,6.0,1.1,20.0,65.0) +P2.2Suspensions <- c(14.6,7.0,8.0) +Prob2.2Object <- Abrem(fail=P2.2Failures, + susp=P2.2Suspensions, col="red", label="Problem 2-2") +Prob2.2Object <- abrem.fit(Prob2.2Object) +plot.abrem(list(Prob2.2Object, Prob2.1Object), + main="Chapter 2 Problems") +@ +\par\end{center} + +In this script some color is added to each problem object and identifying +labels are added for the legend. The data for the second problem has +been prepared for an alternate input means for object initialization +by the \texttt{Abrem()} function. By this convention two vectors of +time values are prepared; one for complete failures, the other for +suspension data. These two vectors must then be assigned to the specifically +named arguments \texttt{``fail''} and \texttt{``susp''} in the +\texttt{Abrem()} function call. + +The procedure for placing two or more abrem objects into a chart is +to explicitly pass the objects in a list to the \texttt{.abrem} method +of the generic \texttt{plot()} fucntion: this is a special case for +plotting abrem objects since the generic R \texttt{plot()} function +cannot recognize such a list of abrem objects. + + +\subsection*{Adding confidence intervals for B--lives to abrem objects } + +Continuing a convention of recalculating examples from our reference +text, two examples are provided below plotting recalculations for +Figure C-1 and Figure 7-5 (from Section 7.5.1) from the \textquotedblleft{}New +Weibull Handbook, Fifth Edition\textquotedblright{} below: + +\begin{center} +<>= +library(abrem) +## Figure C-1 +AppBCData <- data.frame( + time =c(1500, 1750, 2250, 4000, 4300, 5000, 7000), + event=c( 1, 0, 1, 1, 1, 0, 1)) +FigC.1 <- Abrem(AppBCData, pch=2) +FigC.1 <- abrem.fit(FigC.1, col="blue") +FigC.1 <- abrem.conf(FigC.1, col="darkgreen") +FigC.1 <- abrem.conf(FigC.1, method.conf.blives="bbb", + lty=2, lwd=1, col="black") +plot(FigC.1, xlim=c(1,5e5), ylim=c(.01,.90), + main="Figure C-1 with added Beta Binomial Bounds") +@ +\par\end{center} + +The Figure C-1 script demonstrates graphic control over the individual +elements of a single abrem object. During the object construction +using \texttt{Abrem()} the point character, \texttt{pch}, is set to +a value of \texttt{2} to generate upright triangles. By entering \texttt{?points} +in the R console a help page will appear in which all of the \texttt{pch} +options are defined. As the default fit element is created it is assigned +the color blue. As the default MCpivotal confidence bounds are created +they are assigned a dark green color. Then, as a second confidence +interval element is added as the beta binomial bounds, the line type, +\texttt{lty}, and line width, lwd, are customized according to options +that can be explored by entering \texttt{?par} in the R console. So, +it is demonstrated that each abrem object can have multiple elements +and each of those elements can be graphically controlled. + +A final observation can be noted in that the ranges for X and Y axes +for Figure C-1 were altered to provide a more pleasing presentation +than default code would have given. + + +\subsubsection*{Comparing datasets} + +\begin{center} + +\par\end{center} + +\begin{center} +<>= +library(abrem) +## Figure 7-5 (Section 7.5.1) +#set.seed(1234) +Obj1 <- params.to.ob(dist="weibull",beta=2.09,eta=985.8,n=8,ppos="beta") +Obj2 <- params.to.ob(dist="weibull",beta=3.666,eta=1984,n=8,ppos="beta") +Obj1 <- abrem.conf(abrem.fit(Abrem(Obj1, pch=2, + col="orange")), unrel=0.1) +Obj2 <- abrem.conf(abrem.fit(Abrem(Obj2, pch=6, + col="forestgreen")), unrel=0.1) +plot.abrem(list(Obj1,Obj2), + main="Figure 7-5: Compare Data Sets at B10 with 95% Confidence\n") +@ +\par\end{center} + + + +The script for Figure 7-5 begins with a seed setting establising a +uniform state for the random number generator (RNG) as the function +is called. This setting assures consistent numeric output of the Monte +Carlo output regardless of any previous RNG acctivity. The next two +lines establish exact observation times that will result in the given +Weibull parameter pairs after fitting, using abrem's \texttt{params.to.ob()} +function.The latter is a wrapper for R's Weibull quantile function +\texttt{qweibull()} that returns the time value for a given Weibull +parameter pair at a given unreliability percentage. + +A nested approach was used to expand each object with all desired +elements and their graphical features. Notice that the B-lives output +was limited to a single report at B10 for this case. The \texttt{unrel} +option could have been set at any point in the construction of the +objects in this case; association with the \texttt{abrem.conf()} just +seemed most appropriate. Finally, as seen before, these objects were +passed in a list to the object method \texttt{plot.abrem()} with a +definition for the main chart title. + +The point of this figure is to demonstrate a method of determining +statistical distinctness of two data sets similar to a method presented +in 1994 by Jim Lempke of Ford Motor Company. It is hard to see the +distinction that exists at the B10 level (10\% unreliability) between +the upper bound of the orange set and the lower bound of the green +set by the graph itself. It is hard to see this in the book also. +However, the B--lives display in the legend clearly shows the 95\% +upper-bound B10 value to be \Sexpr{signif(subset(Obj1$fit[[1]]$conf$blives[[1]]$bounds,signif(unrel)==0.1)$upper,4)}\texttt{ +}for the orange set, while the 95\% lower-bound B10 value is \Sexpr{signif(subset(Obj2$fit[[1]]$conf$blives[[1]]$bounds,signif(unrel)==0.1)$lower,4)} +\texttt{}for the green set. Since the confidence limit, \texttt{cl},has +been defined as a 90\% double--sided interval, the single--sided confidence +levels are each 95\%. + +At this time preparation of pivotal confidence interval bounds for +data sets with suspensions is a subject of some on--going study. The +user is cautioned that bounds and B-Life determinations prepared for +such abrem objects may be different than other software presentations. +An open invitation is made here for participation in this study (according +to \cite[p. 218]{lawless_2003}, the applied techniques should be +sufficient for arbitrary censoring schemes when the sample size is +'moderately large'). + + +\subsection*{Models beyond two--parameter Weibull} + +One of the few lognormal data examples from the text is provided in +Figure 3-13.With the appropriate data entered into R as a vector +named \texttt{F3.13da} the following script can be run to generate +the lognomal fit on the log-log canvas for depicting linear regression. + +\begin{center} +<>= +library(abrem) +F3.13da <- c(3.46623, 3.732711, 4.052996, 4.628703, 4.8157, 5.84517, + 5.888313,5.892967, 8.168362, 10.02799, 10.06062, 10.49785, + 11.11493, 11.87369, 12.21122, 12.51854, 12.91357, 18.04246, + 18.20712,19.57305, 21.20873, 30.03917, 34.88001, 36.87355, + 53.91168) +F3.13ln2 <- abrem.fit(Abrem(F3.13da),dist="lognormal",col="magenta") +plot(F3.13ln2,log="xy",main="Lognormal Plot") +@ +\par\end{center} + +Continuing the active R session with vector \texttt{F3.13da} and abrem +object \texttt{F3.13ln2} in memory, the lognormal fit we just made +can be re-plotted on the Weibull canvas along with two-parameter and +three-parameter Weibull fitted objects constructed from the same data. + +\begin{center} +<>= +F3.13w2 <- abrem.fit(Abrem(F3.13da), col="blue") +F3.13w3 <- abrem.fit(Abrem(F3.13da), dist="weibull3p", col="red") +plot.abrem(list(F3.13w2,F3.13ln2,F3.13w3), + main=("Multi-distribution Plot")) +@ +\par\end{center} + +This is a somewhat unexpected display. The example was specifically +provided as a lognormal representation, but we observe a much more +convincing curve fit with the three-parameter Weibull. A little examination +with SuperSMITH reveals that use of the PVE for fit selection will +favor the lognormal (due to a somewhat penalized 3p Weibull MC pivotal +determination), but if fitting is performed using MLE, the likelihood +ratio test reveals a 90\% significance for the 3p Weibull over the +lognormal. The \texttt{abrem} package is neutral on this topic as +of yet. For now, there is no \texttt{CCC\textasciicircum{}2} nor \texttt{prr} +value displayed for the three-parmameter model. In the future provision +of a log-likelihood value is expected as a goodness--of--fit measure +even for models fit by rank regression. + + +\subsection*{Exploratory work using third parameter optimization } + +At times sets of component life data can be too large for simple navigation +in a spreadsheet. + +\begin{center} +<>= +library(abrem) +data(abrem_mix1) +da <- abrem_mix1 +dafit <- abrem.fit(Abrem(da),col="red") +plot(dafit, main="Bathtub Life Data") +abline(v=107, col="orange", lty=5, lwd=2) +abline(v=1750, col="orange", lty=5, lwd=2) +@ +\par\end{center} + +This complete set of life data ('Diaphragm life data of an acid gas +compressor') fits poorly as a Weibull. Not only is the \texttt{AbPval} +less than \texttt{10}, it is even zero!We observe, however that the +data takes on somewhat expected behavior at the extremes of life: +an early period of \emph{infant mortality} appears to give way to +a more uniform life profile. Finally, at some later point it appears +that wear--out is ultimately overwhelming. + +The orange vertical markers have been placed simply by eye, with some +added help of examination of the data in the vicinity of the markers. +In accordance with this observation, the data set is split into three +parts for examination as three-parameter Weibulls. Each segment of +life data, so divided, is treated as representing a \emph{mode} of +failure; the other segments of data are considered as suspension data +for that particular mode. + +\begin{center} +<>= +earlyda <- da[1:10] +midda <- da[11:131] +endda <- da[132:200] +Abrem_early <- Abrem(fail=earlyda,susp=c(midda,endda), + col="orange",label="Early data") +Abrem_mid <- Abrem(fail=midda, susp=c(earlyda,endda), + col="magenta",label="Middata") +Abrem_end <- Abrem(fail=endda, susp=c(earlyda,midda), + col="blue",label="Enddata") +earlyfit <- abrem.fit(Abrem_early,dist="weibull3p") +midfit <- abrem.fit(Abrem_mid, dist="weibull3p") +endfit <- abrem.fit(Abrem_end, dist="weibull3p") +plot.abrem(list(earlyfit,midfit,endfit), legend.text.size=0.5, + xlim=c(0.5,2e5), + main="Division of Life Data Using Three-Parameter Weibull") +@ +\par\end{center} + +The resulting plot demonstrates much closer fitting to the segmented +data. + +It is possible to have \texttt{plot.abrem()} subtract a threshold +value from the plot positions and the subsequent fits by passing the +option \texttt{threshold=TRUE} in the call to \texttt{abrem.fit()}. +This will subtract the \texttt{t0} value calculated by \texttt{abrem.fit(..., +dist=''weibull3p'')} from the data and fits before plotting. + +Just add \texttt{threshold=TRUE} to the \texttt{abrem.fit()} lines +in the previous code block: + +\begin{center} +<>= +earlyfit <- abrem.fit(Abrem_early,dist="weibull3p",threshold=TRUE) +midfit <- abrem.fit(Abrem_mid, dist="weibull3p",threshold=TRUE) +endfit <- abrem.fit(Abrem_end, dist="weibull3p",threshold=TRUE) +plot.abrem(list(earlyfit,midfit,endfit), + legend.text.size=0.5,xlim=c(5,2e4), + main="Linear Three-Parameter Weibull Fits by t0 Data Adjustment") +@ +\par\end{center} + + + + +\subsection*{Maximum Likelihood Estimation (MLE) fitting with \texttt{abrem}} + +At this time MLE fitting has only rudimentarily been implemented in +the \texttt{abrem} package. It is possible to assign \texttt{method.fit} +options of either \texttt{\textquotedblleft{}mle\textquotedblright{}} +or \texttt{\textquotedblleft{}mle-rba\textquotedblright{}} to a two-parameter +Weibull abrem object, only. A simple example is drawn from Appendix +C of the text using the very small example data set found there. + +\begin{center} +<>= +library(abrem) +mle_ex <- Abrem(fail=c(1500,2250,4000,4300,7000),susp=c(1750,5000)) +mle_ex <- abrem.fit(mle_ex, method.fit=c("rr","xony")) +mle_ex <- abrem.fit(mle_ex, method.fit="mle", col="blue") +mle_ex <- abrem.fit(mle_ex, method.fit="mle-rba", col="steelblue") +plot(mle_ex, xlim=c(1000,20000), main="RR, MLE and MLE-RBA Fitting") +@ +\par\end{center} + + +\subsection*{Model Selection and Goodness--of--Fit} + +One of the criteria in deciding on which model to go for is the goodness--of--fit. +The handbook discusses many methods like $r^{2}-CCC^{2}$, Pve (in +Abrem this is replaced by the superior AbPval) and Likelihood Ratio +test. + + +\subsubsection*{Likelihood Ratio Test} + +When Maximum Likelihood Estimation is used for fitting the model parameters, +the log--likelihood is calculated and stored in the abrem object. +Alternatively---for Rank Regression fits---on can use functions \texttt{LLw()} +and \texttt{LLln()} from the \texttt{debias} package to calculate +a log--likelihood value. + +When used to establish acceptance of a three--parameter model over +its two-parameter counterpart on the same data, the three--parameter +model is considered the \emph{null--hypothesis} model (H0) while the +two--parameter model is taken as the \emph{alternate hypothesis} model +(H1. Degrees of freedom (df) is \texttt{1} for this constrained model +test. + + + +Acceptance of the 3p model should require an Likelihood Ratio Test +P-value greater than $50[\%]$. Comparison between Weibull and lognormal +for failure data would place Weibull as alternate and lognormal as +null whilDegrees of freedom would be 2. + +\begin{center} +<>= +da <- rweibull(25,2,1000) + 1000 +da <- Abrem(da) +da <- abrem.fit(da,dist="weibull",method.fit="mle-rba",col="red") + # alternative model +da <- abrem.fit(da,dist="weibull3p",method.fit="mle-rba") + # null model +LLalt <- da$fit[[1]]$gof$loglik +LLnull <- da$fit[[2]]$gof$loglik + +LRT_P <- function(LLalt,LLnull,df){ + pchisq(-2*(LLalt-LLnull),df) + # original +} +lrtp <- LRT_P(LLalt,LLnull,df=1) +plot.abrem(da) +@ +\par\end{center} + +The Likelihood Ratio Test percentile value is here \Sexpr{signif(lrtp*100,4)} +$[\%]$. If the percentile is larger than $50[\%]$, the three--parameter +null--hypothesis model is not rejected and thus preferred. If smaller, +the alternative two--parameter model is preferred. + + +\subsection*{MLE Contours} + +Further development with MLE contours and likelihood ratio bounds +is available in the underlying debias package. This is one of the +frontiers of the Abernethy Reliability Methods project. + +\selectlanguage{american}% +\clearpage{}\bibliographystyle{IEEEtranSN} +\nocite{*} +\bibliography{abernethy_bibtex} +\selectlanguage{english}% + +\end{document}