[R] Logistic regression/Cut point? predict ??

[Simon Knapp]

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.