n <- 10; # Sample size
fitglm <- function(sigma,tau){
x <- rnorm(n,0,sigma)
intercept <- 0
beta <- 0
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]
pred<-predict(model)
result<-ifelse(pred>.5,1,0)
accuracy<-length(which(result==z))/length(z)
accuracy
rocpreds<-prediction(result,z)
auc<-performance(rocpreds,"auc")@y.values
fp<-performance(rocpreds,"sens")
sentiv<-slot(fp,"y.values")[[1]]
sentiv<-sentiv[2]
sentiv
fp2<-performance(rocpreds,"spec")
specs<-slot(fp2,"y.values")[[1]]
specs<-specs[2]
specs
output<-c(int,slope,.5,accuracy,auc,sentiv,specs)
names(output)<-c("Intercept","Slope","CutPoint","Accuracy","AUC","Sentivity","Specificity")
return(output)
}
y<-fitglm(2,1)
y
I corrected the code. I am put everything in a loop because I am running
monte carlo reps outside of this later.
It works but the value returned is wrong.
Here is the manual manipulation:
*> x <- rnorm(10,0,2)*
* *
*> intercept <- 0*
* *
*> beta <- 5*
* *
*> ystar <- intercept+beta*x*
* *
*> ystar*
[1] 16.5436337 7.7740329 -10.1653928 -2.8338118 -21.5410780 2.6902171
5.1156558 5.0729933 -10.8556430 0.9208434
*> test <- plogis(ystar)*
* *
*> test*
[1] 9.999999e-01 9.995797e-01 3.847772e-05 5.552417e-02 4.413963e-10
9.364469e-01 9.940338e-01 9.937753e-01 1.929504e-05 7.152139e-01
*> z <- rbinom(10,1,plogis(ystar))*
* *
*> z*
[1] 1 1 0 0 0 1 1 1 0 1
*> xerr <- x + rnorm(10,0,1) *
* *
*> xerr*
[1] 0.5610573 3.1741687 -2.3915066 -0.2546224 -4.1790037 -1.4387786
1.4211227
-1.1141176 -1.6230087 0.7595021
*> model<-glm(z ~ xerr, family=binomial(logit))*
* *
*> model*
Call: glm(formula = z ~ xerr, family = binomial(logit))
Coefficients:
(Intercept) xerr
1.500 1.309
*> int<-coef(model)[1]*
* *
*> slope<-coef(model)[2]*
* *
*> pred1<-predict(model)*
* *
*> pred2<-predict(model,type="response")*
* *
*> pred1*
1 2 3 4 5
6
7 8 9 10
2.23478077 5.65499178 -1.62972799 1.16716581 -3.96932102 -0.38273530
3.36049056 0.04220225 -0.62386760 2.49451832
*> pred2*
1 2 3 4 5 6
7
8 9 10
0.90332965 0.99651221 0.16386763 0.76263234 0.01853617 0.40546735
0.96644669 0.51054900 0.34890234 0.92375664
*> result<-ifelse(pred2>.5,1,0) *
* *
*> result*
1 2 3 4 5 6 7 8 9 10
1 1 0 1 0 0 1 1 0 1
*> accuracy<-length(which(result==z))/length(z)*
* *
*> accuracy*
[1] 0.8
*> rocpreds<-prediction(result,z)*
* *
*> rocpreds*
*> auc<-performance(rocpreds,"auc")@y.values*
* *
*> auc*
[[1]]
[1] 0.7916667
* > fp<-performance(rocpreds,"sens")*
* *
*> sentiv<-slot(fp,"y.values")[[1]]*
* *
*> sentiv<-sentiv[2]*
* *
*> sentiv*
[1] 0.8333333
* *
*> fp2<-performance(rocpreds,"spec")*
* *
*> specs<-slot(fp2,"y.values")[[1]]*
* *
*> specs*
[1] 1.00 0.75 0.00
*> specs<-specs[2]*
* *
*> specs*
[1] 0.75
On Thu, Oct 25, 2012 at 3:43 PM, Berend Hasselman <b...@xs4all.nl> wrote:
>
> On 25-10-2012, at 21:28, Adel Powell wrote:
>
> > I am running my code in a loop and it does not work but when I run it
> > outside the loop I get the values I want.
> >
> > n <- 1000; # Sample size
> >
> > fitglm <- function(sigma,tau){
> > x <- rnorm(n,0,sigma)
> > intercept <- 0
> > beta <- 0
> > 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]
> > pred<-predict(model)
> >
> > result<-ifelse(pred>.5,1,0)
> >
> > accuracy<-length(which(result==z))/length(z)
> > accuracy
> >
> > rocpreds<-prediction(result,z)
> > auc<-performance(rocpreds,"auc")@y.values
> > sentiv<-performance(rocpreds,"sens")@y.values
> > sentiv<-slot(fp,"y.values")[[1]]
> > sentiv<-sentiv[2]
> > sentiv
> > specs<-performance(rocpreds,"spec")@y.values
> > specs<-slot(fp2,"y.values")[[1]]
> > specs<-specs[2]
> > specs
> > output<-c(int,slope,.5,accuracy,auc,sentiv,specs)
> >
> >
> names(output)<-c("Intercept","Slope","CutPoint","Accuracy","AUC","Sentivity",Specificity")
>
> A missing " before Specificity?
>
> > return(output)
> >
> > }
> >
> > y<-fitglm(.05,1)
> > y
> >
>
> Running this after correction of the missing " one gets en error
>
> Error in fitglm(0.05, 1) : could not find function "prediction"
>
> How are you using a loop?
> Your example is not reproducible.
>
> Berend
>
>
> >
> > The code runs without the sentiv and specs but when I remove the loop i
> can
> > get the sensitivity and spec. values ???
> >
> > [[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.
>
>
[[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.
