R/robustRegH.R

In robustreg: Robust Regression Functions

robustRegH<-function(formula,data,tune=1.345,m=TRUE,max.it=1000,tol=1e-5,anova.table=FALSE){

#psiHuber<-function(r,c){
#middle<-abs(r)<=c
#high<- r>c
#low<-r<(-c) 
#h<-middle*r + high*c + low*-c
#return(h)}

bi<-FALSE
if(m==FALSE){bi<-TRUE}

modelFrame=model.frame(formula,data)
X=model.matrix(formula,data)
y=model.extract(modelFrame,"response")

model=lm.fit(X,y)
b=model$coefficients
n<-length(y)
p<-length(b)
if(bi){
 tune<-(tune*sqrt(2*p*n))/(n-2*p)
 hii<-lm.influence(model)$hat
 pi<-(1-hii)/sqrt(hii)}

convergence<-FALSE
for(i in 1:max.it){
 b_temp<-b
 r<-y-fit_rcpp(X,b) #replaced r<-y-X%*%b
 s<-mad_rcpp(r) #replaced s<-median(abs(r-median(r)))/.6745
 
 #rstar<-ifelse(m,r/s,r/(s*pi))
 if(m){rstar<-(r/s)}else{rstar<-r/(s*pi)}


 psiH<-psiHuber_rcpp(rstar,tune)
 w<-psiH/rstar
 b<-lm.wfit(x=X,y=y,w=w[,1])$coefficients

 
 if(i>4){
  if(sum(abs(b-b_temp))<tol){
   cat("\nRobust Regression with Huber Function\n")
   cat("Convergence achieved after:",i,"iterations\n")
   convergence<-TRUE
   break}}
}

#if(convergence==FALSE){b<-NULL;w<-NULL}

#MSE Calc
if(convergence){
 ybarw<-sum(y*w)/sum(w)
 ytild=fit_rcpp(X,b) #replaced ytild<-X%*%b
 sserr<-sum(w*(y-ytild)^2)
 dfr<-length(b)-1 
 dferr<-length(y)-dfr-1
 mse<-sserr/dferr}
else{
 b<-NULL;w<-NULL;mse<-NULL}


if(convergence && anova.table){

derivPsiHuber<-function(r,c){
true<-abs(r)<=c
false<-(r<c*-1 || r>c)
dph<-true*1 +false*0
return(dph)
}
	
r3<-function(x){return(round(x,3))}
r2<-function(x){return(round(x,2))}

 ssreg<-sum(w*(ytild-ybarw)^2)
 sstot<-sum(w*(y-ybarw)^2)
 dftot<-length(y)-1
 msr<-ssreg/dfr
 sbsq<-(s^2*(n^2/(n-length(b)))*sum(psiH^2))/sum(derivPsiHuber(rstar,tune))^2
 F<-msr/sbsq
 W<-Diagonal(x=w[1:length(w),]) 
 c<-diag(x=solve(a=t(X)%*%W%*%X))
 sec<-sqrt(sbsq*c)
 t<-b/sec
 pval<-sapply(t, function(t) if(t>0){2*(1-pt(q=t,df=n-p))}else{2*(1-pt(q=t,df=n-p,lower.tail = FALSE))})

 cat("source  ","\t","SS","\t","\t","df","\t","MS","\t","\t","F","\n")
 cat("model   ","\t",r2(ssreg),"\t",dfr,"\t",r2(msr),"\t",r2(F),"\n")
 cat("error   ","\t",r2(sserr),"\t",dferr,"\n")
 cat("tot     ","\t",r2(sstot),"\t",dftot,"\n")
 cat("rsquared","\t",r3(ssreg/sstot),"\n")
 cat("mse     ","\t",mse,"\n\n")
 cat("Coefficients:\n")
 estimates<-cbind(b,sec,r2(t),round(pval,5))
 colnames(estimates)<-c("estimate","std error","t value","p value") 
 print(estimates)
 }

 object=list("coefficients"=b,"weights"=w, "mse"=mse)
 return(invisible(object)) 
}