#  This routine fits a full model and a reduced model (using
#  rrega with method=function(u) wt.huber(u,c=1.5)) and then
#  computes the approximate F statistic of equation (5.6) in 
#  Birkes and Dodge.
#
#  Any subset of variables can be specified to be in the full model
#  and any subset of those can be specified to be in the reduced model.
#  The variables are specified by listing the corresponding
#  columns of x in vful and vred

FMfittest  <- function(x,y,vful=c(1:ncol(x)),vred=c(1))
{
	# x <- regression predictors
	# y <- regression dependent
	# nrow(x) must == dim(y) 
	# you can choose any subset of x as full variables and reduced
	# variables in x as you wish, for example
	# vred <- c(1,2,4) means choosing the colums 1,2,4 of x as your 
	# reduced variables
	rrful  <- rrega(x[,vful],y,method=function(u) wt.huber(u,c=1.5))
  	sighatf    <- 1.483*median(abs(rrful$res))
  	rrared <- rrega(x[,vred],y,method=function(u) wt.huber(u,c=1.5),
                  fix.scale=sighaf)
	resful <- rrful$res
	resred <- rrared$res
        dful <- length(vful)
        dred <- length(vred)
	for(i 1:dred)
	{
		if(min(abs(vful-vred[i]))==0)
		stop("vred is wrong,the variables are not for reduced model!") 
	}
	n <- length(resful)
	if(n != length(resred))
	   stop("residual vectors not the same length")
   
	#  Full model
	mineif  <- ifelse(resful<1.5*sighatf,resful,1.5*sighatf)
	fif     <- ifelse(-1.5*sighatf>mineif,-1.5*sighatf,mineif)
	gif     <- ifelse(0>abs(resful)-1.5*sighatf,0,abs(resful)-1.5*sighatf)
	s1      <- sum(fif*fif)
	strful  <- s1+3*sighatf*sum(gif)
	# Reduced model
	mineir  <- ifelse(resred<1.5*sighatf,resred,1.5*sighatf)
	fir     <- ifelse(-1.5*sighatf>mineir,-1.5*sighatf,mineir)
	gir     <- ifelse(0>abs(resred)-1.5*sighatf,
           0,abs(resred)-1.5*sighatf)
	strred  <- sum(fir*fir)+3*sighatf*sum(gir)

	if(strred <= strful)
 	stop("sum of transformed residuals for the reduced model is not greater 
		than STR for full model")
	mzero  <- sum(gif[gif==0]+1)
	lambda <- n*s1/(mzero*(n-dful-1))
	FM <- (strred-strful)/((dful-dred)*lambda)
	return(FM)
}