robust<-function(y,x){ yall=as.matrix(cbind(y,x)) newdata=na.omit(yall) nk=dim(newdata) nx=dim(as.matrix(x)) n=nk[1] k=nx[2] ind=k+1 ynew=newdata[,1] y=ynew xnew=as.matrix(newdata[,2:ind]) x=xnew output<-lm(y~x) X=as.matrix(cbind(1,x)) resids=output$resid Xp=t(X) I=diag(n) XX=Xp%*%X Xinv=solve(XX) u2=resids^2 Omega=matrix(0,ind,ind) W=t(t(apply(as.matrix(u2),1,rep,ind))*X)%*%X for(i in 1:n){Omega1=Omega+as.matrix(u2[i]*as.matrix(Xp[,i])%*%as.matrix(t(X[i,]))) Omega=Omega1} varcovr<-(n/(n-k))*Xinv%*%Omega%*%Xinv varcovr<-(n/(n-k))*Xinv%*%W%*%Xinv return(list(output=summary(output),RobustSE=diag(varcovr)^(1/2))) }