Problems getting slope and intercept to change when do multiple reps.

classic Classic list List threaded Threaded
1 message Options
Reply | Threaded
Open this post in threaded view
|

Problems getting slope and intercept to change when do multiple reps.

Adel Powell
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]]

______________________________________________
[hidden email] 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.