R/ddcomp.R

In musto101/wilcox_R: What the package does (short line)

ddcomp <-
function(x, steps = 5)
{
# Makes 4 DD plots using the FMCD and MBA estimators.
# Click left mouse button to identify points.
# Click right mouse button to end the function.
# Unix systems turn on graphics device eg enter
#  command "X11()" or "motif()" before using.
# R users need to type "library(lqs)" before using.
	p <- dim(x)[2]
	par(mfrow = c(2, 2))
	center <- apply(x, 2, mean)
	cov <- var(x)
	md2 <- mahalanobis(x, center, cov)	
	# MD is the classical and RD the robust distance
	MD <- sqrt(md2)	#DGK start
	md2 <- mahalanobis(x, center, cov)
	medd2 <- median(md2)	## get the start
	mns <- center
	covs <- cov	## concentrate
	for(i in 1:steps) {
		md2 <- mahalanobis(x, mns, covs)
		medd2 <- median(md2)
		mns <- apply(x[md2 <= medd2,  ], 2, mean)
		covs <- var(x[md2 <= medd2,  ])
	}
	rd2 <- mahalanobis(x, mns, covs)
	rd <- sqrt(rd2)	#Scale the RD so the plot follows the 0-1 line
#if the data is multivariate normal.
	const <- sqrt(qchisq(0.5, p))/median(rd)
	RDdgk <- const * rd
	plot(MD, RDdgk)
	abline(0, 1)
	identify(MD, RDdgk)
	title("DGK DD Plot")	#MBA
	out <- covmba(x)
	center <- out$center
	cov <- out$cov
	rd2 <- mahalanobis(x, center, cov)
	rd <- sqrt(rd2)	#Scale the RD so the plot follows the identity line
#if the data is multivariate normal.
	const <- sqrt(qchisq(0.5, p))/median(rd)
	RDm <- const * rd
	plot(MD, RDm)
	abline(0, 1)
	identify(MD, RDm)
	title("MBA DD Plot")	#FMCD
	out <- cov.mcd(x)
	center <- out$center
	cov <- out$cov
	rd2 <- mahalanobis(x, center, cov)
	rd <- sqrt(rd2)	#Scale the RD so the plot follows the 0-1 line
#if the data is multivariate normal.
	const <- sqrt(qchisq(0.5, p))/median(rd)
	RDf <- const * rd
	plot(MD, RDf)
	abline(0, 1)
	identify(MD, RDf)
	title("FMCD DD Plot")	#Median Ball start
	covv <- diag(p)
	med <- apply(x, 2, median)
	md2 <- mahalanobis(x, center = med, covv)
	medd2 <- median(md2)	## get the start
	mns <- apply(x[md2 <= medd2,  ], 2, mean)
	covs <- var(x[md2 <= medd2,  ])	## concentrate
	for(i in 1:steps) {
		md2 <- mahalanobis(x, mns, covs)
		medd2 <- median(md2)
		mns <- apply(x[md2 <= medd2,  ], 2, mean)
		covs <- var(x[md2 <= medd2,  ])
	}
	rd2 <- mahalanobis(x, mns, covs)
	rd <- sqrt(rd2)	#Scale the RD so the plot follows the 0-1 line
#if the data is multivariate normal.
	const <- sqrt(qchisq(0.5, p))/median(rd)
	RDmb <- const * rd
	plot(MD, RDmb)
	abline(0, 1)
	identify(MD, RDmb)
	title("Med Ball DD Plot")
}