wilcox_R: What the package does (short line)

kerreg <-
function(x,y,pyhat=FALSE,pts=NA,plotit=TRUE,theta=50,phi=25,expand=.5,
scale=FALSE,zscale=FALSE,eout=FALSE,xout=FALSE,outfun=out,np=100,xlab="X",ylab="Y",zlab="Z",
varfun=pbvar,e.pow=TRUE,pr=TRUE,ticktype="simple",pch='.',...){
#
# Compute local weighted regression with Epanechnikov kernel
#
# See Fan, Annals of Statistics, 1993, 21, 196-217.
# cf. Bjerve and Doksum, Annals of Statistics, 1993, 21, 890-902
#
# With a single predictor, this function calls locreg
# See locreg for information about np and plotting
#
library(akima)
x<-as.matrix(x)
xx<-cbind(x,y)
xx<-elimna(xx)
x<-xx[,1:ncol(x)]
x<-as.matrix(x)
y<-xx[,ncol(x)+1]
d<-ncol(x)
np1<-d+1
m<-elimna(cbind(x,y))
if(xout && eout)stop("Can't have eout=xout=T")
if(eout){
flag<-outfun(m,plotit=FALSE,...)$keep
m<-m[flag,]
}
if(xout){
flag<-outfun(x,plotit=FALSE,...)$keep
m<-m[flag,]
}
if(zscale){
for(j in 1:np1){
m[,j]<-(m[,j]-median(m[,j]))/mad(m[,j])
}}
x<-m[,1:d]
x<-as.matrix(x)
y<-m[,np1]
n<-nrow(x)
if(d>1){
xrem<-x
pi<-gamma(.5)^2
cd<-c(2,pi)
if(d==2)A<-1.77
if(d==3)A<-2.78
if(d>2){
for(j in 3:d)cd[j]<-2*pi*cd[j-2]/j  # p. 76
}
if(d>3)A<-(8*d*(d+2)*(d+4)*(2*sqrt(pi))^d)/((2*d+1)*cd[d])  # p. 87
hval<-A*(1/n)^(1/(d+4))  # p. 86
for(j in 1:d){
sig<-sqrt(var(x[,j]))
temp<-idealf(x[,j])
iqr<-(temp$qu-temp$ql)/1.34
A<-min(c(sig,iqr))
x[,j]<-x[,j]/A
}
xx<-cbind(rep(1,nrow(x)),x)
yhat<-NA
for(j in 1:n){
yhat[j]<-NA
temp1<-t(t(x)-x[j,])/(hval)
temp1<-temp1^2
temp1<-apply(temp1,1,FUN="sum")
temp<-.5*(d+2)*(1-temp1)/cd[d]
epan<-ifelse(temp1<1,temp,0) # Epanechnikov kernel, p. 76
chkit<-sum(epan!=0)
if(chkit >= np1){
vals<-lsfit(x,y,wt=epan)$coef
yhat[j]<-xx[j,]%*%vals
}}
if(plotit  && d==2){
if(pr){
if(!scale){
print("scale=F is specified")
print("If there is dependence, might use scale=T")
}}
m<-elimna(cbind(xrem,yhat))
xrem<-m[,1:d]
yhat<-m[,np1]
fitr<-yhat
iout<-c(1:length(fitr))
nm1<-length(fitr)-1
for(i in 1:nm1){
ip1<-i+1
for(k in ip1:length(fitr))if(sum(xrem[i,]==xrem[k,])==2)iout[k]<-0
}
fitr<-fitr[iout>=1]
mkeep<-xrem[iout>=1,]
fit<-interp(mkeep[,1],mkeep[,2],fitr)
persp(fit,theta=theta,phi=phi,expand=expand,xlab=xlab,ylab=ylab,zlab=zlab,
scale=scale,ticktype=ticktype)
}}
if(d==1){
yhat<-locreg(x[,1],y,pyhat=TRUE,np=np,plotit=plotit,pts=pts,
xlab=xlab,ylab=ylab,pch=pch)
yhat2<-locreg(x[,1],y,pyhat=TRUE,np=0,plotit=FALSE)
}
if(d>1)yhat2<-yhat
m<-NULL
#E.pow<-varfun(yhat2[!is.na(yhat2)])/varfun(y)
# Estimate of explanatory power performs poorly.
if(pyhat)m<-yhat
#list(Strength.Assoc=sqrt(E.pow),Explanatory.Power=E.pow,yhat=m)
m
}