R version 3.2.3 (2015-12-10) -- "Wooden Christmas-Tree"
Copyright (C) 2015 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> x <- array(list(1,10403816,5.709088235,1,10407940,5.768005384,1,10403821,5.74483159,1,10403825,5.948173502,1,10548415,6.019050695,1,10414953,5.844783999,1,10485622,6.161365478,1,10548497,6.198138128,1,10549647,6.233888197,1,10548422,6.233675026,0,10412211,6.642731362,0,10420308,6.376033012,0,10538979,6.128666093,0,10428576,6.117440006,0,10538993,5.931002217,0,10414885,6.14089668,0,10607738,6.231257962,0,10408935,6.38834449,0,10415015,5.998148415,0,10403978,6.078712279),dim=c(3,20),dimnames=list(c('Y','X1','X2'),1:20))
>  y <- array(NA,dim=c(3,20),dimnames=list(c('Y','X1','X2'),1:20))
>  for (i in 1:dim(x)[1])
+  {
+  	for (j in 1:dim(x)[2])
+  	{
+  		y[i,j] <- as.numeric(x[i,j])
+  	}
+  }
> #'GNU S' R Code compiled by R2WASP v. 1.2.327 (Mon, 27 Jul 2015 15:23:40 +0100)
> #Author: root
> #To cite this work: Wessa P., (2015), Bias Reduced Logistic Regression (v1.0.5) in Free Statistics Software (v$_version), Office for Research Development and Education, URL http://www.wessa.net/rwasp_logisticregression.wasp/
> #Source of accompanying publication: Office for Research, Development, and Education
> #
> library(brglm)
Loading required package: profileModel
> roc.plot <- function (sd, sdc, newplot = TRUE, ...)
+ {
+ sall <- sort(c(sd, sdc))
+ sens <- 0
+ specc <- 0
+ for (i in length(sall):1) {
+ sens <- c(sens, mean(sd >= sall[i], na.rm = T))
+ specc <- c(specc, mean(sdc >= sall[i], na.rm = T))
+ }
+ if (newplot) {
+ plot(specc, sens, xlim = c(0, 1), ylim = c(0, 1), type = 'l',
+ xlab = '1-specificity', ylab = 'sensitivity', main = 'ROC plot', ...)
+ abline(0, 1)
+ }
+ else lines(specc, sens, ...)
+ npoints <- length(sens)
+ area <- sum(0.5 * (sens[-1] + sens[-npoints]) * (specc[-1] -
+ specc[-npoints]))
+ lift <- (sens - specc)[-1]
+ cutoff <- sall[lift == max(lift)][1]
+ sensopt <- sens[-1][lift == max(lift)][1]
+ specopt <- 1 - specc[-1][lift == max(lift)][1]
+ list(area = area, cutoff = cutoff, sensopt = sensopt, specopt = specopt)
+ }
> roc.analysis <- function (object, newdata = NULL, newplot = TRUE, ...)
+ {
+ if (is.null(newdata)) {
+ sd <- object$fitted[object$y == 1]
+ sdc <- object$fitted[object$y == 0]
+ }
+ else {
+ sd <- predict(object, newdata, type = 'response')[newdata$y ==
+ 1]
+ sdc <- predict(object, newdata, type = 'response')[newdata$y ==
+ 0]
+ }
+ roc.plot(sd, sdc, newplot, ...)
+ }
> hosmerlem <- function (y, yhat, g = 10)
+ {
+ cutyhat <- cut(yhat, breaks = quantile(yhat, probs = seq(0,
+ 1, 1/g)), include.lowest = T)
+ obs <- xtabs(cbind(1 - y, y) ~ cutyhat)
+ expect <- xtabs(cbind(1 - yhat, yhat) ~ cutyhat)
+ chisq <- sum((obs - expect)^2/expect)
+ P <- 1 - pchisq(chisq, g - 2)
+ c('X^2' = chisq, Df = g - 2, 'P(>Chi)' = P)
+ }
> x <- as.data.frame(t(y))
> r <- brglm(x)
> summary(r)

Call:
brglm(formula = x)


Coefficients:
              Estimate Std. Error z value Pr(>|z|)  
(Intercept) -4.290e+01  7.647e+01  -0.561   0.5748  
X1           7.174e-06  7.924e-06   0.905   0.3652  
X2          -5.295e+00  3.082e+00  -1.718   0.0858 .
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 20.884  on 19  degrees of freedom
Residual deviance: 21.462  on 17  degrees of freedom
Penalized deviance: -1.38545 
AIC:  27.462 

> rc <- summary(r)$coeff
> try(hm <- hosmerlem(y[1,],r$fitted.values),silent=T)
> try(hm,silent=T)
       X^2         Df    P(>Chi) 
11.4494380  8.0000000  0.1775101 
> postscript(file="/var/wessaorg/rcomp/tmp/190cq1457420710.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) 
> ra <- roc.analysis(r)
> dev.off()
null device 
          1 
> te <- array(0,dim=c(2,99))
> for (i in 1:99) {
+ threshold <- i / 100
+ numcorr1 <- 0
+ numfaul1 <- 0
+ numcorr0 <- 0
+ numfaul0 <- 0
+ for (j in 1:length(r$fitted.values)) {
+ if (y[1,j] > 0.99) {
+ if (r$fitted.values[j] >= threshold) numcorr1 = numcorr1 + 1 else numfaul1 = numfaul1 + 1
+ } else {
+ if (r$fitted.values[j] < threshold) numcorr0 = numcorr0 + 1 else numfaul0 = numfaul0 + 1
+ }
+ }
+ te[1,i] <- numfaul1 / (numfaul1 + numcorr1)
+ te[2,i] <- numfaul0 / (numfaul0 + numcorr0)
+ }
> postscript(file="/var/wessaorg/rcomp/tmp/2plom1457420710.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) 
> op <- par(mfrow=c(2,2))
> plot((1:99)/100,te[1,],xlab='Threshold',ylab='Type I error', main='1 - Specificity')
> plot((1:99)/100,te[2,],xlab='Threshold',ylab='Type II error', main='1 - Sensitivity')
> plot(te[1,],te[2,],xlab='Type I error',ylab='Type II error', main='(1-Sens.) vs (1-Spec.)')
> plot((1:99)/100,te[1,]+te[2,],xlab='Threshold',ylab='Sum of Type I & II error', main='(1-Sens.) + (1-Spec.)')
> par(op)
> dev.off()
null device 
          1 
> 
> #Note: the /var/wessaorg/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/wessaorg/rcomp/createtable")
> 
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a,'Coefficients of Bias-Reduced Logistic Regression',5,TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a,'Variable',header=TRUE)
> a<-table.element(a,'Parameter',header=TRUE)
> a<-table.element(a,'S.E.',header=TRUE)
> a<-table.element(a,'t-stat',header=TRUE)
> a<-table.element(a,'2-sided p-value',header=TRUE)
> a<-table.row.end(a)
> for (i in 1:length(rc[,1])) {
+ a<-table.row.start(a)
+ a<-table.element(a,labels(rc)[[1]][i],header=TRUE)
+ a<-table.element(a,rc[i,1])
+ a<-table.element(a,rc[i,2])
+ a<-table.element(a,rc[i,3])
+ a<-table.element(a,2*(1-pt(abs(rc[i,3]),r$df.residual)))
+ a<-table.row.end(a)
+ }
> a<-table.end(a)
> table.save(a,file="/var/wessaorg/rcomp/tmp/3h85k1457420710.tab") 
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a,'Summary of Bias-Reduced Logistic Regression',2,TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a,'Deviance',1,TRUE)
> a<-table.element(a,r$deviance)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a,'Penalized deviance',1,TRUE)
> a<-table.element(a,r$penalized.deviance)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a,'Residual Degrees of Freedom',1,TRUE)
> a<-table.element(a,r$df.residual)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a,'ROC Area',1,TRUE)
> a<-table.element(a,ra$area)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a,'Hosmer–Lemeshow test',2,TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a,'Chi-square',1,TRUE)
> phm <- array('NA',dim=3)
> for (i in 1:3) { try(phm[i] <- hm[i],silent=T) }
> a<-table.element(a,phm[1])
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a,'Degrees of Freedom',1,TRUE)
> a<-table.element(a,phm[2])
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a,'P(>Chi)',1,TRUE)
> a<-table.element(a,phm[3])
> a<-table.row.end(a)
> a<-table.end(a)
> table.save(a,file="/var/wessaorg/rcomp/tmp/41sf51457420710.tab") 
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a,'Fit of Logistic Regression',4,TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a,'Index',1,TRUE)
> a<-table.element(a,'Actual',1,TRUE)
> a<-table.element(a,'Fitted',1,TRUE)
> a<-table.element(a,'Error',1,TRUE)
> a<-table.row.end(a)
> for (i in 1:length(r$fitted.values)) {
+ a<-table.row.start(a)
+ a<-table.element(a,i,1,TRUE)
+ a<-table.element(a,y[1,i])
+ a<-table.element(a,r$fitted.values[i])
+ a<-table.element(a,y[1,i]-r$fitted.values[i])
+ a<-table.row.end(a)
+ }
> a<-table.end(a)
> table.save(a,file="/var/wessaorg/rcomp/tmp/566re1457420710.tab") 
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a,'Type I & II errors for various threshold values',3,TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a,'Threshold',1,TRUE)
> a<-table.element(a,'Type I',1,TRUE)
> a<-table.element(a,'Type II',1,TRUE)
> a<-table.row.end(a)
> for (i in 1:99) {
+ a<-table.row.start(a)
+ a<-table.element(a,i/100,1,TRUE)
+ a<-table.element(a,te[1,i])
+ a<-table.element(a,te[2,i])
+ a<-table.row.end(a)
+ }
> a<-table.end(a)
> table.save(a,file="/var/wessaorg/rcomp/tmp/6wzpd1457420710.tab") 
> 
> try(system("convert tmp/190cq1457420710.ps tmp/190cq1457420710.png",intern=TRUE))
character(0)
> try(system("convert tmp/2plom1457420710.ps tmp/2plom1457420710.png",intern=TRUE))
character(0)
> 
> 
> proc.time()
   user  system elapsed 
  0.962   0.145   1.116