R/disker.R

disker <-
function(x,y,z=NA,plotit=FALSE,op=1){
#
#  Estimate apparent effect size
#  using probability of correct classification based on values in
#  first group.
#
#  A "CORRECT" classification is the event of deciding
#  that an observation
#  from the first group did indeed come from the first group based
#  on a kernel density estimate of the distributions.
#  The function returns the
#  proportion of correctly classified observations (phat).
#
#  The function also returns a vector of 0s and 1s (in zhat)
#  indicating  whether values in z would be
#  classified as coming from the first group.
#
# op=1, use Rosenblatt's shifted histogram version of kernel estimate
# op=2, use adaptive kernel estimate with initial estimate based
#       on expected frequency curve.
#
xsort<-sort(x)
ysort<-sort(y)
xhat<-0
yhat<-0
yyhat<-0
if(op==1){
for(i in 1:length(xsort))xhat[i]<-kerden(x,0,xsort[i])
for(i in 1:length(xsort))yhat[i]<-kerden(y,0,xsort[i])
}
if(op==2){
xhat<-akerd(x,pts=xsort,pyhat=TRUE,plotit=FALSE)
yhat<-akerd(y,pts=xsort,pyhat=TRUE,plotit=FALSE)
}
yhat[is.na(yhat)]<-0
if(plotit){
if(op==1){
for(i in 1:length(ysort))yyhat[i]<-kerden(y,0,ysort[i])
}
if(op==2)yyhat<-akerd(y,pts=ysort,plotit=FALSE,pyhat=T)
plot(c(xsort,ysort),c(xhat,yyhat),type="n",xlab="",ylab="")
lines(xsort,xhat)
lines(ysort,yyhat)
}
#
# Compute apparent error
#
phat<-sum(xhat>yhat)/length(x)
zhat<-NA
if(!is.na(z[1])){
#
#  Make decisions for the data in z,
#  set zhat=1 if decide it came from
#  group 1.
#
zxhat<-0
zyhat<-0
zhat<-0
if(op==2){
zxhat<-akerd(x,pts=z,pyhat=TRUE,plotit=FALSE)
zyhat<-akerd(y,pts=z,pyhat=TRUE,plotit=FALSE)
}
for(i in 1:length(z)){
if(op==1){
zxhat[i]<-kerden(x,0,z[i])
zyhat[i]<-kerden(y,0,z[i])
}
zhat[i]<-1
if(is.na(zxhat[i]) || is.na(zyhat[i])){
# Missing values,
# data can't be used to make a decision,
# so make a random decision about whether a value
# came from first group.
arb<-runif(1)
zhat[i]<-1
if(arb < .5)zhat[i]<-0
}
else
if(zxhat[i]<zyhat[i])zhat[i]<-0
}
}
list(phat=phat,zhat=zhat)
#phat is the apparent probability  of a correct classification
}
musto101/wilcox_R documentation built on May 23, 2019, 10:52 a.m.