R/mbalata.R

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

mbalata <-
function(x, y, k=6, nsamp = 7)
{
#gets the median ball fit with 7 centers, med resid crit, 7 ball sizes
        x <- as.matrix(x)
	n <- dim(x)[1]
	q <- dim(x)[2]	
	# q + 1 is number of predictors including intercept
	vals <- c(q + 3 + floor(n/100), q + 3 + floor(n/40), q + 3 +
		floor(n/20), q + 3 + floor(n/10), q + 3 + floor(n/5), q +
		3 + floor(n/3), q + 3 + floor(n/2))
	covv <- diag(q)
	centers <- sample(n, nsamp)
	temp <- lsfit(x, y)
	mbaf <- temp$coef	## get LATA criterion
	res <- temp$residuals
	crit <- k^2*median(res^2)
	cn <- sum(res^2 <= crit)
	absres <- sort(abs(res))
	critf <- sum(absres[1:cn])	##
	for(i in 1:nsamp) {
		md2 <- mahalanobis(x, center = x[centers[i],  ], covv)
		smd2 <- sort(md2)
		for(j in 1:7) {
			temp <- lsfit(x[md2 <= smd2[vals[j]],  ], y[md2 <=
				smd2[vals[j]]])	
	#Use OLS on rows with md2 <= cutoff = smd2[vals[j]]
			res <- y - temp$coef[1] - x %*% temp$coef[-1]	
	## get LATA criterion
			crit <- k^2*median(res^2)
			cn <- sum(res^2 <= crit)
			absres <- sort(abs(res))
			crit <- sum(absres[1:cn])	##
			if(crit < critf) {
				critf <- crit
				mbaf <- temp$coef
			}
		}
	}
	list(coef = mbaf, critf = critf)
}