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