R/lplotv2.R

lplotv2 <-
function(x,y,span=.75,pyhat=FALSE,eout=FALSE,xout=FALSE,outfun=out,plotit=TRUE,
expand=.5,low.span=2/3,varfun=pbvar,cor.op=FALSE,cor.fun=pbcor,ADJ=FALSE,nboot=20,
scale=FALSE,xlab="X",ylab="Y",zlab="",theta=50,phi=25,family="gaussian",
duplicate="error",pr=TRUE,SEED=TRUE,ticktype="simple"){
#
# Plot regression surface using LOESS
#
# low.span is the span when lowess is used and there is one predictor
# span is the span when loess is used with two or more predictors
# pyhat=T will return Y hat values
# eout=T will eliminate outliers
# xout=T  will eliminate points where X is an outliers
# family="gaussian"; see the description of the built-in function loess
#
# duplicate="error"
# In some situations where duplicate values occur, when plotting with
# two predictors, it is necessary to set duplicate="strip"
#
st.adj=NULL
e.adj=NULL
if(ADJ){
if(SEED)set.seed(2)
}
si=1
library(stats)
x<-as.matrix(x)
if(!is.matrix(x))stop("x is not a matrix")
d<-ncol(x)
if(d>=2){
library(akima)
if(ncol(x)==2 && !scale){
if(pr){
print("scale=F is specified.")
print("If there is dependence, might use scale=T")
}}
m<-elimna(cbind(x,y))
x<-m[,1:d]
y<-m[,d+1]
if(eout && xout)stop("Can't have both eout and xout = F")
if(eout){
flag<-outfun(m,plotit=FALSE)$keep
m<-m[flag,]
}
if(xout){
flag<-outfun(x,plotit=FALSE)$keep
m<-m[flag,]
}
x<-m[,1:d]
y<-m[,d+1]
if(d==2)fitr<-fitted(loess(y~x[,1]*x[,2],span=span,family=family))
if(d==3)fitr<-fitted(loess(y~x[,1]*x[,2]*x[,3],span=span,family=family))
if(d==4)fitr<-fitted(loess(y~x[,1]*x[,2]*x[,3]*x[,4],span=span,family=family))
if(d>4)stop("Can have at most four predictors")
last<-fitr
if(d==2 && plotit){
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 R function interp
mkeep<-x[iout>=1,]
fitr<-interp(mkeep[,1],mkeep[,2],fitr,duplicate=duplicate)
persp(fitr,theta=theta,phi=phi,xlab=xlab,ylab=ylab,zlab=zlab,expand=expand,
scale=scale,ticktype=ticktype)
}}
if(d==1){
m<-elimna(cbind(x,y))
x<-m[,1:d]
y<-m[,d+1]
if(eout && xout)stop("Can't have both eout and xout = F")
if(eout){
flag<-outfun(m)$keep
m<-m[flag,]
}
if(xout){
flag<-outfun(x)$keep
m<-m[flag,]
}
x<-m[,1:d]
y<-m[,d+1]
if(plotit){
plot(x,y,xlab=xlab,ylab=ylab)
lines(lowess(x,y,f=low.span))
}
yyy<-lowess(x,y)$y
xxx<-lowess(x,y)$x
if(d==1){
ordx=order(xxx)
yord=yyy[ordx]
flag=NA
for (i in 2:length(yyy))flag[i-1]=sign(yord[i]-yord[i-1])
if(sum(flag)<0)si=-1
}
last<-yyy
chkit<-sum(duplicated(x))
if(chkit>0){
last<-rep(1,length(y))
for(j in 1:length(yyy)){
for(i in 1:length(y)){
if(x[i]==xxx[j])last[i]<-yyy[j]
}}
}
}
E.power<-1
if(!cor.op)E.power<-varfun(last[!is.na(last)])/varfun(y)
if(cor.op || E.power>=1){
if(d==1){
xord<-order(x)
E.power<-cor.fun(last,y[xord])$cor^2
}
if(d>1)E.power<-cor.fun(last,y)$cor^2
}
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=lplot.sub(x[data1[i,],],y[data2[i,]],plotit=FALSE,pr=FALSE)
val[i]=temp$Explanatory.power
}
vindt=median(val)
v2indt=median(sqrt(val))
st.adj=(sqrt(E.power)-max(c(0,v2indt)))/(1-max(c(0,v2indt)))
e.adj=(E.power-max(c(0,vindt)))/(1-max(c(0,vindt)))
st.adj=max(c(0,st.adj))
e.adj=max(c(0,e.adj))
}
if(!pyhat)last <- NULL
list(Strength.Assoc=si*sqrt(E.power),Explanatory.power=E.power,
Strength.Adj=st.adj,Explanatory.Adj=e.adj,yhat.values=last)
}
musto101/wilcox_R documentation built on May 23, 2019, 10:52 a.m.