#' @title Sequential Normal Scores
#' @description Transform a vector \code{X} into SNS using initial observations \code{Y} if available
#' SNS follow the order of \code{X}.
#' @section Comments:
#' If ties occur, average ranks are used.
#' @seealso \code{\link{NS}} for normal scores
#' @inheritParams NS
#' @inheritParams getRL
#' @param X.id vector. The id of the vector \code{X}.
#' @param snsRaw logical. If \code{TRUE} return also the sns for each observation in vector \code{X}.
#' @param omit.id vector. Elements of the vector are the id which are omitted in the analysis.
#' @param auto.omit.alarm logical. Determine if OC signals are added (or not) to reference sample. By default is set to TRUE.
#' @return Multiple output. Select by \code{output$}
#' \itemize{
#' \item \code{coefficients}: list. Three elements: \code{n} the number of observation per group in \code{X}, \code{chart} the selected chart to perform the analysis, and \code{chart.par} the parameters of the selected chart.
#' \item \code{R}: vector. Ranks for the new observations (Monitoring sample).
#' \item \code{X}: vector. New observations (Monitoring sample) to obtain the SNS.
#' \item \code{Z}: vector. SNS of the \code{X} monitoring sample.
#' \item \code{X.id}: vector. The id of each column (variable) of the matrix \code{X}.
#' \item \code{UCL}: vector. Upper control limit for each group in \code{X}.
#' \item \code{LCL}: vector. Lower control limit for each group in \code{X}.
#' \item \code{scoring}: string. Selected score to evaluate SNS.
#' }
#' @export
#' @examples
#' # EXAMPLE CONDITIONAL WITH REFERENCE SAMPLE
#' Y <- c(10, 20, 30, 40, 50, 60, 70, 80, 90, 100)
#' X <- c(30, 35, 45)
#' theta <- 40
#' Ftheta <- 0.5
#' sample.id <- c("a", "b", "c")
#' SNS(X = X, X.id = sample.id, Y = Y, theta = theta, Ftheta = Ftheta)
#'
#' # EXAMPLE CONDITIONAL WITH REFERENCE SAMPLE
#' Y <- c(10, 20, 30, 40, 50, 60, 70, 80, 90, 100)
#' X <- c(30, 35, 45)
#' theta <- 40
#' Ftheta <- 0.5
#' sample.id <- c("a", "b", "c")
#' SNS(X = X, X.id = sample.id, Y = Y, theta = theta, Ftheta = Ftheta)
#'
#' # EXAMPLE UNCONDITIONAL WITH REFERENCE SAMPLE
#' Y <- c(10, 20, 30, 40, 50, 60, 70, 80, 90, 100)
#' X <- c(30, 35, 45)
#' theta <- NULL
#' Ftheta <- NULL
#' sample.id <- c("a", "b", "c")
#' SNS(X = X, X.id = sample.id, Y = Y, theta = theta, Ftheta = Ftheta)
#'
#' # EXAMPLE CONDITIONAL WITHOUT REFERENCE SAMPLE
#' Y <- NULL # c(10,20,30,40,50,60,70,80,90,100)
#' X <- c(30, 35, 45)
#' theta <- 40
#' Ftheta <- 0.5
#' sample.id <- c("a", "b", "c")
#' SNS(X = X, X.id = sample.id, Y = Y, theta = theta, Ftheta = Ftheta)
#'
#' # EXAMPLE UNCONDITIONAL WITHOUT REFERENCE SAMPLE
#' Y <- NULL
#' X <- c(30, 35, 45)
#' theta <- NULL
#' Ftheta <- NULL
#' sample.id <- c("a", "b", "c")
#' SNS(X = X, X.id = sample.id, Y = Y, theta = theta, Ftheta = Ftheta)
SNS <- function(X, X.id, Y = NULL, theta = NULL, Ftheta = NULL,
scoring = "Z", Chi2corrector="None",
alignment = "unadjusted", constant = NULL, absolute = FALSE,
chart="Shewhart", chart.par=c(3),
snsRaw = FALSE, isFixed = FALSE,
omit.id = NULL, auto.omit.alarm = TRUE) {
if (is.null(theta) != is.null(Ftheta)) { # in case one is NULL and not the other
stop("theta or Ftheta missing")
return()
} else if (length(X) != length(X.id)) {
stop("observations (X) have different length of the observations id (X.id)")
return()
}
omit.id.found = NULL
if(!is.null(omit.id)){#check which groups are omitted
ids = unique(X.id)
omit.id.found = which(ids %in% omit.id)
omit.id.missing = omit.id[!(omit.id %in% ids)]
if(!is.null(omit.id.missing) || length(omit.id.missing) > 0 ){
warning("ids to omit not found:", omit.id.missing, "\n")
}
if(auto.omit.alarm){
auto.omit.alarm = FALSE
warning("auto.omit.alarm = FALSE make omit.id = NULL to enable.")
}
if(is.null(omit.id.found) || length(omit.id.found) == 0){
warning("omitted ids not found, OC signals not added to reference sample (auto.omit.alarm = TRUE).")
}
}
# detect the changes in the observation id vector
changes.in.X.id = c(1, as.numeric(X.id[1:(length(X.id) - 1)] != X.id[2:(length(X.id))]))
#change the observation id
Xb.id = cumsum(changes.in.X.id)
#get the different groups of the id
groups = unique(Xb.id)
z = rep(NA, length(groups)) # preallocate memory to initialize the SNS (one for group)
r = rep(NA, length(groups)) # preallocate memory to initialize the ranks (one for group)
if(snsRaw){
Zraw = rep(NA, length(X))
Rraw = rep(NA, length(X))
}
i = 1 # initialize the group index of the observation id vector
Yb = Y
if(!is.null(Yb)){
Yb = Yb[!is.na(Yb)] # initialize reference sample (remove na values)
}
UCL = rep(NA, length(groups))
LCL = rep(NA, length(groups))
switch(chart,
Shewhart = {
k <- chart.par[1]
},
CUSUM = {
k <- chart.par[1]
h <- chart.par[2]
type <- chart.par[3]
Cplus = rep(NA, length(groups))
Cminus = rep(NA, length(groups))
cplus <- 0
cminus <- 0
},
EWMA = {
lambda <- chart.par[1]
L <- chart.par[2]
E = rep(NA, length(groups))
e <- 0
}
)
ucl = 0
if (scoring == "Z-SQ"){
inf = 1000000 #infinite value to better approximation
alpha = 0.005 #confindent interval
vec <- rchisq(inf, length(X[which(Xb.id == groups[i])])) #chi-sq random generator numbers according to the "infinite value"
ucl <- quantile(vec , 1-alpha) #control limit
}
while (i <= length(groups)) { # repeat until the total groups are analized
Xb = X[which(Xb.id == groups[i])] # get the observations to evalute from the positions
n = length(Xb)
ad = SNSchart::dataAlignment(Xb, Yb, alignment = alignment)
Xb = ad$X
Yb = ad$Y
ns = SNSchart::NS(X = Xb, Y = Yb, theta = theta, Ftheta = Ftheta, scoring = scoring, Chi2corrector = Chi2corrector, alignment = alignment, constant = constant) # calculate the normal score
r[i] = mean(ns$R)
if(snsRaw){#save raw data
Zraw[(1+n*(i-1)):(n+n*(i-1))] = ns$Z
Rraw[(1+n*(i-1)):(n+n*(i-1))] = ns$R
}
ns = ns$Z
switch (scoring,
"Z" = {# it is a vector with a subgroup size so it is needed to average them
z[i] = sum(ns)/sqrt(n)
},
"Z-SQ" = {# it is a vector with a subgroup size so it is needed to sum them
z[i] = sum(ns)
}
)
Z = z[i]
# check if the subgroup is in control according to each scheme
# the reference sample is updated
updateSample <- FALSE
switch(chart,
Shewhart = {
if (scoring == "Z"){
ucl = k
}
if (abs(Z) < ucl) updateSample <- TRUE
},
CUSUM = {
switch(type,
"1" = {
cplus <- max(c(0, cplus + Z - k))
},
"2" = {
cminus <- min(c(0, cminus + Z + k))
},
"3" = {
cplus <- max(c(0, cplus + Z - k))
cminus <- min(c(0, cminus + Z + k))
}
)
Cplus[i] <- cplus
Cminus[i] <- cminus
ucl = h
if (cplus < ucl || cminus > -ucl) updateSample <- TRUE
},
EWMA = {
e <- lambda * Z + (1 - lambda) * e
E[i] = e
ucl <- L * sqrt(lambda / (2 - lambda) * (1 - (1 - lambda)^(2 * i)))
if (abs(e) < ucl) updateSample <- TRUE
}
)
UCL[i] = ucl
LCL[i] = -ucl
if (scoring == "Z-SQ"){
LCL[i] = 0
}
if (auto.omit.alarm){
if (updateSample && !isFixed){# if the subgroup is in control (updateSample change to TRUE)
Yb = c(Yb, Xb) # add to reference sample the new observations
}
}else{
if (!(i %in% omit.id.found) && !isFixed){#and if the id of the group is not omitted
Yb = c(Yb, Xb) # add to reference sample the new observations
}
}
i = i + 1 # continue with the next group
}
output = list(
coefficients = list(
n=n,
chart = chart,
chart.par = chart.par
),
R = r,
Z = z,
X.id = X.id,
UCL = UCL,
LCL = LCL,
scoring = scoring
)
if(snsRaw){
output$Zraw = Zraw
output$Rraw = Rraw
}
switch(chart,
CUSUM = {
output$Cplus = Cplus
output$Cminus = Cminus
},
EWMA = {
output$E = E
}
)
class(output)="SNS" # Class definition
return(output) # return the sequential normal score
}
#' @import graphics
#' @export
plot.SNS <- function(x, ...){
#Recommended margins
#par(mar = c(6,6,4,2))
args = list(...)
names.args = names(args)
Z = x$Z
n = x$n
o.id = unique(x$X.id) # original id
chart = coef(x)$chart
chart.par = coef(x)$chart.par
UCL = x$UCL
LCL = x$LCL
Cplus = x$Cplus
Cminus = x$Cminus
E = x$E
scoring = x$scoring
CEX = 0.75
switch(chart,
Shewhart = {
difMaxZ = abs(max(Z) - max(UCL))
difMinZ = abs(min(Z) - min(LCL))
},
CUSUM = {
difMaxZ = abs(max(c(Cplus, Cminus)) - max(UCL))
difMinZ = abs(min(c(Cplus, Cminus)) - min(LCL))
},
EWMA = {
difMaxZ = abs(max(E) - max(UCL))
difMinZ = abs(min(E) - min(LCL))
}
)
ymax = max(UCL)
ymin = min(LCL)
if(difMaxZ > difMinZ){
ymax = ymax + difMaxZ
ymin = ymin - difMaxZ
}else{
ymax = ymax + difMinZ
ymin = ymin - difMinZ
}
switch(chart,
Shewhart = {
ylab = "Z"
if (scoring == "Z-SQ"){
ymin = 0
ylab = expression(Z^2)
}
plot(o.id, Z, pch=19, ylim=c(ymin, ymax), xlab = "Batch",ylab=ylab,cex.lab=CEX*2, cex.axis=CEX*1.25, cex=CEX, col = ifelse((LCL > Z) | (Z > UCL), "red", "black"))
lines(o.id, Z, lt=2, lwd=CEX*1.5)
},
CUSUM = {
type = chart.par[3]
switch(type,
"1" = {
plot(o.id, Cplus, pch=19, ylim=c(ymin, ymax), xlab = "Batch",ylab=expression(C^"+"),cex.lab=CEX*2, cex.axis=CEX*1.25, cex=CEX, col = ifelse(Cplus > UCL, "red", "black"))
lines(o.id, Cplus, lt=2, lwd=CEX*1.5)
},
"2" = {
plot(o.id, Cminus, pch=19, ylim=c(ymin, ymax), xlab = "Batch",ylab=expression(C^"-"),cex.lab=CEX*2, cex.axis=CEX*1.25, cex=CEX, col = ifelse(Cminus < LCL, "red", "black"))
lines(o.id, Cminus, lt=2, lwd=CEX*1.5)
},
"3" = {
plot(o.id, Cplus, pch=19, ylim=c(ymin, ymax), xlab = "Batch",ylab=expression(C^"+"/C^"-"),cex.lab=CEX*2, cex.axis=CEX*1.25, cex=CEX, col = ifelse(Cplus > UCL, "red", "black"))
points(o.id, Cminus, pch=15, cex=CEX, col = ifelse(Cminus < LCL, "red", "black"))
lines(o.id, Cplus, lt=2, lwd=CEX*1.5)
lines(o.id, Cminus, lt=2, lwd=CEX*1.5)
legend("topleft", c(expression(C^"+"), expression(C^"-")), pch=c(19, 15))
}
)
},
EWMA = {
plot(o.id, E, pch=19, ylim=c(ymin, ymax), xlab = "Batch",ylab="E",cex.lab=CEX*2, cex.axis=CEX*1.25, cex=CEX, col = ifelse(E > UCL, "red", "black"))
lines(o.id, E, lt=2, lwd=CEX*1.5)
}
)
change = 1
if(o.id[1] > o.id[length(o.id)]){
change = -1
}
lines(o.id, UCL,lt=4, lwd=CEX*1.5, col = "gray48")
if(sum(LCL) != 0) lines(o.id, LCL,lt=4, lwd=CEX*1.5, col = "gray48")
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.