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]]
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
