Processing math: 72%

CUSUM chart based on linear regression model

The following is an example of an application to a CUSUM chart based on a linear regression model.

Assume we have n past in-control data (Yn,Xn),,(Y1,X1), where Yi is a response variable and Xi is a corresponding vector of covariates. The parameters β of a linear model EY=Xβ are estimated using e.g. the lm function. The corresponding risk adjusted CUSUM chart to detect a shift of Δ>0 in the mean of the response for new observations (Y1,X1),,(Yn,Xn) is then defined by S0=0,St=max

The following generates a data set of past observations (replace this with your observed past data) from the model \mbox{E}Y=2+x_1+x_2+x_3 with standard normal noise and distribution of the covariate values as specified below.

n <- 1000
Xlinreg <- data.frame(x1= rbinom(n,1,0.4), x2= runif(n,0,1), x3= rnorm(n))
Xlinreg$y <- 2 + Xlinreg$x1 + Xlinreg$x2 + Xlinreg$x3 + rnorm(n)

Next, we initialise the chart and compute the estimate needed for running the chart - in this case \hat \beta.

library(spcadjust)
chartlinreg <- new("SPCCUSUM",model=SPCModellm(Delta=1,formula="y~x1+x2+x3"))
xihat <- xiofdata(chartlinreg,Xlinreg)
xihat
## 
## Call:
## lm(formula = formula, data = P)
## 
## Coefficients:
## (Intercept)           x1           x2           x3  
##       2.003        1.091        1.003        1.016

Calibrating the chart to a given average run length (ARL)

We now compute a threshold that with roughly 90% probability results in an average run length of at least 100 in control. In this regression model this is computed by non-parametric bootstrapping. The number of bootstrap replications (the argument nrep) shoud be increased for real applications.

cal <- SPCproperty(data=Xlinreg,
             nrep=100,
             property="calARL",chart=chartlinreg,params=list(target=100),quiet=TRUE)
cal
## 90 % CI: A threshold of 3.06 gives an in-control ARL of at least
##   100. 
## Unadjusted result:  2.817 
## Based on  100 bootstrap repetitions.

Run the chart

Next, we run the chart with new observations that are in-control.

n <- 100
newXlinreg <- data.frame(x1= rbinom(n,1,0.4), x2= runif(n,0,1), x3=rnorm(n))
newXlinreg$y <- 2 + newXlinreg$x1 + newXlinreg$x2 + newXlinreg$x3 + rnorm(n)
S <- runchart(chartlinreg, newdata=newXlinreg,xi=xihat)

plot of chunk unnamed-chunk-7

In the next example, the chart is run with data that that are out-of-control from time 51 and onwards.

n <- 100
newXlinreg <- data.frame(x1= rbinom(n,1,0.4), x2= runif(n,0,1),x3=rnorm(n))
outind <- c(rep(0,50),rep(1,50))
newXlinreg$y <- 2 + newXlinreg$x1 + newXlinreg$x2 + newXlinreg$x3 + rnorm(n)+outind
S <- runchart(chartlinreg, newdata=newXlinreg,xi=xihat)

plot of chunk unnamed-chunk-9