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 <bhh@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]]


