What do you mean by "at x equal zero"? On Sun, Oct 21, 2012 at 8:37 AM, Adel Powell <powella...@gmail.com> wrote: > 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.
______________________________________________ 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.