R/runm3d.R

runm3d <-
function(x,y,theta=50,phi=25,fr=.8,tr=.2,plotit=TRUE,pyhat=FALSE,nmin=0,
expand=.5,scale=FALSE,zscale=FALSE,xout=FALSE,outfun=out,eout=FALSE,xlab="X",ylab="Y",zlab="",
pr=TRUE,SEED=TRUE,ticktype="simple"){
#
# running mean using interval method
#
# fr controls amount of smoothing
# tr is the amount of trimming
# x is an n by p matrix of predictors.
#
#  Rows of data with missing values are automatically removed.
#
# When plotting, theta and phi can be used to change
# the angle at which the plot is viewed.
#
#  theta is the azimuthal direction and phi the colatitude
#   expand controls relative length of z-axis
#
library(MASS)
library(akima)
if(plotit){
if(pr){
print("Note: when there is independence, scale=F is probably best")
print("When there is dependence, scale=T is probably best")
}}
if(!is.matrix(x))stop("x should be a matrix")
if(nrow(x) != length(y))stop("number of rows of x should equal length of y")
temp<-cbind(x,y)
p<-ncol(x)
p1<-p+1
temp<-elimna(temp) # Eliminate any rows with missing values.
if(xout){
keepit<-rep(T,nrow(x))
flag<-outfun(x,plotit=FALSE)$out.id
keepit[flag]<-F
x<-x[keepit,]
y<-y[keepit]
}
if(zscale){
for(j in 1:p1){
temp[,j]<-(temp[,j]-median(temp[,j]))/mad(temp[,j])
}}
x<-temp[,1:p]
y<-temp[,p1]
pyhat<-as.logical(pyhat)
plotit<-as.logical(plotit)
if(SEED)set.seed(12)
m<-cov.mve(x)
iout<-c(1:nrow(x))
rmd<-1 # Initialize rmd
nval<-1
for(i in 1:nrow(x))rmd[i]<-mean(y[near3d(x,x[i,],fr,m)],tr)
for(i in 1:nrow(x))nval[i]<-length(y[near3d(x,x[i,],fr,m)])
if(plotit){
if(ncol(x)!=2)stop("When plotting, x must be an n by 2 matrix")
fitr<-rmd[nval>nmin]
y<-y[nval>nmin]
x<-x[nval>nmin,]
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,]
fit<-interp(mkeep[,1],mkeep[,2],fitr)
persp(fit,theta=theta,phi=phi,xlab=xlab,ylab=ylab,zlab=zlab,expand=expand,
scale=scale,ticktype=ticktype)
}
last<-"Done"
if(pyhat)last<-rmd
last
}
musto101/wilcox_R documentation built on May 23, 2019, 10:52 a.m.