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(predcutmid,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.
