I am new to R and I am trying to do a monte carlo simulation where I
generate data and interject error then test various cut points; however, my
output was garbage (at x equal zero, I did not get .50)
I am basically testing the performance of classifiers.

Here is the code:
n <- 1000; # Sample size

fitglm <- function(sigma,tau){
    x <- rnorm(n,0,sigma)
    intercept <- 0
    beta <- 5
   * ystar <- intercept+beta*x*
   * z <- rbinom(n,1,plogis(ystar))*    *# I believe plogis accepts the a
+bx augments and return the  e^x/(1+e^x)  which is then used to generate 0
and 1 data*
    xerr <- x + rnorm(n,0,tau)    # error is added here
    model<-glm(z ~ xerr, family=binomial(logit))
    int<-coef(model)[1]
    slope<-coef(model)[2]
    pred<-predict(model)  #this gives me the a+bx data for new error?  I
know I can add type= response to get the probab. but only e^x not *e^x/(1+e^x)
*

    pi1hat<-length(z[which(z==1)]/length(z)) My cut point is calculated  is
the proportion of 0s to 1.
    pi0hat<-length(z[which(z==0)]/length(z))

    cutmid <- log(pi0hat/pi1hat)
    pred<-ifelse(pred>cutmid,1,0) * I am not sure if I need to compare
these two. I think this is an error.
*
    accuracy<-length(which(pred==z))/length(z)
    accuracy

    rocpreds<-prediction(pred,z)
    auc<-performance(rocpreds,"auc")@y.values

    output<-c(int,slope,cutmid,accuracy,auc)
    names(output)<-c("Intercept","Slope","CutPoint","Accuracy","AUC")
    return(output)

}

y<-fitglm(.05,1)
y

nreps <- 500;
output<-data.frame(matrix(rep(NA),nreps,6,ncol=6))

mysigma<-.5
mytau<-.1

i<-1

for(j in 1:nreps) {
    output[j,1:5]<-fitglm(mysigma,mytau)
    output[j,6]<-j
}

names(output)<-c("Intercept","Slope","CutPoint","Accuracy","AUC","Iteration")

apply(output,2, mean)
apply(output,2, var)

	[[alternative HTML version deleted]]

