Markovchart (version 2.1.5)

Markovchart: Cost-efficient X-bar control charts with fixed/random shift size, random repair and random sampling time.

Description

Wrapper for Markov chain-based cost optimal control charts. Includes cost calculation methods for different shift size distributions and optimisation with respect to the average cost and cost standard deviation where the free parameters are the sampling interval (h) and the control limit/critical value (k).

Usage

Markovchart(statdist, h = NULL, k = NULL,
            OPTIM = FALSE, p = 1, constantr = FALSE,
            ooc_rep = 0, cs = NULL, cofun = cofun_default,
            coparams = NULL, crfun = crfun_default, crparams = NULL,
            cf = crparams, vcofun = vcofun_default,
            vcoparams = c(0, 0), vcrfun = vcrfun_default,
            vcrparams = c(0, 0), method = c("L-BFGS-B", "Nelder-Mead",
            "BFGS", "CG", "SANN", "Brent"), parallel_opt = NULL,
            silent = TRUE, ...)

Arguments

statdist

The stationary distribution of the Markov chain. Must be an object of class Markov_stationary, preferably created by Markovstat.

h

The time between samplings. Must be a positive value, can be a numeric vector. For optimisation, this is the initial value. Inherited from statdist if not given.

k

The control limit (critical value). Must be a positive value, can be a numeric vector. For optimisation, this is the initial value. Only one sided shifts are allowed, thus there is only one control limit. Inherited from statdist if not given.

OPTIM

Logical. Should the resulting G-value (weighted average of the expected cost and cost standard deviation) be optimised by finding the adequate value of h and k.

p

The weight of the cost expectation in the calculation of the G-value; should be between 0 and 1.

constantr

Logical. Should the repair cost be assumed to constantly occur over time (TRUE) or assumed to only occur when there is a repair due to an alarm (FALSE, default)? If TRUE, then the repair cost should be given per unit time.

ooc_rep

Numeric value between 0 and 1. The percentage of repair cost ocurring during out-of-control operation. Default is 0. If a value greater than 0 is set, then constantr should be TRUE, but it is not forced.

cs

Sampling cost per sampling.

cofun

A function describing the relationship between the distance from the target value and the resulting out-of-control costs. Default is calculated using a base and a distance-scaling out-of-control parameter. See "Details".

coparams

Numeric vector. Parameters of cofun.

crfun

A function describing the relationship between the distance from the target value and the resulting repair costs. The default function assumes a linear relationship between the repair cost and the distance, and uses a base and a distance-scaling repair cost parameter. See "Details".

crparams

Numeric vector. Parameters of crfun.

cf

Numeric. The false alarm cost. Only relevant when statdist was created using shiftfun="deg".

vcofun

A function describing the relationship between the distance from the target value and the resulting out-of-control cost variance. For the default function see "Details".

vcoparams

Numeric vector. Parameters of vcofun.

vcrfun

A function describing the relationship between the distance from the target value and the resulting repair cost variance. For the default function see "Details".

vcrparams

Numeric vector. Parameters of vcrfun.

method

Method used for optimisation. Same as the argument of optim, but the default here is "L-BFGS-B", because it turned out to be more robust in testing.

parallel_opt

A list of parallel options. See e.g. the argument parallel in the documentation of optimParallel. Can be left empty, in this case the number of cores (threads) is automatically detected and all but one is used. (Single-core computers simply use one core.)

silent

Should the call be returned? Default is FALSE.

…

Further arguments to be passed down to optimParallel.

Value

The value depends on the parameters:

If either h or k have length greater than 1, then the G-value (weighted average of average cost and cost standard deviation) is calculated for all given values without optimisation. The value of the function in this case is a data frame of class codeMarkov_grid with length(h)*length(k) number of rows and three columns for h, k and the G-value.

If h and k are both of length 1 (they may be inherited from statdist), then the value of the function is a Markov_chart object, which is a list of length 4, detailing the properties of the control chart setup.

Results

Vector of G-value, expected cost, cost standard deviation and further process moments. Note that these further moments only take into account the process variation (i.e. the standard deviation of the process itself), while the "Total cost std. dev." takes into account all sources of variance (e.g. the different costs that can occur due to being out-of-control). The "Total cost std. dev." is only relevant and calculated for non-degenerate distributions.

Subcosts

Vector of sub-costs that are parts of the total expected cost.

Parameters

A vector that contains the time between samplings (h) and critical value (k) which was used in the control chart setup.

Stationary_distribution

The stationary distribution of the Markov chain. Further information about the stationary distribution can be calculated using the Markovstat function.

Details

The constantr parameter is used for different repair assumptions. In traditional control chart theory, repair cost only occurs in case of an alarm signal. This is represented by constantr=FALSE, which is the default. In this case the repair is just a momentary cost, occurring at the time of the sampling. However this model is inappropriate in several cases in healthcare. For example there are chronic diseases that require constant medication (repair in the sense of the model). In this approach (constantr=TRUE) the repair cost still depends on the state of the process during sampling, but occurs even if there is no alarm and is divided by h to represent the constant repair through the whole sampling interval. Thus the repair cost should be given in a way which corresponds to the model used.

The default cofun calculates the out-of-control (OOC) cost using a base and a distance-scaling OOC parameter:

$$c_{o} = c_{ob} + c_{os} A^2(v),$$

where \(c_{o}\) is the total OOC cost, \( c_{ob}\) is the base OOC cost (even without shift), \(c_{os}\) is the shift-scaling cost and \(A^2(v)\) is the squared distance from the target value. This latter part is defined like this because a Taguchi-type loss function is used. This \(A^2(v)\) incorporates the distances (the base of the losses) incurred not just at the time of the sampling, but also between samplings (hence it dependens on h). Even if the user defines a custom cost function for the OOC cost, this \(A^2(v)\) term must be included, as a closed form solution has been developed for the calculation of the squared distances in case of exponential shifts, considerably decreasing run times. Thus the arguments of the OOC cost function should look like this: function(\(A^2(v)\), other parameters contained in a vector). \(A^2(v)\) is fed to the cost function as a vector, thus the function should vectorised with respect to this argument. The default function looks like this:

cofun_default	<-	function(sqmudist,coparams)
{
  sqmudist=sqmudist
  cob=coparams[1]
  cos=coparams[2]
  co		<-	cob + cos*sqmudist
  return(co)
}

The default vcofun also uses a Taguchi-type loss function and has identical parts and requirements as cofun. The final standard deviation itself is calculated using the law of total variance. The default vcofun is:

vcofun_default	<-	function(sqmudist,vcoparams)
{
  sqmudist=sqmudist
  vcob=vcoparams[1]
  vcos=vcoparams[2]
  vco		<-	vcob + vcos*sqmudist
  return(vco)
}

The defaults for the repair cost and cost variance are simple linear functions. For crfun it is

$$c_{r} = c_{rb} + c_{rs} v,$$

where the notation are the same as before and "r" denotes repair. A custom function can be defined more freely here, but the first argument should be \(v\) and the second a vector of further parameters.

The default function are:

crfun_default	<-	function(mudist,crparams)
{
  mudist=mudist
  crb=crparams[1]
  crs=crparams[2]
  cr		<-	crb + crs*mudist
  return(cr)
}

vcrfun_default	<-	function(mudist,vcrparams)
{
  mudist=mudist
  vcrb=vcrparams[1]
  vcrs=vcrparams[2]
  vcr		<-	vcrb + vcrs*mudist;
  return(vcr)
}

References

Zempleni A, Veber M, Duarte B and Saraiva P. (2004) Control charts: a cost-optimization approach for processes with random shifts. Applied Stochastic Models in Business and Industry, 20(3), 185-200.

Dobi B and Zempleni A. (2019) Markov chain-based cost-optimal control charts for health care data. Quality and Reliability Engineering International, 35(5), 1379-1395.

Dobi B and Zempleni A. (2019) Markov chain-based cost-optimal control charts with different shift size distributions. Annales Univ. Sci. Budapest., Sect. Comp., 49, 129-146.

Examples

# NOT RUN {
#Defining parallel_opt parallel settings.
#parallel_opt can also be left empty to be defined automatically by the function.
require(parallel)
num_workers <-  min(c(detectCores(),2))
parall	    <-  list(cl=makeCluster(num_workers), forward=FALSE, loginfo=TRUE)


#Fixed shift size (essentially Duncan's cycle model) - no optimisation.
stat_deg <- Markovstat(shiftfun="deg", h=1, k=1, sigma=1, s=0.2, delta=2.5)
res1     <- Markovchart(statdist=stat_deg, cs=1, crparams=20, coparams=50)
res1
# }
# NOT RUN {
#Fixed shift size (essentially Duncan's cycle model) - with optimisation.
res2 <-  Markovchart(statdist=stat_deg, OPTIM=TRUE, cs=1, crparams=20, coparams=50,
                    lower = c(0.01,0.01), upper = c(5,5),
                    parallel_opt=parall)
res2

#Exponential shift - no optimisation - default cost functions.
stat_exp <- Markovstat(shiftfun="exp", h=0.5, k=2, sigma=1, s=0.2, delta=2,
                       RanRep=FALSE, Vd=30, V=18)
res3     <- Markovchart(stat_exp, p=0.9, cs=1, coparams=c(10,3), crparams=c(1,2))
res3
# }
# NOT RUN {
#Exponential shift - with optimisation - default cost functions.
stat_exp2 <- Markovstat(shiftfun="exp", h=1, k=1, sigma=1, s=0.2, delta=2,
                        RanRep=TRUE, alpha=1, beta=3, Vd=30, V=18)
parall    <- list(cl=makeCluster(num_workers), forward=FALSE, loginfo=TRUE)
res4      <- Markovchart(statdist=stat_exp2, OPTIM=TRUE, p=0.9, cs=1,
                         coparams=c(10,3), crparams=c(1,2), vcoparams=c(8,1.5),
                         vcrparams=c(5,2), parallel_opt=parall)
res4
# }
# NOT RUN {
#Exponential-geometric mixture shift - no optimisation -
#random sampling - custom repair variance function.
stat_expgeo <- Markovstat(shiftfun="exp-geo",h=1.5, k=2, sigma=1,
                          s=0.2, delta=1.2, probmix=0.7, probnbin=0.8,
                          disj=2, RanRep=TRUE, alpha=1, beta=3, RanSam=TRUE,
                          StateDep=TRUE, a=1, b=15, Vd=100, V=8)

vcrfun_new	<-	function(mudist,vcrparams)
{
  mudist=mudist
  vcrb=vcrparams[1]
  vcrs=vcrparams[2]
  vcrs2=vcrparams[3]

  vcr		<-	vcrb + vcrs/(mudist + vcrs2)
  return(vcr)
}

res5 <- Markovchart(statdist=stat_expgeo, p=0.9, cs=1,
                   coparams=c(10,6), crparams=c(20,3),
                   vcoparams=c(10000,100), vcrfun=vcrfun_new,
                   vcrparams=c(50000,-600000,1.5))
res5

#Exponential shift - no optimisation  - vectorised.
parall <- list(cl=makeCluster(num_workers), forward=FALSE, loginfo=TRUE)
Gmtx   <-	Markovchart(statdist=stat_exp2, h=seq(1,10,by=(10-1)/5),
                      k=seq(0.1,5,by=(5-0.1)/5), p=0.9, cs=1,
                      coparams=c(10,3), crparams=c(1,2),
                      vcoparams=c(8,1.5), vcrparams=c(5,2),
                      V=18, parallel_opt=parall)
Gmtx
# }