wilcox_R: What the package does (short line)

gamplotv2 <-
function(x,y,sop=FALSE,pyhat=FALSE,eout=FALSE,xout=FALSE,outfun=out,plotit=TRUE,
varfun=pbvar,xlab="X",ylab="",zlab="",theta=50,phi=25,expand=.5,SCALE=FALSE,
cor.fun=pbcor,ADJ=FALSE,nboot=20,pr=TRUE,SEED=TRUE,ticktype="simple"){
#
# Plot regression surface using generalized additive model
#
# sop=F, use lowess
# sop=T, use splines
#
if(ADJ){
if(SEED)set.seed(2)
}
if(pr){
if(!ADJ){
print("To get adjusted estimates of strength of association, use ADJ=T")
print("The strength of association is estimated under independence")
print(" and then rescaled")
}}
library(akima)
library(mgcv)
x<-as.matrix(x)
np<-ncol(x)
np1<-np+1
if(ncol(x)>4)stop("x should have at most four columns of data")
m<-elimna(cbind(x,y))
if(xout && eout)stop("Can't have xout=eout=T")
if(eout){
flag<-outfun(m)$keep
m<-m[flag,]
}
if(xout){
flag<-outfun(x,plotit=FALSE)$keep
m<-m[flag,]
}
x<-m[,1:np]
x=as.matrix(x)
y<-m[,np1]
if(!sop){
if(ncol(x)==1)fitr<-fitted(gam(y~x[,1]))
if(ncol(x)==2)fitr<-fitted(gam(y~x[,1]+x[,2]))
if(ncol(x)==3)fitr<-fitted(gam(y~x[,1]+x[,2]+x[,3]))
if(ncol(x)==4)fitr<-fitted(gam(y~x[,1]+x[,2]+x[,3]+x[,4]))
}
if(sop){
if(ncol(x)==1)fitr<-fitted(gam(y~s(x[,1])))
if(ncol(x)==2)fitr<-fitted(gam(y~s(x[,1])+s(x[,2])))
if(ncol(x)==3)fitr<-fitted(gam(y~s(x[,1])+s(x[,2])+s(x[,3])))
if(ncol(x)==4)fitr<-fitted(gam(y~s(x[,1])+s(x[,2])+s(x[,3])+s(x[,4])))
}
last<-fitr
if(plotit){
if(ncol(x)==1){
plot(x,fitr,xlab=xlab,ylab=ylab)
}
if(ncol(x)==2){
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(x[i,]==x[k,])==2)iout[k]<-0
}
fitr<-fitr[iout>=1] # Eliminate duplicate points in the x-y plane
#                 This is necessary when doing three dimensional plots
#                 with the S-PLUS function interp
mkeep<-x[iout>=1,]
fitr<-interp(mkeep[,1],mkeep[,2],fitr)
persp(fitr,theta=theta,phi=phi,expand=expand,xlab="x1",ylab="x2",zlab="",
scale=scale,ticktype=ticktype)
}
}
top=varfun(last)
ep=top/varfun(y)
if(ep>=1)ep=cor.fun(last,y)$cor^2
eta=sqrt(ep)
st.adj=NULL
e.adj=NULL
if(ADJ){
x=as.matrix(x)
val=NA
n=length(y)
data1<-matrix(sample(n,size=n*nboot,replace=TRUE),nrow=nboot)
data2<-matrix(sample(n,size=n*nboot,replace=TRUE),nrow=nboot)
for(i in 1:nboot){
temp=gamplotv2.sub(x[data1[i,],],y[data2[i,]],plotit=FALSE)
val[i]=temp$Explanatory.power
}
vindt=median(val)
v2indt=median(sqrt(val))
st.adj=(sqrt(ep)-max(c(0,v2indt)))/(1-max(c(0,v2indt)))
e.adj=(ep-max(c(0,vindt)))/(1-max(c(0,vindt)))
st.adj=max(c(0,st.adj))
e.adj=max(c(0,e.adj))
}
eta=as.matrix(eta)
ep=as.matrix(ep)
dimnames(eta)=NULL
dimnames(ep)=NULL
eta=eta[1]
ep=ep[1]
list(Strength.Assoc=eta,Explanatory.power=ep,
Strength.Adj=st.adj,Explanatory.Adj=e.adj)
}