Package 'weibulltools'

Description

Provides statistical methods and visualizations that are often used in reliability engineering. Comprises a compact and easily accessible set of methods and visualization tools that make the examination and adjustment as well as the analysis and interpretation of field data (and bench tests) as simple as possible.

Besides the well-known Weibull analysis, the package supports multiple lifetime distributions and also contains Monte Carlo methods for the correction and completion of imprecisely recorded or unknown lifetime characteristics.

Plots are created statically (ggplot2) or interactively (plotly) and can be customized with functions of the respective visualization package.

Description

A dataset containing the number of cycles of fatigue life for Alloy T7987 specimens.

Usage

alloy

Format

A tibble with 72 rows and 2 variables:

cycles

Number of cycles (in thousands).

status

If specimen failed before 300 thousand cycles 1 else 0.

Source

Meeker, William Q; Escobar, Luis A., Statistical Methods for Reliability Data, New York: Wiley series in probability and statistics (1998, p.131)

Description

This function computes the non-parametric beta binomial confidence bounds (BB) for quantiles and failure probabilities.

Usage

confint_betabinom(x, ...)

## S3 method for class 'wt_model'
confint_betabinom(
  x,
  b_lives = c(0.01, 0.1, 0.5),
  bounds = c("two_sided", "lower", "upper"),
  conf_level = 0.95,
  direction = c("y", "x"),
  ...
)

Arguments

Details

The procedure is similar to the Median Ranks method but with the difference that instead of finding the probability for the j-th rank at the 50% level the probability (probabilities) has (have) to be found at the given confidence level.

Value

A tibble with class wt_confint containing the following columns:

• x : An ordered sequence of the lifetime characteristic regarding the failed units, starting at min(x) and ending up at max(x). With b_lives = c(0.01, 0.1, 0.5) the 1%, 10% and 50% quantiles are additionally included in x, but only if the specified probabilities are in the range of the estimated probabilities.

• rank : Interpolated ranks as a function of probabilities, computed with the converted approximation formula of Benard.

• prob : An ordered sequence of probabilities with specified b_lives included.

• lower_bound : Provided, if bounds is one of "two_sided" or "lower". Lower confidence limits with respect to direction, i.e. limits for quantiles or probabilities.

• upper_bound : Provided, if bounds is one of "two_sided" or "upper". Upper confidence limits with respect to direction, i.e. limits for quantiles or probabilities.

• cdf_estimation_method : Method for the estimation of failure probabilities which was specified in estimate_cdf.

Further information is stored in the attributes of this tibble:

• distribution : Distribution which was specified in rank_regression.

• bounds : Specified bound(s).

• direction : Specified direction.

• model_estimation : Input list with class wt_model.

Examples

# Reliability data preparation:
## Data for two-parametric model:
data_2p <- reliability_data(
  shock,
  x = distance,
  status = status
)

## Data for three-parametric model:
data_3p <- reliability_data(
  alloy,
  x = cycles,
  status = status
)

# Probability estimation:
prob_tbl_2p <- estimate_cdf(
  data_2p,
  methods = "johnson"
)

prob_tbl_3p <- estimate_cdf(
  data_3p,
  methods = "johnson"
)

prob_tbl_mult <- estimate_cdf(
  data_3p,
  methods = c("johnson", "mr")
)

# Model estimation with rank_regression():
rr_2p <- rank_regression(
  prob_tbl_2p,
  distribution = "weibull"
)

rr_3p <- rank_regression(
  prob_tbl_3p,
  distribution = "lognormal3",
  conf_level = 0.90
)

rr_lists <- rank_regression(
  prob_tbl_mult,
  distribution = "loglogistic3",
  conf_level = 0.90
)

# Example 1 - Two-sided 95% confidence interval for probabilities ('y'):
conf_betabin_1 <- confint_betabinom(
  x = rr_2p,
  bounds = "two_sided",
  conf_level = 0.95,
  direction = "y"
)

# Example 2 - One-sided lower/upper 90% confidence interval for quantiles ('x'):
conf_betabin_2_1 <- confint_betabinom(
  x = rr_2p,
  bounds = "lower",
  conf_level = 0.90,
  direction = "x"
)

conf_betabin_2_2 <- confint_betabinom(
  x = rr_2p,
  bounds = "upper",
  conf_level = 0.90,
  direction = "x"
)

# Example 3 - Two-sided 90% confidence intervals for both directions using
# a three-parametric model:
conf_betabin_3_1 <- confint_betabinom(
  x = rr_3p,
  bounds = "two_sided",
  conf_level = 0.90,
  direction = "y"
)

conf_betabin_3_2 <- confint_betabinom(
  x = rr_3p,
  bounds = "two_sided",
  conf_level = 0.90,
  direction = "x"
)

# Example 4 - Confidence intervals if multiple methods in estimate_cdf, i.e.
# "johnson" and "mr", were specified:

conf_betabin_4 <- confint_betabinom(
  x = rr_lists,
  bounds = "two_sided",
  conf_level = 0.99,
  direction = "y"
)

Description

This function computes the non-parametric beta binomial confidence bounds (BB) for quantiles and failure probabilities.

Usage

## Default S3 method:
confint_betabinom(
  x,
  status,
  dist_params,
  distribution = c("weibull", "lognormal", "loglogistic", "sev", "normal", "logistic",
    "weibull3", "lognormal3", "loglogistic3", "exponential", "exponential2"),
  b_lives = c(0.01, 0.1, 0.5),
  bounds = c("two_sided", "lower", "upper"),
  conf_level = 0.95,
  direction = c("y", "x"),
  ...
)

Arguments

Details

The procedure is similar to the Median Ranks method but with the difference that instead of finding the probability for the j-th rank at the 50% level the probability (probabilities) has (have) to be found at the given confidence level.

Value

A tibble with class wt_confint containing the following columns:

• x : An ordered sequence of the lifetime characteristic regarding the failed units, starting at min(x) and ending up at max(x). With b_lives = c(0.01, 0.1, 0.5) the 1%, 10% and 50% quantiles are additionally included in x, but only if the specified probabilities are in the range of the estimated probabilities.

• rank : Interpolated ranks as a function of probabilities, computed with the converted approximation formula of Benard.

• prob : An ordered sequence of probabilities with specified b_lives included.

• lower_bound : Provided, if bounds is one of "two_sided" or "lower". Lower confidence limits with respect to direction, i.e. limits for quantiles or probabilities.

• upper_bound : Provided, if bounds is one of "two_sided" or "upper". Upper confidence limits with respect to direction, i.e. limits for quantiles or probabilities.

• cdf_estimation_method : A character that is always NA_character. Only needed for internal use.

Further information is stored in the attributes of this tibble:

• distribution : Distribution which was specified in rank_regression.

• bounds : Specified bound(s).

• direction : Specified direction.

See Also

confint_betabinom

Examples

# Vectors:
obs <- seq(10000, 100000, 10000)
status_1 <- c(0, 1, 1, 0, 0, 0, 1, 0, 1, 0)

cycles <- alloy$cycles
status_2 <- alloy$status

# Probability estimation:
prob_tbl <- estimate_cdf(
  x = obs,
  status = status_1,
  method = "johnson"
)

prob_tbl_2 <- estimate_cdf(
  x = cycles,
  status = status_2,
  method = "johnson"
)

# Model estimation with rank_regression():
rr <- rank_regression(
  x = prob_tbl$x,
  y = prob_tbl$prob,
  status = prob_tbl$status,
  distribution = "weibull",
  conf_level = 0.9
)

rr_2 <- rank_regression(
  x = prob_tbl_2$x,
  y = prob_tbl_2$prob,
  status = prob_tbl_2$status,
  distribution = "lognormal3"
)

# Example 1 - Two-sided 95% confidence interval for probabilities ('y'):
conf_betabin_1 <- confint_betabinom(
  x = prob_tbl$x,
  status = prob_tbl$status,
  dist_params = rr$coefficients,
  distribution = "weibull",
  bounds = "two_sided",
  conf_level = 0.95,
  direction = "y"
)

# Example 2 - One-sided lower/upper 90% confidence interval for quantiles ('x'):
conf_betabin_2_1 <- confint_betabinom(
  x = prob_tbl$x,
  status = prob_tbl$status,
  dist_params = rr$coefficients,
  distribution = "weibull",
  bounds = "lower",
  conf_level = 0.9,
  direction = "x"
)

conf_betabin_2_2 <- confint_betabinom(
  x = prob_tbl$x,
  status = prob_tbl$status,
  dist_params = rr$coefficients,
  distribution = "weibull",
  bounds = "upper",
  conf_level = 0.9,
  direction = "x"
)

# Example 3 - Two-sided 90% confidence intervals for both directions using
# a three-parametric model:

conf_betabin_3_1 <- confint_betabinom(
  x = prob_tbl_2$x,
  status = prob_tbl_2$status,
  dist_params = rr_2$coefficients,
  distribution = "lognormal3",
  bounds = "two_sided",
  conf_level = 0.9,
  direction = "y"
)

conf_betabin_3_2 <- confint_betabinom(
  x = prob_tbl_2$x,
  status = prob_tbl_2$status,
  dist_params = rr_2$coefficients,
  distribution = "lognormal3",
  bounds = "two_sided",
  conf_level = 0.9,
  direction = "x"
)

Description

This function computes normal-approximation confidence intervals for quantiles and failure probabilities.

Usage

confint_fisher(x, ...)

## S3 method for class 'wt_model'
confint_fisher(
  x,
  b_lives = c(0.01, 0.1, 0.5),
  bounds = c("two_sided", "lower", "upper"),
  conf_level = 0.95,
  direction = c("y", "x"),
  ...
)

Arguments

Details

The basis for the calculation of these confidence bounds are the standard errors obtained by the delta method.

The bounds on the probability are determined by the z-procedure. See 'References' for more information on this approach.

Value

A tibble with class wt_confint containing the following columns:

• x : An ordered sequence of the lifetime characteristic regarding the failed units, starting at min(x) and ending up at max(x). With b_lives = c(0.01, 0.1, 0.5) the 1%, 10% and 50% quantiles are additionally included in x, but only if the specified probabilities are in the range of the estimated probabilities.

• prob : An ordered sequence of probabilities with specified b_lives included.

• std_err : Estimated standard errors with respect to direction.

• lower_bound : Provided, if bounds is one of "two_sided" or "lower". Lower confidence limits with respect to direction, i.e. limits for quantiles or probabilities.

• upper_bound : Provided, if bounds is one of "two_sided" or "upper". Upper confidence limits with respect to direction, i.e. limits for quantiles or probabilities.

• cdf_estimation_method : A character that is always NA_character. Only needed for internal use.

Further information is stored in the attributes of this tibble:

• distribution : Distribution which was specified in ml_estimation.

• bounds : Specified bound(s).

• direction : Specified direction.

• model_estimation : Input list with classes wt_model and wt_ml_estimation.

References

Meeker, William Q; Escobar, Luis A., Statistical methods for reliability data, New York: Wiley series in probability and statistics, 1998

Examples

# Reliability data preparation:
## Data for two-parametric model:
data_2p <- reliability_data(
  shock,
  x = distance,
  status = status
)

## Data for three-parametric model:
data_3p <- reliability_data(
  alloy,
  x = cycles,
  status = status
)

# Model estimation with ml_estimation():
ml_2p <- ml_estimation(
  data_2p,
  distribution = "weibull"
)

ml_3p <- ml_estimation(
  data_3p,
  distribution = "lognormal3",
  conf_level = 0.90
)


# Example 1 - Two-sided 95% confidence interval for probabilities ('y'):
conf_fisher_1 <- confint_fisher(
  x = ml_2p,
  bounds = "two_sided",
  conf_level = 0.95,
  direction = "y"
)

# Example 2 - One-sided lower/upper 90% confidence interval for quantiles ('x'):
conf_fisher_2_1 <- confint_fisher(
  x = ml_2p,
  bounds = "lower",
  conf_level = 0.90,
  direction = "x"
)

conf_fisher_2_2 <- confint_fisher(
  x = ml_2p,
  bounds = "upper",
  conf_level = 0.90,
  direction = "x"
)

# Example 3 - Two-sided 90% confidence intervals for both directions using
# a three-parametric model:

conf_fisher_3_1 <- confint_fisher(
  x = ml_3p,
  bounds = "two_sided",
  conf_level = 0.90,
  direction = "y"
)

conf_fisher_3_2 <- confint_fisher(
  x = ml_3p,
  bounds = "two_sided",
  conf_level = 0.90,
  direction = "x"
)

Description

This function computes normal-approximation confidence intervals for quantiles and failure probabilities.

Usage

## Default S3 method:
confint_fisher(
  x,
  status,
  dist_params,
  dist_varcov,
  distribution = c("weibull", "lognormal", "loglogistic", "sev", "normal", "logistic",
    "weibull3", "lognormal3", "loglogistic3", "exponential", "exponential2"),
  b_lives = c(0.01, 0.1, 0.5),
  bounds = c("two_sided", "lower", "upper"),
  conf_level = 0.95,
  direction = c("y", "x"),
  ...
)

Arguments

Details

The basis for the calculation of these confidence bounds are the standard errors obtained by the delta method.

The bounds on the probability are determined by the z-procedure. See 'References' for more information on this approach.

Value

A tibble with class wt_confint containing the following columns:

• x : An ordered sequence of the lifetime characteristic regarding the failed units, starting at min(x) and ending up at max(x). With b_lives = c(0.01, 0.1, 0.5) the 1%, 10% and 50% quantiles are additionally included in x, but only if the specified probabilities are in the range of the estimated probabilities.

• prob : An ordered sequence of probabilities with specified b_lives included.

• std_err : Estimated standard errors with respect to direction.

• lower_bound : Provided, if bounds is one of "two_sided" or "lower". Lower confidence limits with respect to direction, i.e. limits for quantiles or probabilities.

• upper_bound : Provided, if bounds is one of "two_sided" or "upper". Upper confidence limits with respect to direction, i.e. limits for quantiles or probabilities.

• cdf_estimation_method : A character that is always NA_character. Only needed for internal use.

Further information is stored in the attributes of this tibble:

• distribution : Distribution which was specified in ml_estimation.

• bounds : Specified bound(s).

• direction : Specified direction.

References

Meeker, William Q; Escobar, Luis A., Statistical methods for reliability data, New York: Wiley series in probability and statistics, 1998

See Also

confint_fisher

Examples

# Vectors:
obs <- seq(10000, 100000, 10000)
status_1 <- c(0, 1, 1, 0, 0, 0, 1, 0, 1, 0)

cycles <- alloy$cycles
status_2 <- alloy$status


# Model estimation with ml_estimation():
ml <- ml_estimation(
  x = obs,
  status = status_1,
  distribution = "weibull",
  conf_level = 0.90
)

ml_2 <- ml_estimation(
  x = cycles,
  status = status_2,
  distribution = "lognormal3"
)

# Example 1 - Two-sided 95% confidence interval for probabilities ('y'):
conf_fisher_1 <- confint_fisher(
  x = obs,
  status = status_1,
  dist_params = ml$coefficients,
  dist_varcov = ml$varcov,
  distribution = "weibull",
  bounds = "two_sided",
  conf_level = 0.95,
  direction = "y"
)

# Example 2 - One-sided lower/upper 90% confidence interval for quantiles ('x'):
conf_fisher_2_1 <- confint_fisher(
  x = obs,
  status = status_1,
  dist_params = ml$coefficients,
  dist_varcov = ml$varcov,
  distribution = "weibull",
  bounds = "lower",
  conf_level = 0.90,
  direction = "x"
)

conf_fisher_2_2 <- confint_fisher(
  x = obs,
  status = status_1,
  dist_params = ml$coefficients,
  dist_varcov = ml$varcov,
  distribution = "weibull",
  bounds = "upper",
  conf_level = 0.90,
  direction = "x"
)

# Example 3 - Two-sided 90% confidence intervals for both directions using
# a three-parametric model:

conf_fisher_3_1 <- confint_fisher(
  x = cycles,
  status = status_2,
  dist_params = ml_2$coefficients,
  dist_varcov = ml_2$varcov,
  distribution = "lognormal3",
  bounds = "two_sided",
  conf_level = 0.90,
  direction = "y"
)

conf_fisher_3_2 <- confint_fisher(
  x = cycles,
  status = status_2,
  dist_params = ml_2$coefficients,
  dist_varcov = ml_2$varcov,
  distribution = "lognormal3",
  bounds = "two_sided",
  conf_level = 0.90,
  direction = "x"
)

Description

This function applies the delta method to a parametric lifetime distribution.

Usage

delta_method(
  x,
  dist_params,
  dist_varcov,
  distribution = c("weibull", "lognormal", "loglogistic", "sev", "normal", "logistic",
    "weibull3", "lognormal3", "loglogistic3", "exponential", "exponential2"),
  direction = c("y", "x"),
  p = deprecated()
)

Arguments

Details

The delta method estimates the standard errors for quantities that can be written as non-linear functions of ML estimators. Hence, the parameters as well as the variance-covariance matrix of these quantities have to be estimated with maximum likelihood.

The estimated standard errors are used to calculate Fisher's (normal approximation) confidence intervals. For confidence bounds on the probability, standard errors of the standardized quantiles (direction = "y") have to be computed (z-procedure) and for bounds on quantiles, standard errors of quantiles (direction = "x") are required. For more information see confint_fisher.

Value

A numeric vector of estimated standard errors for quantiles or standardized quantiles (z-values).

References

Meeker, William Q; Escobar, Luis A., Statistical methods for reliability data, New York: Wiley series in probability and statistics, 1998

Examples

# Reliability data preparation:
data <- reliability_data(
  shock,
  x = distance,
  status = status
)

# Parameter estimation using maximum likelihood:
mle <- ml_estimation(
  data,
  distribution = "weibull",
  conf_level = 0.95
)

# Example 1 - Standard errors of standardized quantiles:
delta_y <- delta_method(
  x = shock$distance,
  dist_params = mle$coefficients,
  dist_varcov = mle$varcov,
  distribution = "weibull",
  direction = "y"
)

# Example 2 - Standard errors of quantiles:
delta_x <- delta_method(
  x = seq(0.01, 0.99, 0.01),
  dist_params = mle$coefficients,
  dist_varcov = mle$varcov,
  distribution = "weibull",
  direction = "x"
)

Description

This function models a delay (in days) random variable (e.g. in logistic, registration, report) using a supposed continuous distribution. First, the row-wise differences in days of the related date columns are calculated and then the parameter(s) of the assumed distribution is (are) estimated with maximum likelihood. See 'Details' for more information.

Usage

dist_delay(...)

## S3 method for class 'wt_mcs_delay_data'
dist_delay(..., x, distribution = c("lognormal", "exponential"))

Arguments

Details

The distribution parameter(s) is (are) determined on the basis of complete cases, i.e. there is no NA (row-wise) in one of the related date columns. Time differences less than or equal to zero are not considered as well.

Value

A list with class wt_delay_estimation which contains:

• coefficients : A named vector of estimated parameter(s).

• delay : A numeric vector of element-wise computed differences in days.

• distribution : Specified distribution.

If more than one delay was considered in mcs_delay_data, the resulting output is a list with class wt_delay_estimation_list. In this case each list element has class wt_delay_estimation and the items listed above, are included.

Examples

# MCS data preparation:
## Data for delay in registration:
mcs_tbl_1 <- mcs_delay_data(
  field_data,
  date_1 = production_date,
  date_2 = registration_date,
  time = dis,
  status = status,
  id = vin
)

## Data for delay in report:
mcs_tbl_2 <- mcs_delay_data(
  field_data,
  date_1 = repair_date,
  date_2 = report_date,
  time = dis,
  status = status,
  id = vin
)

## Data for both delays:
mcs_tbl_both <- mcs_delay_data(
  field_data,
  date_1 = c(production_date, repair_date),
  date_2 = c(registration_date, report_date),
  time = dis,
  status = status,
  id = vin
)

# Example 1 - Delay in registration:
params_delay_regist  <- dist_delay(
  x = mcs_tbl_1,
  distribution = "lognormal"
)

# Example 2 - Delay in report:
params_delay_report  <- dist_delay(
  x = mcs_tbl_2,
  distribution = "exponential"
)

# Example 3 - Delays in registration and report with same distribution:
params_delays  <- dist_delay(
  x = mcs_tbl_both,
  distribution = "lognormal"
)

# Example 4 - Delays in registration and report with different distributions:
params_delays_2  <- dist_delay(
  x = mcs_tbl_both,
  distribution = c("lognormal", "exponential")
)

Description

dist_delay_register() is no longer under active development, switching to dist_delay is recommended.

Usage

dist_delay_register(date_prod, date_register, distribution = "lognormal")

Arguments

Details

This function introduces a delay random variable by calculating the time difference between the registration and production date for the sample units and afterwards estimates the parameter(s) of a supposed distribution, using maximum likelihood.

Value

A named vector of estimated parameters for the specified distribution.

Examples

date_of_production   <- c("2014-07-28", "2014-02-17", "2014-07-14",
                          "2014-06-26", "2014-03-10", "2014-05-14",
                          "2014-05-06", "2014-03-07", "2014-03-09",
                          "2014-04-13", "2014-05-20", "2014-07-07",
                          "2014-01-27", "2014-01-30", "2014-03-17",
                          "2014-02-09", "2014-04-14", "2014-04-20",
                          "2014-03-13", "2014-02-23", "2014-04-03",
                          "2014-01-08", "2014-01-08")
date_of_registration <- c(NA, "2014-03-29", "2014-12-06", "2014-09-09",
                          NA, NA, "2014-06-16", NA, "2014-05-23",
                          "2014-05-09", "2014-05-31", NA, "2014-04-13",
                          NA, NA, "2014-03-12", NA, "2014-06-02",
                          NA, "2014-03-21", "2014-06-19", NA, NA)

params_delay_regist  <- dist_delay_register(
  date_prod = date_of_production,
  date_register = date_of_registration,
  distribution = "lognormal"
)

Description

dist_delay_report()is no longer under active development, switching to dist_delay is recommended.

Usage

dist_delay_report(date_repair, date_report, distribution = "lognormal")

Arguments

Details

This function introduces a delay random variable by calculating the time difference between the report and repair date for the sample units and afterwards estimates the parameter(s) of a supposed distribution, using maximum likelihood.

Value

A named vector of estimated parameters for the specified distribution.

Examples

date_of_repair <- c(NA, "2014-09-15", "2015-07-04", "2015-04-10", NA,
                    NA, "2015-04-24", NA, "2015-04-25", "2015-04-24",
                    "2015-06-12", NA, "2015-05-04", NA, NA,
                    "2015-05-22", NA, "2015-09-17", NA, "2015-08-15",
                    "2015-11-26", NA, NA)

date_of_report <- c(NA, "2014-10-09", "2015-08-28", "2015-04-15", NA,
                    NA, "2015-05-16", NA, "2015-05-28", "2015-05-15",
                    "2015-07-11", NA, "2015-08-14", NA, NA,
                    "2015-06-05", NA, "2015-10-17", NA, "2015-08-21",
                    "2015-12-02", NA, NA)

params_delay_report  <- dist_delay_report(
  date_repair = date_of_repair,
  date_report = date_of_report,
  distribution = "lognormal"
)

Description

This function models a delay (in days) random variable (e.g. in logistic, registration, report) using a supposed continuous distribution. First, the element-wise differences in days of both vectors date_1 and date_2 are calculated and then the parameter(s) of the assumed distribution is (are) estimated with maximum likelihood. See 'Details' for more information.

Usage

## Default S3 method:
dist_delay(..., date_1, date_2, distribution = c("lognormal", "exponential"))

Arguments

Details

The distribution parameter(s) is (are) determined on the basis of complete cases, i.e. there is no NA in one of the related vector elements c(date_1[i], date_2[i]). Time differences less than or equal to zero are not considered as well.

Value

A list with class wt_delay_estimation which contains:

• coefficients : A named vector of estimated parameter(s).

• delay : A numeric vector of element-wise computed differences in days.

• distribution : Specified distribution.

If more than one delay was considered, the resulting output is a list with class wt_delay_estimation_list. In this case each list element has class wt_delay_estimation and the items listed above, are included.

See Also

dist_delay

Examples

# Example 1 - Delay in registration:
params_delay_regist  <- dist_delay(
  date_1 = field_data$production_date,
  date_2 = field_data$registration_date,
  distribution = "lognormal"
)

# Example 2 - Delay in report:
params_delay_report  <- dist_delay(
  date_1 = field_data$repair_date,
  date_2 = field_data$report_date,
  distribution = "exponential"
)

# Example 3 - Delays in registration and report with same distribution:
params_delays  <- dist_delay(
  date_1 = list(field_data$production_date, field_data$repair_date),
  date_2 = list(field_data$registration_date, field_data$report_date),
  distribution = "lognormal"
)

# Example 4 - Delays in registration and report with different distributions:
params_delays_2  <- dist_delay(
  date_1 = list(field_data$production_date, field_data$repair_date),
  date_2 = list(field_data$registration_date, field_data$report_date),
  distribution = c("lognormal", "exponential")
)

Description

This function models a mileage random variable on an annual basis with respect to a supposed continuous distribution. First, the distances are calculated for one year (365 days) using a linear relationship between the distance and operating time. Second, the parameter(s) of the assumed distribution are estimated with maximum likelihood. See 'Details' for more information.

Usage

dist_mileage(x, ...)

## S3 method for class 'wt_mcs_mileage_data'
dist_mileage(x, distribution = c("lognormal", "exponential"), ...)

Arguments

Details

The distribution parameter(s) is (are) determined on the basis of complete cases, i.e. there is no NA (row-wise) in one of the related columns mileage and time. Distances and operating times less than or equal to zero are not considered as well.

Assumption of linear relationship: Imagine a component in a vehicle has endured a distance of 25000 kilometers (km) in 500 days (d), the annual distance of this unit is

$25000 km \cdot (\frac{365 d} {500 d}) = 18250 km$

Value

A list with class wt_mileage_estimation which contains:

• coefficients : A named vector of estimated parameter(s).

• miles_annual : A numeric vector of element-wise computed annual distances using the linear relationship described in 'Details'.

• distribution : Specified distribution.

Examples

# MCS data preparation:
mcs_tbl <- mcs_mileage_data(
  field_data,
  mileage = mileage,
  time = dis,
  status = status,
  id = vin
)

# Example 1 - Assuming lognormal annual mileage distribution:
params_mileage_annual <- dist_mileage(
  x = mcs_tbl,
  distribution = "lognormal"
)

# Example 2 - Assuming exponential annual mileage distribution:
params_mileage_annual_2 <- dist_mileage(
  x = mcs_tbl,
  distribution = "exponential"
)

Description

This function models a mileage random variable on an annual basis with respect to a supposed continuous distribution. First, the distances are calculated for one year (365 days) using a linear relationship between the distance and operating time. Second, the parameter(s) of the assumed distribution are estimated with maximum likelihood. See 'Details' for more information.

Usage

## Default S3 method:
dist_mileage(x, time, distribution = c("lognormal", "exponential"), ...)

Arguments

Details

The distribution parameter(s) is (are) determined on the basis of complete cases, i.e. there is no NA in one of the related vector elements c(mileage[i], time[i]). Distances and operating times less than or equal to zero are not considered as well.

Assumption of linear relationship: Imagine a component in a vehicle has endured a distance of 25000 kilometers (km) in 500 days (d), the annual distance of this unit is

$25000 km \cdot (\frac{365 d} {500 d}) = 18250 km$

Value

A list with class wt_mileage_estimation which contains:

• coefficients : A named vector of estimated parameter(s).

• miles_annual : A numeric vector of element-wise computed annual distances using the linear relationship described in 'Details'.

• distribution : Specified distribution.

See Also

dist_mileage

Examples

# Example 1 - Assuming lognormal annual mileage distribution:
params_mileage_annual <- dist_mileage(
  x = field_data$mileage,
  time = field_data$dis,
  distribution = "lognormal"
)

# Example 2 - Assuming exponential annual mileage distribution:
params_mileage_annual_2 <- dist_mileage(
  x = field_data$mileage,
  time = field_data$dis,
  distribution = "exponential"
)

Description

This function applies a non-parametric method to estimate the failure probabilities of complete data taking (multiple) right-censored observations into account.

Usage

estimate_cdf(x, ...)

## S3 method for class 'wt_reliability_data'
estimate_cdf(
  x,
  methods = c("mr", "johnson", "kaplan", "nelson"),
  options = list(),
  ...
)

Arguments

Details

One or multiple techniques can be used for the methods argument:

• "mr" : Method Median Ranks is used to estimate the failure probabilities of failed units without considering censored items. Tied observations can be handled in three ways (See 'Options'):

  • "max" : Highest observed rank is assigned to tied observations.

  • "min" : Lowest observed rank is assigned to tied observations.

  • "average" : Mean rank is assigned to tied observations.

  Two formulas can be used to determine cumulative failure probabilities F(t) (See 'Options'):

  • "benard" : Benard's approximation for Median Ranks.

  • "invbeta" : Exact Median Ranks using the inverse beta distribution.

• "johnson" : The Johnson method is used to estimate the failure probabilities of failed units, taking censored units into account. Compared to complete data, correction of probabilities is done by the computation of adjusted ranks. Two formulas can be used to determine cumulative failure probabilities F(t) (See 'Options'):

  • "benard" : Benard's approximation for Median Ranks.

  • "invbeta" : Exact Median Ranks using the inverse beta distribution.

• "kaplan" : The method of Kaplan and Meier is used to estimate the survival function S(t) with respect to (multiple) right censored data. The complement of S(t), i.e. F(t), is returned. In contrast to the original Kaplan-Meier estimator, one modification is made (see 'References').

• "nelson" : The Nelson-Aalen estimator models the cumulative hazard rate function in case of (multiple) right censored data. Equating the formal definition of the hazard rate with that according to Nelson-Aalen results in a formula for the calculation of failure probabilities.

Value

A tibble with class wt_cdf_estimation containing the following columns:

• id : Identification for every unit.

• x : Lifetime characteristic.

• status : Binary data (0 or 1) indicating whether a unit is a right censored observation (= 0) or a failure (= 1).

• rank : The (computed) ranks. Determined for methods "mr" and "johnson", filled with NA for other methods or if status = 0.

• prob : Estimated failure probabilities, NA if status = 0.

• cdf_estimation_method : Specified method for the estimation of failure probabilities.

Options

Argument options is a named list of options:

References

NIST/SEMATECH e-Handbook of Statistical Methods, 8.2.1.5. Empirical model fitting - distribution free (Kaplan-Meier) approach, NIST SEMATECH, December 3, 2020

Examples

# Reliability data:
data <- reliability_data(
  alloy,
  x = cycles,
  status = status
)

# Example 1 - Johnson method:
prob_tbl <- estimate_cdf(
  x = data,
  methods = "johnson"
)

# Example 2 - Multiple methods:
prob_tbl_2 <- estimate_cdf(
  x = data,
  methods = c("johnson", "kaplan", "nelson")
)

# Example 3 - Method 'mr' with options:
prob_tbl_3 <- estimate_cdf(
  x = data,
  methods = "mr",
  options = list(
    mr_method = "invbeta",
    mr_ties.method = "average"
  )
)

# Example 4 - Multiple methods and options:
prob_tbl_4 <- estimate_cdf(
  x = data,
  methods = c("mr", "johnson"),
  options = list(
    mr_ties.method = "max",
    johnson_method = "invbeta"
  )
)

Description

This function applies a non-parametric method to estimate the failure probabilities of complete data taking (multiple) right-censored observations into account.

Usage

## Default S3 method:
estimate_cdf(
  x,
  status,
  id = NULL,
  method = c("mr", "johnson", "kaplan", "nelson"),
  options = list(),
  ...
)

Arguments

Details

The following techniques can be used for the method argument:

• "mr" : Method Median Ranks is used to estimate the failure probabilities of failed units without considering censored items. Tied observations can be handled in three ways (See 'Options'):

  • "max" : Highest observed rank is assigned to tied observations.

  • "min" : Lowest observed rank is assigned to tied observations.

  • "average" : Mean rank is assigned to tied observations.

  Two formulas can be used to determine cumulative failure probabilities F(t) (See 'Options'):

  • "benard" : Benard's approximation for Median Ranks.

  • "invbeta" : Exact Median Ranks using the inverse beta distribution.

• "johnson" : The Johnson method is used to estimate the failure probabilities of failed units, taking censored units into account. Compared to complete data, correction of probabilities is done by the computation of adjusted ranks. Two formulas can be used to determine cumulative failure probabilities F(t) (See 'Options'):

  • "benard" : Benard's approximation for Median Ranks.

  • "invbeta" : Exact Median Ranks using the inverse beta distribution.

• "kaplan" : The method of Kaplan and Meier is used to estimate the survival function S(t) with respect to (multiple) right censored data. The complement of S(t), i.e. F(t), is returned. In contrast to the original Kaplan-Meier estimator, one modification is made (see 'References').

• "nelson" : The Nelson-Aalen estimator models the cumulative hazard rate function in case of (multiple) right censored data. Equating the formal definition of the hazard rate with that according to Nelson-Aalen results in a formula for the calculation of failure probabilities.

Value

A tibble with class wt_cdf_estimation containing the following columns:

• id : Identification for every unit.

• x : Lifetime characteristic.

• status : Binary data (0 or 1) indicating whether a unit is a right censored observation (= 0) or a failure (= 1).

• rank : The (computed) ranks. Determined for methods "mr" and "johnson", filled with NA for other methods or if status = 0.

• prob : Estimated failure probabilities, NA if status = 0.

• cdf_estimation_method : Specified method for the estimation of failure probabilities.

Options

Argument options is a named list of options:

References

NIST/SEMATECH e-Handbook of Statistical Methods, 8.2.1.5. Empirical model fitting - distribution free (Kaplan-Meier) approach, NIST SEMATECH, December 3, 2020

See Also

estimate_cdf

Examples

# Vectors:
cycles <- alloy$cycles
status <- alloy$status

# Example 1 - Johnson method:
prob_tbl <- estimate_cdf(
  x = cycles,
  status = status,
  method = "johnson"
)


# Example 2 - Method 'mr' with options:
prob_tbl_2 <- estimate_cdf(
  x = cycles,
  status = status,
  method = "mr",
  options = list(
    mr_method = "invbeta",
    mr_ties.method = "average"
  )
)

Description

An illustrative field dataset that contains a variety of variables commonly collected in the automotive sector.

The dataset has complete information about failed and incomplete information about intact vehicles. See 'Format' and 'Details' for further insights.

Usage

field_data

Format

A tibble with 10,684 rows and 20 variables:

vin

Vehicle identification number.

dis

Days in service.

mileage

Distances covered, which are unknown for censored units.

status

1 for failed and 0 for censored units.

production_date

Date of production.

registration_date

Date of registration. Known for all failed units and for a few intact units.

repair_date

The date on which the failure was repaired. It is assumed that the repair date is equal to the date of failure occurrence.

report_date

The date on which lifetime information about the failure were available.

country

Delivering country.

region

The region within the country of delivery. Known for registered vehicles, NA for units with a missing registration_date.

climatic_zone

Climatic zone based on "Köppen-Geiger" climate classification. Known for registered vehicles, NA for units with a missing registration_date.

climatic_subzone

Climatic subzone based on "Köppen-Geiger" climate classification. Known for registered vehicles, NA for units with a registration_date.

brand

Brand of the vehicle.

vehicle_model

Model of the vehicle.

engine_type

Type of the engine.

engine_date

Date where the engine was installed.

gear_type

Type of the gear.

gear_date

Date where the gear was installed.

transmission

Transmission of the vehicle.

fuel

Vehicle fuel.

Details

All vehicles were produced in 2014 and an analysis of the field data was made at the end of 2015. At the date of analysis, there were 684 failed and 10,000 intact vehicles.

Censored vehicles:

For censored units the service time (dis) was computed as the difference of the date of analysis "2015-12-31" and the registration_date.

For many units the latter date is unknown. For these, the difference of the analysis date and production_date was used to get a rough estimation of the real service time. This uncertainty has to be considered in the subsequent analysis (see delay in registration in the section 'Details' of mcs_delay).

Furthermore, due to the delay in report, the computed service time could also be inaccurate. This uncertainty should be considered as well (see delay in report in the section 'Details' of mcs_delay).

The lifetime characteristic mileage is unknown for all censored units. If an analysis is to be made for this lifetime characteristic, covered distances for these units have to be estimated (see mcs_mileage).

Failed vehicles: For failed units the service time (dis) is computed as the difference of repair_date and registration_date, which are known for all of them.

See Also

mcs_mileage_data

Description

johnson_method() is no longer under active development, switching to estimate_cdf is recommended.

Usage

johnson_method(x, status, id = NULL, method = c("benard", "invbeta"))

Arguments

Details

This non-parametric approach is used to estimate the failure probabilities in terms of uncensored or (multiple) right censored data. Compared to complete data the correction is done by calculating adjusted ranks which takes non-defective units into account.

Value

A tibble containing the following columns:

• id : Identification for every unit.

• x : Lifetime characteristic.

• status : Binary data (0 or 1) indicating whether a unit is a right censored observation (= 0) or a failure (= 1).

• rank : Adjusted ranks, NA if status = 0.

• prob : Estimated failure probabilities, NA if status = 0.

• cdf_estimation_method : Specified method for the estimation of failure probabilities (always 'johnson').

Examples

# Vectors:
obs   <- seq(10000, 100000, 10000)
state <- c(0, 1, 1, 0, 0, 0, 1, 0, 1, 0)
uic   <- c("3435", "1203", "958X", "XX71", "abcd", "tz46",
           "fl29", "AX23", "Uy12", "kl1a")

# Example 1 - Johnson method for intact and failed units:
tbl_john <- johnson_method(
  x = obs,
  status = state,
  id = uic
)

# Example 2 - Johnson's method works also if only defective units are considered:
tbl_john_2 <- johnson_method(
  x = obs,
  status = rep(1, length(obs))
)

Description

kaplan_method() is no longer under active development, switching to estimate_cdf is recommended.

Usage

kaplan_method(x, status, id = NULL)

Arguments

Details

Whereas the non-parametric Kaplan-Meier estimator is used to estimate the survival function S(t) in terms of (multiple) right censored data, the complement is an estimate of the cumulative distribution function F(t). One modification is made in contrast to the original Kaplan-Meier estimator (see 'References').

Value

A tibble containing the following columns:

• id : Identification for every unit.

• x : Lifetime characteristic.

• status : Binary data (0 or 1) indicating whether a unit is a right censored observation (= 0) or a failure (= 1).

• rank : Filled with NA.

• prob : Estimated failure probabilities, NA if status = 0.

• cdf_estimation_method : Specified method for the estimation of failure probabilities (always 'kaplan').

References

NIST/SEMATECH e-Handbook of Statistical Methods, 8.2.1.5. Empirical model fitting - distribution free (Kaplan-Meier) approach, NIST SEMATECH, December 3, 2020

Examples

# Vectors:
obs   <- seq(10000, 100000, 10000)
state <- c(0, 1, 1, 0, 0, 0, 1, 0, 1, 0)
state_2 <- c(0, 1, 1, 0, 0, 0, 1, 0, 0, 1)
uic   <- c("3435", "1203", "958X", "XX71", "abcd", "tz46",
           "fl29", "AX23","Uy12", "kl1a")

# Example 1 - Observation with highest characteristic is an intact unit:
tbl_kap <- kaplan_method(
  x = obs,
  status = state,
  id = uic
)

# Example 2 - Observation with highest characteristic is a defective unit:
tbl_kap_2 <- kaplan_method(
  x = obs,
  status = state_2
)

Description

This function computes the log-likelihood value with respect to a given set of parameters. In terms of Maximum Likelihood Estimation this function can be optimized (optim) to estimate the parameters and variance-covariance matrix of the parameters.

Usage

loglik_function(x, ...)

## S3 method for class 'wt_reliability_data'
loglik_function(
  x,
  wts = rep(1, nrow(x)),
  dist_params,
  distribution = c("weibull", "lognormal", "loglogistic", "sev", "normal", "logistic",
    "weibull3", "lognormal3", "loglogistic3", "exponential", "exponential2"),
  ...
)

Arguments

Value

Returns the log-likelihood value for the parameters in dist_params given the data.

Distributions

The following table summarizes the available distributions and their parameters

• location parameter $\mu$,

• scale parameter $\sigma$ or $\theta$ and

• threshold parameter $\gamma$.

The order within dist_params is given in the table header.

References

Meeker, William Q; Escobar, Luis A., Statistical methods for reliability data, New York: Wiley series in probability and statistics, 1998

Examples

# Reliability data preparation:
data <- reliability_data(
  alloy,
  x = cycles,
  status = status
)

# Example 1 - Evaluating Log-Likelihood function of two-parametric weibull:
loglik_weib <- loglik_function(
  x = data,
  dist_params = c(5.29, 0.33),
  distribution = "weibull"
)

# Example 2 - Evaluating Log-Likelihood function of three-parametric weibull:
loglik_weib3 <- loglik_function(
  x = data,
  dist_params = c(4.54, 0.76, 92.99),
  distribution = "weibull3"
)

Description

This function computes the log-likelihood value with respect to a given set of parameters. In terms of Maximum Likelihood Estimation this function can be optimized (optim) to estimate the parameters and variance-covariance matrix of the parameters.

Usage

## Default S3 method:
loglik_function(
  x,
  status,
  wts = rep(1, length(x)),
  dist_params,
  distribution = c("weibull", "lognormal", "loglogistic", "sev", "normal", "logistic",
    "weibull3", "lognormal3", "loglogistic3", "exponential", "exponential2"),
  ...
)

Arguments

Value

Returns the log-likelihood value for the parameters in dist_params given the data.

Distributions

The following table summarizes the available distributions and their parameters

• location parameter $\mu$,

• scale parameter $\sigma$ or $\theta$ and

• threshold parameter $\gamma$.

The order within dist_params is given in the table header.

References

Meeker, William Q; Escobar, Luis A., Statistical methods for reliability data, New York: Wiley series in probability and statistics, 1998

See Also

loglik_function

Examples

# Vectors:
cycles <- alloy$cycles
status <- alloy$status

# Example 1 - Evaluating Log-Likelihood function of two-parametric weibull:
loglik_weib <- loglik_function(
  x = cycles,
  status = status,
  dist_params = c(5.29, 0.33),
  distribution = "weibull"
)

# Example 2 - Evaluating Log-Likelihood function of three-parametric weibull:
loglik_weib3 <- loglik_function(
  x = cycles,
  status = status,
  dist_params = c(4.54, 0.76, 92.99),
  distribution = "weibull3"
)

Description

This function evaluates the log-likelihood with respect to a given threshold parameter of a parametric lifetime distribution. In terms of Maximum Likelihood Estimation this function can be optimized (optim) to estimate the threshold parameter.

Usage

loglik_profiling(x, ...)

## S3 method for class 'wt_reliability_data'
loglik_profiling(
  x,
  wts = rep(1, nrow(x)),
  thres,
  distribution = c("weibull3", "lognormal3", "loglogistic3", "exponential2"),
  ...
)

Arguments

Value

Returns the log-likelihood value for the threshold parameter thres given the data.

References

Meeker, William Q; Escobar, Luis A., Statistical methods for reliability data, New York: Wiley series in probability and statistics, 1998

Examples

# Reliability data preparation:
data <- reliability_data(
  alloy,
  x = cycles,
  status = status
)

# Determining the optimal loglikelihood value:
## Range of threshold parameter must be smaller than the first failure:
threshold <- seq(
  0,
  min(data$x[data$status == 1]) - 0.1,
  length.out = 50
)

## loglikelihood value with respect to threshold values:
profile_logL <- loglik_profiling(
  x = data,
  thres = threshold,
  distribution = "weibull3"
)

## Threshold value (among the candidates) that maximizes the
## loglikelihood:
threshold[which.max(profile_logL)]

## plot:
plot(
  threshold,
  profile_logL,
  type = "l"
)
abline(
  v = threshold[which.max(profile_logL)],
  h = max(profile_logL),
  col = "red"
)

Description

This function evaluates the log-likelihood with respect to a given threshold parameter of a parametric lifetime distribution. In terms of Maximum Likelihood Estimation this function can be optimized (optim) to estimate the threshold parameter.

Usage

## Default S3 method:
loglik_profiling(
  x,
  status,
  wts = rep(1, length(x)),
  thres,
  distribution = c("weibull3", "lognormal3", "loglogistic3", "exponential2"),
  ...
)

Arguments

Value

Returns the log-likelihood value for the threshold parameter thres given the data.

References

Meeker, William Q; Escobar, Luis A., Statistical methods for reliability data, New York: Wiley series in probability and statistics, 1998

See Also

loglik_profiling

Examples

# Vectors:
cycles <- alloy$cycles
status <- alloy$status

# Determining the optimal loglikelihood value:
## Range of threshold parameter must be smaller than the first failure:
threshold <- seq(
  0,
  min(cycles[status == 1]) - 0.1,
  length.out = 50
)

## loglikelihood value with respect to threshold values:
profile_logL <- loglik_profiling(
  x = cycles,
  status = status,
  thres = threshold,
  distribution = "weibull3"
)

## Threshold value (among the candidates) that maximizes the
## loglikelihood:
threshold[which.max(profile_logL)]

## plot:
plot(
  threshold,
  profile_logL,
  type = "l"
)
abline(
  v = threshold[which.max(profile_logL)],
  h = max(profile_logL),
  col = "red"
)

Description

In general, the amount of available information about units in the field is very different. During the warranty period, there are only a few cases with complete data (mainly failed units) but lots of cases with incomplete data (usually censored units). As a result, the operating time of units with incomplete information is often inaccurate and must be adjusted by delays.

This function reduces the operating times of incomplete observations by simulated delays (in days). A unit is considered as incomplete if the later of the related dates is unknown. See 'Details' for some practical examples.

Random delay numbers are drawn from the distribution determined by complete cases (described in 'Details' of dist_delay).

Usage

mcs_delay(...)

## S3 method for class 'wt_mcs_delay_data'
mcs_delay(..., x, distribution = c("lognormal", "exponential"))

Arguments

Details

In field data analysis time-dependent characteristics (e.g. time in service) are often imprecisely recorded. These inaccuracies are caused by unconsidered delays.

For a better understanding of the MCS application in the context of field data, two cases are described below.

• Delay in registration: It is common that a supplier, which provides parts to the manufacturing industry does not know when the unit, in which its parts are installed, were put in service (due to unknown registration or sales date (date_2)). Without taking the described delay into account, the time in service of the failed units would be the difference between the repair date and the production date (date_1) and for intact units the difference between the present date and the production date. But the real operating times are (much) shorter, since the stress on the components have not started until the whole systems were put in service. Hence, units with incomplete data (missing date_2) must be reduced by the delays.

• Delay in report:: Authorized repairers often do not immediately notify the manufacturer or OEM of repairs that were made during the warranty period, but instead pass the information about these repairs in collected forms e.g. weekly, monthly or quarterly. The resulting time difference between the reporting (date_2) of the repair in the guarantee database and the actual repair date (date_1), which is often assumed to be the failure date, is called the reporting delay. For a given date where the analysis is made there could be units which had a failure but the failure isn't reported and therefore they are treated as censored units. In order to take this into account and according to the principle of equal opportunities, the lifetime of units with missing report date (date_2[i] = NA) is reduced by simulated reporting delays.

Value

A list with class wt_mcs_delay containing the following elements:

• data : A tibble returned by mcs_delay_data where two modifications has been made:

  • If the column status exists, the tibble has additional classes wt_mcs_data and wt_reliability_data. Otherwise, the tibble only has the additional class wt_mcs_data (which is not supported by estimate_cdf).

  • The column time is renamed to x (to be in accordance with reliability_data) and contains the adjusted operating times for incomplete observations and input operating times for the complete observations.

• sim_data : A tibble with column sim_delay that holds the simulated delay-specific numbers for incomplete cases and 0 for complete cases. If more than one delay was considered multiple columns with names sim_delay_1, sim_delay_2, ..., sim_delay_i and corresponding delay-specific random numbers are presented.

• model_estimation : A list returned by dist_delay.

References

Verband der Automobilindustrie e.V. (VDA); Qualitätsmanagement in der Automobilindustrie. Zuverlässigkeitssicherung bei Automobilherstellern und Lieferanten. Zuverlässigkeits-Methoden und -Hilfsmittel.; 4th Edition, 2016, ISSN:0943-9412

See Also

dist_delay for the determination of a parametric delay distribution and estimate_cdf for the estimation of failure probabilities.

Examples

# MCS data preparation:
## Data for delay in registration:
mcs_tbl_1 <- mcs_delay_data(
  field_data,
  date_1 = production_date,
  date_2 = registration_date,
  time = dis,
  status = status,
  id = vin
)

## Data for delay in report:
mcs_tbl_2 <- mcs_delay_data(
  field_data,
  date_1 = repair_date,
  date_2 = report_date,
  time = dis,
  status = status,
  id = vin
)

## Data for both delays:
mcs_tbl_both <- mcs_delay_data(
  field_data,
  date_1 = c(production_date, repair_date),
  date_2 = c(registration_date, report_date),
  time = dis,
  status = status,
  id = vin
)

# Example 1 - MCS for delay in registration:
mcs_regist <- mcs_delay(
  x = mcs_tbl_1,
  distribution = "lognormal"
)

# Example 2 - MCS for delay in report:
mcs_report <- mcs_delay(
  x = mcs_tbl_2,
  distribution = "exponential"
)

# Example 3 - Reproducibility of random numbers:
set.seed(1234)
mcs_report_reproduce <- mcs_delay(
  x = mcs_tbl_2,
  distribution = "exponential"
)

# Example 4 - MCS for delays in registration and report with same distribution:
mcs_delays <- mcs_delay(
  x = mcs_tbl_both,
  distribution = "lognormal"
)

# Example 5 - MCS for delays in registration and report with different distributions:
## Assuming lognormal registration and exponential reporting delays.
mcs_delays_2 <- mcs_delay(
  x = mcs_tbl_both,
  distribution = c("lognormal", "exponential")
)

Description

Create consistent mcs_delay_data based on an existing data.frame (preferred) or on multiple equal length vectors.

Usage

mcs_delay_data(
  data = NULL,
  date_1,
  date_2,
  time,
  status = NULL,
  id = NULL,
  .keep_all = FALSE
)

Arguments

Value

A tibble with class wt_mcs_delay_data that is formed for the downstream Monte Carlo method mcs_delay. It contains the following columns (if .keep_all = FALSE):

• Column(s) preserving the input of date_1. For the vector-based approach with unnamed input, column name(s) is (are) date_1 (date_1.1, date_1.2, ..., date_1.i).

• Column(s) preserving the input of date_2. For the vector-based approach with unnamed input, column name(s) is (are) date_2 (date_2.1, date_2.2, ..., date_2.i).

• time : Input operating times.

• status (optional) :

  • If is.null(status) column status does not exist.

  • If status is provided the column contains the entered binary data (0 or 1).

• id : Identification for every unit.

If .keep_all = TRUE, the remaining columns of data are also preserved.

The attributes mcs_start_dates and mcs_end_dates hold the name(s) of the column(s) that preserve the input of date_1 and date_2.

See Also

dist_delay for the determination of a parametric delay distribution and mcs_delay for the Monte Carlo method with respect to delays.

Examples

# Example 1 -  Based on an existing data.frame/tibble and column names:
mcs_tbl <- mcs_delay_data(
  data = field_data,
  date_1 = production_date,
  date_2 = registration_date,
  time = dis,
  status = status
)

# Example 2 - Based on an existing data.frame/tibble and column positions:
mcs_tbl_2 <- mcs_delay_data(
  data = field_data,
  date_1 = 7,
  date_2 = 8,
  time = 2,
  id = 1
)

# Example 3 - Keep all variables of the tibble/data.frame entered to argument data:
mcs_tbl_3 <- mcs_delay_data(
  data = field_data,
  date_1 = production_date,
  date_2 = registration_date,
  time = dis,
  status = status,
  id = vin,
  .keep_all = TRUE
)

# Example 4 - For multiple delays (data-based):
mcs_tbl_4 <- mcs_delay_data(
  data = field_data,
  date_1 = c(production_date, repair_date),
  date_2 = c(registration_date, report_date),
  time = dis,
  status = status
)

# Example 5 - Based on vectors:
mcs_tbl_5 <- mcs_delay_data(
  date_1 = field_data$production_date,
  date_2 = field_data$registration_date,
  time = field_data$dis,
  status = field_data$status,
  id = field_data$vin
)

# Example 6 - For multiple delays (vector-based):
mcs_tbl_6 <- mcs_delay_data(
  date_1 = list(field_data$production_date, field_data$repair_date),
  date_2 = list(field_data$registration_date, field_data$report_date),
  time = field_data$dis,
  status = field_data$status,
  id = field_data$vin
)

# Example 7 - For multiple delays (vector-based with named dates):
mcs_tbl_7 <- mcs_delay_data(
  date_1 = list(d11 = field_data$production_date, d12 = field_data$repair_date),
  date_2 = list(d21 = field_data$registration_date, d22 = field_data$report_date),
  time = field_data$dis,
  status = field_data$status,
  id = field_data$vin
)

Description

mcs_delay_register() is no longer under active development, switching to mcs_delay is recommended.

Usage

mcs_delay_register(
  date_prod,
  date_register,
  time,
  status,
  distribution = "lognormal",
  details = FALSE
)

Arguments

Details

In general the amount of information about units in the field, that have not failed yet, are rare. For example it is common that a supplier, who provides parts to the automotive industry does not know when a vehicle was put in service and therefore does not know the exact operating time of the supplied parts. This function uses a Monte Carlo approach for simulating the operating times of (multiple) right censored observations, taking account of registering delays. The simulation is based on the distribution of operating times that were calculated from complete data (see dist_delay_register).

Value

A numeric vector of corrected operating times for the censored units and the input operating times for the failed units if details = FALSE. If details = TRUE the output is a list which consists of the following elements:

• time : Numeric vector of corrected operating times for the censored observations and input operating times for failed units.

• x_sim : Simulated random numbers of specified distribution with estimated parameters. The length of x_sim is equal to the number of censored observations.

• coefficients : Estimated coefficients of supposed distribution.

Examples

date_of_production   <- c("2014-07-28", "2014-02-17", "2014-07-14",
                          "2014-06-26", "2014-03-10", "2014-05-14",
                          "2014-05-06", "2014-03-07", "2014-03-09",
                          "2014-04-13", "2014-05-20", "2014-07-07",
                          "2014-01-27", "2014-01-30", "2014-03-17",
                          "2014-02-09", "2014-04-14", "2014-04-20",
                          "2014-03-13", "2014-02-23", "2014-04-03",
                          "2014-01-08", "2014-01-08")
date_of_registration <- c(NA, "2014-03-29", "2014-12-06", "2014-09-09",
                          NA, NA, "2014-06-16", NA, "2014-05-23",
                          "2014-05-09", "2014-05-31", NA, "2014-04-13",
                          NA, NA, "2014-03-12", NA, "2014-06-02",
                          NA, "2014-03-21", "2014-06-19", NA, NA)

op_time <- rep(1000, length(date_of_production))
status <- c(0, 1, 1, 1, 0, 0, 1, 0, 1, 1, 1, 0, 1, 0, 0, 1, 0, 1, 0, 1, 1, 0, 0)

# Example 1 - Simplified vector output:
x_corrected <- mcs_delay_register(
  date_prod = date_of_production,
  date_register = date_of_registration,
  time = op_time,
  status = status,
  distribution = "lognormal",
  details = FALSE
)

# Example 2 - Detailed list output:
list_detail <- mcs_delay_register(
  date_prod = date_of_production,
  date_register = date_of_registration,
  time = op_time,
  status = status,
  distribution = "lognormal",
  details = TRUE
)

Description

mcs_delay_report() is no longer under active development, switching to mcs_delay is recommended.

Usage

mcs_delay_report(
  date_repair,
  date_report,
  time,
  status,
  distribution = "lognormal",
  details = FALSE
)

Arguments

Details

The delay in report describes the time between the occurrence of a damage and the registration in the warranty database. For a given date where the analysis is made there could be units which had a failure but are not registered in the database and therefore treated as censored units. To overcome this problem this function uses a Monte Carlo approach for simulating the operating times of (multiple) right censored observations, taking account of reporting delays. The simulation is based on the distribution of operating times that were calculated from complete data, i.e. failed items (see dist_delay_report).

Value

A numeric vector of corrected operating times for the censored units and the input operating times for the failed units if details = FALSE. If details = TRUE the output is a list which consists of the following elements:

• time : Numeric vector of corrected operating times for the censored observations and input operating times for failed units.

• x_sim : Simulated random numbers of specified distribution with estimated parameters. The length of x_sim is equal to the number of censored observations.

• coefficients : Estimated coefficients of supposed distribution.

Examples

date_of_repair <- c(NA, "2014-09-15", "2015-07-04", "2015-04-10", NA,
                   NA, "2015-04-24", NA, "2015-04-25", "2015-04-24",
                    "2015-06-12", NA, "2015-05-04", NA, NA,
                    "2015-05-22", NA, "2015-09-17", NA, "2015-08-15",
                    "2015-11-26", NA, NA)

date_of_report <- c(NA, "2014-10-09", "2015-08-28", "2015-04-15", NA,
                    NA, "2015-05-16", NA, "2015-05-28", "2015-05-15",
                    "2015-07-11", NA, "2015-08-14", NA, NA,
                    "2015-06-05", NA, "2015-10-17", NA, "2015-08-21",
                    "2015-12-02", NA, NA)

op_time <- rep(1000, length(date_of_repair))
status <- c(0, 1, 1, 1, 0, 0, 1, 0, 1, 1, 1, 0, 1, 0, 0, 1, 0, 1, 0, 1, 1, 0, 0)

# Example 1 - Simplified vector output:
x_corrected <- mcs_delay_report(
  date_repair = date_of_repair,
  date_report = date_of_report,
  time = op_time,
  status = status,
  distribution = "lognormal",
  details = FALSE
)

# Example 2 - Detailed list output:
list_detail <- mcs_delay_report(
  date_repair = date_of_repair,
  date_report = date_of_report,
  time = op_time,
  status = status,
  distribution = "lognormal",
  details = TRUE
)

Description

In general, the amount of available information about units in the field is very different. During the warranty period, there are only a few cases with complete data (mainly failed units) but lots of cases with incomplete data (usually censored units). As a result, the operating time of units with incomplete information is often inaccurate and must be adjusted by delays.

This function reduces the operating times of incomplete observations by simulated delays (in days). A unit is considered as incomplete if the later of the related dates is unknown. See 'Details' for some practical examples.

Random delay numbers are drawn from the distribution determined by complete cases (described in 'Details' of dist_delay).

Usage

## Default S3 method:
mcs_delay(
  ...,
  date_1,
  date_2,
  time,
  status = NULL,
  id = paste0("ID", seq_len(length(time))),
  distribution = c("lognormal", "exponential")
)

Arguments

Details

In field data analysis time-dependent characteristics (e.g. time in service) are often imprecisely recorded. These inaccuracies are caused by unconsidered delays.

For a better understanding of the MCS application in the context of field data, two cases are described below.

• Delay in registration: It is common that a supplier, which provides parts to the manufacturing industry does not know when the unit, in which its parts are installed, were put in service (due to unknown registration or sales date (date_2)). Without taking the described delay into account, the time in service of the failed units would be the difference between the repair date and the production date (date_1) and for intact units the difference between the present date and the production date. But the real operating times are (much) shorter, since the stress on the components have not started until the whole systems were put in service. Hence, units with incomplete data (missing date_2) must be reduced by the delays.

• Delay in report:: Authorized repairers often do not immediately notify the manufacturer or OEM of repairs that were made during the warranty period, but instead pass the information about these repairs in collected forms e.g. weekly, monthly or quarterly. The resulting time difference between the reporting (date_2) of the repair in the guarantee database and the actual repair date (date_1), which is often assumed to be the failure date, is called the reporting delay. For a given date where the analysis is made there could be units which had a failure but the failure isn't reported and therefore they are treated as censored units. In order to take this into account and according to the principle of equal opportunities, the lifetime of units with missing report date (date_2[i] = NA) is reduced by simulated reporting delays.

Value

A list with class wt_mcs_delay containing the following elements:

• data : A tibble returned by mcs_delay_data where two modifications has been made:

  • If the column status exists, the tibble has additional classes wt_mcs_data and wt_reliability_data. Otherwise, the tibble only has the additional class wt_mcs_data (which is not supported by estimate_cdf).

  • The column time is renamed to x (to be in accordance with reliability_data) and contains the adjusted operating times for incomplete observations and input operating times for the complete observations.

• sim_data : A tibble with column sim_delay that holds the simulated delay-specific numbers for incomplete cases and 0 for complete cases. If more than one delay was considered multiple columns with names sim_delay_1, sim_delay_2, ..., sim_delay_i and corresponding delay-specific random numbers are presented.

• model_estimation : A list returned by dist_delay.

References

Verband der Automobilindustrie e.V. (VDA); Qualitätsmanagement in der Automobilindustrie. Zuverlässigkeitssicherung bei Automobilherstellern und Lieferanten. Zuverlässigkeits-Methoden und -Hilfsmittel.; 4th Edition, 2016, ISSN:0943-9412

See Also

dist_delay for the determination of a parametric delay distribution and estimate_cdf for the estimation of failure probabilities.

Examples

# Example 1 - MCS for delay in registration:
mcs_regist <- mcs_delay(
  date_1 = field_data$production_date,
  date_2 = field_data$registration_date,
  time = field_data$dis,
  status = field_data$status,
  distribution = "lognormal"
)

# Example 2 - MCS for delay in report:
mcs_report <- mcs_delay(
  date_1 = field_data$repair_date,
  date_2 = field_data$report_date,
  time = field_data$dis,
  status = field_data$status,
  distribution = "exponential"
)

# Example 3 - Reproducibility of random numbers:
set.seed(1234)
mcs_report_reproduce <- mcs_delay(
  date_1 = field_data$repair_date,
  date_2 = field_data$report_date,
  time = field_data$dis,
  status = field_data$status,
  distribution = "exponential"
)

# Example 4 - MCS for delays in registration and report with same distribution:
mcs_delays <- mcs_delay(
  date_1 = list(field_data$production_date, field_data$repair_date),
  date_2 = list(field_data$registration_date, field_data$report_date),
  time = field_data$dis,
  status = field_data$status,
  distribution = "lognormal"
)

# Example 5 - MCS for delays in registration and report with different distributions:
## Assuming lognormal registration and exponential reporting delays.
mcs_delays_2 <- mcs_delay(
  date_1 = list(field_data$production_date, field_data$repair_date),
  date_2 = list(field_data$registration_date, field_data$report_date),
  time = field_data$dis,
  status = field_data$status,
  distribution = c("lognormal", "exponential")
)

Description

mcs_delays() is no longer under active development, switching to mcs_delay is recommended.

Usage

mcs_delays(
  date_prod,
  date_register,
  date_repair,
  date_report,
  time,
  status,
  distribution = "lognormal",
  details = FALSE
)

Arguments

Details

This function is a wrapper that combines both, mcs_delay_register and mcs_delay_report functions for the adjustment of operating times of censored units.

Value

A numerical vector of corrected operating times for the censored units and the input operating times for the failed units if details = FALSE. If details = TRUE the output is a list which consists of the following elements:

• time : A numeric vector of corrected operating times for the censored observations and input operating times for failed units.

• x_sim_regist : Simulated random numbers of specified distribution with estimated parameters for delay in registration. The length of x_sim_regist is equal to the number of censored observations.

• x_sim_report : Simulated random numbers of specified distribution with estimated parameters for delay in report. The length of x_sim_report is equal to the number of censored observations.

• coefficients_regist : Estimated coefficients of supposed distribution for delay in registration.

• coefficients_report : Estimated coefficients of supposed distribution for delay in report

Examples

date_of_production   <- c("2014-07-28", "2014-02-17", "2014-07-14",
                          "2014-06-26", "2014-03-10", "2014-05-14",
                          "2014-05-06", "2014-03-07", "2014-03-09",
                          "2014-04-13", "2014-05-20", "2014-07-07",
                          "2014-01-27", "2014-01-30", "2014-03-17",
                          "2014-02-09", "2014-04-14", "2014-04-20",
                          "2014-03-13", "2014-02-23", "2014-04-03",
                          "2014-01-08", "2014-01-08")
date_of_registration <- c("2014-08-17", "2014-03-29", "2014-12-06",
                          "2014-09-09", "2014-05-14", "2014-07-01",
                          "2014-06-16", "2014-04-03", "2014-05-23",
                          "2014-05-09", "2014-05-31", "2014-08-12",
                          "2014-04-13", "2014-02-15", "2014-07-07",
                          "2014-03-12", "2014-05-27", "2014-06-02",
                          "2014-05-20", "2014-03-21", "2014-06-19",
                          "2014-02-12", "2014-03-27")
date_of_repair <- c(NA, "2014-09-15", "2015-07-04", "2015-04-10", NA,
                   NA, "2015-04-24", NA, "2015-04-25", "2015-04-24",
                    "2015-06-12", NA, "2015-05-04", NA, NA,
                    "2015-05-22", NA, "2015-09-17", NA, "2015-08-15",
                    "2015-11-26", NA, NA)

date_of_report <- c(NA, "2014-10-09", "2015-08-28", "2015-04-15", NA,
                    NA, "2015-05-16", NA, "2015-05-28", "2015-05-15",
                    "2015-07-11", NA, "2015-08-14", NA, NA,
                    "2015-06-05", NA, "2015-10-17", NA, "2015-08-21",
                    "2015-12-02", NA, NA)

op_time <- rep(1000, length(date_of_repair))
status <- c(0, 1, 1, 1, 0, 0, 1, 0, 1, 1, 1, 0, 1, 0, 0, 1, 0, 1, 0, 1, 1, 0, 0)

# Example 1 - Simplified vector output:
x_corrected <- mcs_delays(
  date_prod = date_of_production,
  date_register = date_of_registration,
  date_repair = date_of_repair,
  date_report = date_of_report,
  time = op_time,
  status = status,
  distribution = "lognormal",
  details = FALSE
)

# Example 2 - Detailed list output:
list_detail <- mcs_delays(
  date_prod = date_of_production,
  date_register = date_of_registration,
  date_repair = date_of_repair,
  date_report = date_of_report,
  time = op_time,
  status = status,
  distribution = "lognormal",
  details = TRUE
)

Description

This function simulates distances for units where these are unknown.

First, random numbers of the annual mileage distribution, estimated by dist_mileage, are drawn. Second, the drawn annual distances are converted with respect to the actual operating times (in days) using a linear relationship. See 'Details'.

Usage

mcs_mileage(x, ...)

## S3 method for class 'wt_mcs_mileage_data'
mcs_mileage(x, distribution = c("lognormal", "exponential"), ...)

Arguments

Details

Assumption of linear relationship: Imagine the distance of the vehicle is unknown. A distance of 3500.25 kilometers (km) was drawn from the annual distribution and the known operating time is 200 days (d). So the resulting distance of this vehicle is

$3500.25 km \cdot (\frac{200 d} {365 d}) = 1917.945 km$

Value

A list with class wt_mcs_mileage containing the following elements:

• data : A tibble returned by mcs_mileage_data where two modifications has been made:

  • If the column status exists, the tibble has additional classes wt_mcs_data and wt_reliability_data. Otherwise, the tibble only has the additional class wt_mcs_data (which is not supported by estimate_cdf).

  • The column mileage is renamed to x (to be in accordance with reliability_data) and contains simulated distances for incomplete observations and input distances for the complete observations.

• sim_data : A tibble with column sim_mileage that holds the simulated distances for incomplete cases and 0 for complete cases.

• model_estimation : A list returned by dist_mileage.

See Also

dist_mileage for the determination of a parametric annual mileage distribution and estimate_cdf for the estimation of failure probabilities.

Examples

# MCS data preparation:
mcs_tbl <- mcs_mileage_data(
  field_data,
  mileage = mileage,
  time = dis,
  status = status,
  id = vin
)

# Example 1 - Reproducibility of drawn random numbers:
set.seed(1234)
mcs_distances <- mcs_mileage(
  x = mcs_tbl,
  distribution = "lognormal"
)

# Example 2 - MCS for distances with exponential annual mileage distribution:
mcs_distances_2 <- mcs_mileage(
  x = mcs_tbl,
  distribution = "exponential"
)

# Example 3 - MCS for distances with downstream probability estimation:
## Apply 'estimate_cdf()' to *$data:
prob_estimation <- estimate_cdf(
  x = mcs_distances$data,
  methods = "kaplan"
)

## Apply 'plot_prob()':
plot_prob_estimation <- plot_prob(prob_estimation)

Description

Create consistent mcs_mileage_data based on an existing data.frame (preferred) or on multiple equal length vectors

Usage

mcs_mileage_data(
  data = NULL,
  mileage,
  time,
  status = NULL,
  id = NULL,
  .keep_all = FALSE
)

Arguments

Value

A tibble with class wt_mcs_mileage_data that is formed for the downstream Monte Carlo method mcs_mileage. It contains the following columns (if .keep_all = FALSE):

• mileage : Input mileages.

• time : Input operating times.

• status (optional) :

  • If is.null(status) column status does not exist.

  • If status is provided the column contains the entered binary data (0 or 1).

• id : Identification for every unit.

If .keep_all = TRUE, the remaining columns of data are also preserved.

The attribute mcs_characteristic is set to "mileage".

See Also

dist_mileage for the determination of a parametric annual mileage distribution and mcs_mileage for the Monte Carlo method with respect to unknown distances.

Examples

# Example 1 -  Based on an existing data.frame/tibble and column names:
mcs_tbl <- mcs_mileage_data(
  data = field_data,
  mileage = mileage,
  time = dis,
  status = status
)

# Example 2 - Based on an existing data.frame/tibble and column positions:
mcs_tbl_2 <- mcs_mileage_data(
  data = field_data,
  mileage = 3,
  time = 2,
  id = 1
)

# Example 3 - Keep all variables of the tibble/data.frame entered to argument data:
mcs_tbl_3 <- mcs_mileage_data(
  data = field_data,
  mileage = mileage,
  time = dis,
  status = status,
  id = vin,
  .keep_all = TRUE
)

# Example 4 - Based on vectors:
mcs_tbl_4 <- mcs_mileage_data(
  mileage = field_data$mileage,
  time = field_data$dis,
  status = field_data$status,
  id = field_data$vin
)

Description

This function simulates distances for units where these are unknown.

First, random numbers of the annual mileage distribution, estimated by dist_mileage, are drawn. Second, the drawn annual distances are converted with respect to the actual operating times (in days) using a linear relationship. See 'Details'.

Usage

## Default S3 method:
mcs_mileage(
  x,
  time,
  status = NULL,
  id = paste0("ID", seq_len(length(time))),
  distribution = c("lognormal", "exponential"),
  ...
)

Arguments

Details

Assumption of linear relationship: Imagine the distance of the vehicle is unknown. A distance of 3500.25 kilometers (km) was drawn from the annual distribution and the known operating time is 200 days (d). So the resulting distance of this vehicle is

$3500.25 km \cdot (\frac{200 d} {365 d}) = 1917.945 km$

Value

A list with class wt_mcs_mileage containing the following elements:

• data : A tibble returned by mcs_mileage_data where two modifications has been made:

  • If the column status exists, the tibble has additional classes wt_mcs_data and wt_reliability_data. Otherwise, the tibble only has the additional class wt_mcs_data (which is not supported by estimate_cdf).

  • The column mileage is renamed to x (to be in accordance with reliability_data) and contains simulated distances for incomplete observations and input distances for the complete observations.

• sim_data : A tibble with column sim_mileage that holds the simulated distances for incomplete cases and 0 for complete cases.

• model_estimation : A list returned by dist_mileage.

See Also

dist_mileage for the determination of a parametric annual mileage distribution and estimate_cdf for the estimation of failure probabilities.

Examples

# Example 1 - Reproducibility of drawn random numbers:
set.seed(1234)
mcs_distances <- mcs_mileage(
  x = field_data$mileage,
  time = field_data$dis,
  status = field_data$status,
  id = field_data$vin,
  distribution = "lognormal"
)

# Example 2 - MCS for distances with exponential annual mileage distribution:
mcs_distances_2 <- mcs_mileage(
  x = field_data$mileage,
  time = field_data$dis,
  status = field_data$status,
  id = field_data$vin,
  distribution = "exponential"
)

# Example 3 - MCS for distances with downstream probability estimation:
## Apply 'estimate_cdf()' to *$data:
prob_estimation <- estimate_cdf(
  x = mcs_distances$data,
  methods = "kaplan"
)

## Apply 'plot_prob()':
plot_prob_estimation <- plot_prob(prob_estimation)

Description

This method applies the expectation-maximization (EM) algorithm to estimate the parameters of a univariate Weibull mixture model. See 'Details'.

Usage

mixmod_em(x, ...)

## S3 method for class 'wt_reliability_data'
mixmod_em(
  x,
  distribution = "weibull",
  conf_level = 0.95,
  k = 2,
  method = "EM",
  n_iter = 100L,
  conv_limit = 1e-06,
  diff_loglik = 0.01,
  ...
)

Arguments

Details

The EM algorithm is an iterative algorithm for which starting values must be defined. Starting values can be provided for the unknown parameter vector as well as for the posterior probabilities. This implementation employs initial values for the posterior probabilities. These are assigned randomly by using the Dirichlet distribution, the conjugate prior of a multinomial distribution (see Mr. Gelissen's blog post listed under references).

M-Step : On the basis of the initial posterior probabilities, the parameter vector is estimated with Newton-Raphson.

E-Step : The actual estimated parameter vector is used to perform an update of the posterior probabilities.

This procedure is repeated until the complete log-likelihood has converged.

Value

A list with classes wt_model and wt_mixmod_em. The length of the list depends on the number of specified subgroups k. The first k lists contain information provided by ml_estimation. The values of logL, aic and bic are the results of a weighted log-likelihood, where the weights are the posterior probabilities determined by the algorithm. The last list summarizes further results of the EM algorithm and is therefore called em_results. It contains the following elements:

• a_priori : A vector with estimated prior probabilities.

• a_posteriori : A matrix with estimated posterior probabilities.

• groups : Numeric vector specifying the group membership of every observation.

• logL : The value of the complete log-likelihood.

• aic : Akaike Information Criterion.

• bic : Bayesian Information Criterion.

References

• Doganaksoy, N.; Hahn, G.; Meeker, W. Q., Reliability Analysis by Failure Mode, Quality Progress, 35(6), 47-52, 2002

Examples

# Reliability data preparation:
## Data for mixture model:
data_mix <- reliability_data(
  voltage,
  x = hours,
  status = status
)

# Example 1 - EM algorithm with k = 2:
mix_mod_em <- mixmod_em(
  x = data_mix,
  conf_level = 0.95,
  k = 2,
  n_iter = 150
)

# Example 2 - Maximum likelihood is applied when k = 1:
mix_mod_em_2 <- mixmod_em(
  x = data_mix,
  conf_level = 0.95,
  k = 1,
  n_iter = 150
)

Description

This method applies the expectation-maximization (EM) algorithm to estimate the parameters of a univariate Weibull mixture model. See 'Details'.

Usage

## Default S3 method:
mixmod_em(
  x,
  status,
  distribution = "weibull",
  conf_level = 0.95,
  k = 2,
  method = "EM",
  n_iter = 100L,
  conv_limit = 1e-06,
  diff_loglik = 0.01,
  ...
)

Arguments

Details

The EM algorithm is an iterative algorithm for which starting values must be defined. Starting values can be provided for the unknown parameter vector as well as for the posterior probabilities. This implementation employs initial values for the posterior probabilities. These are assigned randomly by using the Dirichlet distribution, the conjugate prior of a multinomial distribution (see Mr. Gelissen's blog post listed under references).

M-Step : On the basis of the initial posterior probabilities, the parameter vector is estimated with Newton-Raphson.

E-Step : The actual estimated parameter vector is used to perform an update of the posterior probabilities.

This procedure is repeated until the complete log-likelihood has converged.

Value

A list with classes wt_model and wt_mixmod_em. The length of the list depends on the number of specified subgroups k. The first k lists contain information provided by ml_estimation. The values of logL, aic and bic are the results of a weighted log-likelihood, where the weights are the posterior probabilities determined by the algorithm. The last list summarizes further results of the EM algorithm and is therefore called em_results. It contains the following elements:

• a_priori : A vector with estimated prior probabilities.

• a_posteriori : A matrix with estimated posterior probabilities.

• groups : Numeric vector specifying the group membership of every observation.

• logL : The value of the complete log-likelihood.

• aic : Akaike Information Criterion.

• bic : Bayesian Information Criterion.

References

• Doganaksoy, N.; Hahn, G.; Meeker, W. Q., Reliability Analysis by Failure Mode, Quality Progress, 35(6), 47-52, 2002

See Also

mixmod_em

Examples

# Vectors:
hours <- voltage$hours
status <- voltage$status

# Example 1 - EM algorithm with k = 2:
mix_mod_em <- mixmod_em(
  x = hours,
  status = status,
  distribution = "weibull",
  conf_level = 0.95,
  k = 2,
  n_iter = 150
)

#' # Example 2 - Maximum likelihood is applied when k = 1:
mix_mod_em_2 <- mixmod_em(
  x = hours,
  status = status,
  distribution = "weibull",
  conf_level = 0.95,
  k = 1,
  method = "EM",
  n_iter = 150
)

Description

This function uses piecewise linear regression to divide the data into subgroups. See 'Details'.

Usage

mixmod_regression(x, ...)

## S3 method for class 'wt_cdf_estimation'
mixmod_regression(
  x,
  distribution = c("weibull", "lognormal", "loglogistic"),
  conf_level = 0.95,
  k = 2,
  control = segmented::seg.control(),
  ...
)

Arguments

Details

The segmentation process is based on the lifetime realizations of failed units and their corresponding estimated failure probabilities for which intact items are taken into account. It is performed with the support of segmented.lm.

Segmentation can be done with a specified number of subgroups or in an automated fashion (see argument k). The algorithm tends to overestimate the number of breakpoints when the separation is done automatically (see 'Warning' in segmented.lm).

In the context of reliability analysis it is important that the main types of failures can be identified and analyzed separately. These are

• early failures,

• random failures and

• wear-out failures.

In order to reduce the risk of overestimation as well as being able to consider the main types of failures, a maximum of three subgroups (k = 3) is recommended.

Value

A list with classes wt_model and wt_rank_regression if no breakpoint was detected. See rank_regression.

A list with classes wt_model and wt_mixmod_regression if at least one breakpoint was determined. The length of the list depends on the number of identified subgroups. Each list element contains the information provided by rank_regression. In addition, the returned tibble data of each list element only retains information on the failed units and has two more columns:

• q : Quantiles of the standard distribution calculated from column prob.

• group : Membership to the respective segment.

If more than one method was specified in estimate_cdf, the resulting output is a list with classes wt_model and wt_mixmod_regression_list where each list element has classes wt_model and wt_mixmod_regression.

References

Doganaksoy, N.; Hahn, G.; Meeker, W. Q., Reliability Analysis by Failure Mode, Quality Progress, 35(6), 47-52, 2002

Examples

# Reliability data preparation:
## Data for mixture model:
data_mix <- reliability_data(
  voltage,
  x = hours,
  status = status
)

## Data for simple unimodal distribution:
data <- reliability_data(
  shock,
  x = distance,
  status = status
)

# Probability estimation with one method:
prob_mix <- estimate_cdf(
  data_mix,
  methods = "johnson"
)

prob <- estimate_cdf(
  data,
  methods = "johnson"
)

# Probability estimation for multiple methods:
prob_mix_mult <- estimate_cdf(
  data_mix,
  methods = c("johnson", "kaplan", "nelson")
)

# Example 1 - Mixture identification using k = 2 two-parametric Weibull models:
mix_mod_weibull <- mixmod_regression(
  x = prob_mix,
  distribution = "weibull",
  conf_level = 0.99,
  k = 2
)

# Example 2 - Mixture identification using k = 3 two-parametric lognormal models:
mix_mod_lognorm <- mixmod_regression(
  x = prob_mix,
  distribution = "lognormal",
  k = 3
)

# Example 3 - Mixture identification for multiple methods specified in estimate_cdf:
mix_mod_mult <- mixmod_regression(
  x = prob_mix_mult,
  distribution = "loglogistic"
)

# Example 4 - Mixture identification using control argument:
mix_mod_control <- mixmod_regression(
  x = prob_mix,
  distribution = "weibull",
  control = segmented::seg.control(display = TRUE)
)

# Example 5 - Mixture identification performs rank_regression for k = 1:
mod <- mixmod_regression(
  x = prob,
  distribution = "weibull",
  k = 1
)

Description

This function uses piecewise linear regression to divide the data into subgroups. See 'Details'.

Usage

## Default S3 method:
mixmod_regression(
  x,
  y,
  status,
  distribution = c("weibull", "lognormal", "loglogistic"),
  conf_level = 0.95,
  k = 2,
  control = segmented::seg.control(),
  ...
)

Arguments

Details

The segmentation process is based on the lifetime realizations of failed units and their corresponding estimated failure probabilities for which intact items are taken into account. It is performed with the support of segmented.lm.

Segmentation can be done with a specified number of subgroups or in an automated fashion (see argument k). The algorithm tends to overestimate the number of breakpoints when the separation is done automatically (see 'Warning' in segmented.lm).

In the context of reliability analysis it is important that the main types of failures can be identified and analyzed separately. These are

• early failures,

• random failures and

• wear-out failures.

In order to reduce the risk of overestimation as well as being able to consider the main types of failures, a maximum of three subgroups (k = 3) is recommended.

Value

A list with classes wt_model and wt_rank_regression if no breakpoint was detected. See rank_regression. The returned tibble data is of class wt_cdf_estimation and contains the dummy columns cdf_estimation_method and id. The former is filled with NA_character, due to internal usage and the latter is filled with "XXXXXX" to point out that unit identification is not possible when using the vector-based approach.

A list with classes wt_model and wt_mixmod_regression if at least one breakpoint was determined. The length of the list depends on the number of identified subgroups. Each list contains the information provided by rank_regression. The returned tibble data of each list element only retains information on the failed units and has modified and additional columns:

• id : Modified id, overwritten with "XXXXXX" to point out that unit identification is not possible when using the vector-based approach.

• cdf_estimation_method : A character that is always NA_character. Only needed for internal use.

• q : Quantiles of the standard distribution calculated from column prob.

• group : Membership to the respective segment.

References

Doganaksoy, N.; Hahn, G.; Meeker, W. Q., Reliability Analysis by Failure Mode, Quality Progress, 35(6), 47-52, 2002

See Also

mixmod_regression

Examples

# Vectors:
## Data for mixture model:
hours <- voltage$hours
status <- voltage$status

## Data for simple unimodal distribution:
distance <- shock$distance
status_2 <- shock$status

# Probability estimation with one method:
prob_mix <- estimate_cdf(
  x = hours,
  status = status,
  method = "johnson"
)

prob <- estimate_cdf(
  x = distance,
  status = status_2,
  method = "johnson"
)

 # Example 1 - Mixture identification using k = 2 two-parametric Weibull models:
mix_mod_weibull <- mixmod_regression(
   x = prob_mix$x,
   y = prob_mix$prob,
   status = prob_mix$status,
   distribution = "weibull",
   conf_level = 0.99,
   k = 2
)

# Example 2 - Mixture identification using k = 3 two-parametric lognormal models:
mix_mod_lognorm <- mixmod_regression(
   x = prob_mix$x,
   y = prob_mix$prob,
   status = prob_mix$status,
   distribution = "lognormal",
   k = 3
)

# Example 3 - Mixture identification using control argument:
mix_mod_control <- mixmod_regression(
  x = prob_mix$x,
  y = prob_mix$prob,
  status = prob_mix$status,
  distribution = "weibull",
  k = 2,
  control = segmented::seg.control(display = TRUE)
)

# Example 4 - Mixture identification performs rank_regression for k = 1:
mod <- mixmod_regression(
  x = prob$x,
  y = prob$prob,
  status = prob$status,
  distribution = "weibull",
  k = 1
)

Description

This function estimates the parameters of a parametric lifetime distribution for complete and (multiple) right-censored data. The parameters are determined in the frequently used (log-)location-scale parameterization.

For the Weibull, estimates are additionally transformed such that they are in line with the parameterization provided by the stats package (see Weibull).

Usage

ml_estimation(x, ...)

## S3 method for class 'wt_reliability_data'
ml_estimation(
  x,
  distribution = c("weibull", "lognormal", "loglogistic", "sev", "normal", "logistic",
    "weibull3", "lognormal3", "loglogistic3", "exponential", "exponential2"),
  wts = rep(1, nrow(x)),
  conf_level = 0.95,
  start_dist_params = NULL,
  control = list(),
  ...
)

Arguments