library(ROCR)
n <- 1000
fitglm <- function(iteration,intercept,sigma,tau,beta){
x <- rnorm(n,0,sigma)
ystar <- intercept+beta*x
z <- rbinom(n,1,plogis(ystar))
xerr <- x + rnorm(n,0,tau)
model<-glm(z ~ xerr, family=binomial(logit))
*int*<-coef(model)[1]
*slope*<-coef(model)[2] # when add error you are suppose to get slightly
bias slope. However when I change the beta in the original X, I am not
getting the save average slope as output? strange?
pred<-predict(model,type="response")
cutp<-.5
result<-ifelse(pred>cutp,1,0)
rocpreds<-prediction(result,z)
auc<-performance(rocpreds,"auc")@y.values
accuracy<-length(which(result==z))/length(z)
tn<- sum(z==0 & result==0) # True Negative
fp<- sum(z==0 & result==1) # False Positive
tp<- sum(z==1 & result==1) # True Positive
fn<- sum(z==1 & result==0) # False Negative
sensitivity<- tp/(tp+fn)
specificity<-tn/(tn+fp)
output<-c(*int,slope*,cutp,accuracy,auc,sensitivity,specificity,iteration)
names(output)<-c("Intercept","Slope","CutPoint","Accuracy","AUC","Sensitivity","Specificity","iteration")
return(output)
}
y<-fitglm(1,2,1,2,4)
y
Output<-t(sapply(1:10, function(x) fitglm(x,intercept=2, sigma=1,tau=2,beta=
4)))
apply(output,2, function(x) mean(unlist(x)))
apply(output,2, function(x) var(unlist(x)))
Output<-t(sapply(1:500, function(x) fitglm(x,intercept=0, sigma=1,tau=.25,
beta=1)))
apply(output,2, function(x) mean(unlist(x)))
apply(output,2, function(x) var(unlist(x)))
Output<-t(sapply(1:500, function(x) fitglm(x,intercept=2, sigma=1,tau=.25,
beta=6)))
apply(output,2, function(x) mean(unlist(x)))
apply(output,2, function(x) var(unlist(x)))
[[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.
