Repository of Reproducible Computations

Free Statistics

of Irreproducible Research!

Author's title

Author

*The author of this computation has been verified*

R Software Module

rwasp_logisticregression.wasp

Title produced by software

Bias-Reduced Logistic Regression

Date of computation

Tue, 21 Dec 2010 18:56:29 +0000

Cite this page as follows

Statistical Computations at FreeStatistics.org, Office for Research Development and Education, URL https://freestatistics.org/blog/index.php?v=date/2010/Dec/21/t1292957645d1tyecinz69m4ul.htm/, Retrieved Sat, 18 May 2024 14:42:52 +0000

Statistical Computations at FreeStatistics.org, Office for Research Development and Education, URL https://freestatistics.org/blog/index.php?pk=113827, Retrieved Sat, 18 May 2024 14:42:52 +0000

QR Codes:

Paste this QR Code to cite your computation.

Original text written by user:

IsPrivate?

No (this computation is public)

User-defined keywords

Estimated Impact

122

Family? (F = Feedback message, R = changed R code, M = changed R Module, P = changed Parameters, D = changed Data)

-     [Multiple Regression] [ws sleep] [2010-12-13 13:22:04] [df61ce38492c371f14c407a12b3bb2eb]
- RMP   [Bias-Reduced Logistic Regression] [ws sleep] [2010-12-19 10:12:57] [df61ce38492c371f14c407a12b3bb2eb]
-    D      [Bias-Reduced Logistic Regression] [] [2010-12-21 18:56:29] [76f6fcd790878de142f355e7238b5c71] [Current]

Feedback Forum

Post a new message

Dataseries X:

Download CSV

Histogram

Boxplots

-999,00	6654,00	3,00
6,30	1,00	3,00
-999,00	3,39	1,00
-999,00	0,92	3,00
2,10	2547,00	4,00
9,10	10,55	4,00
15,80	0,02	1,00
5,20	160,00	4,00
10,90	3,30	1,00
8,30	52,16	1,00
11,00	0,43	4,00
3,20	465,00	5,00
7,60	0,55	2,00
-999,00	187,10	5,00
6,30	0,08	1,00
8,60	3,00	2,00
6,60	0,79	2,00
9,50	0,20	2,00
4,80	1,41	1,00
12,00	60,00	1,00
-999,00	529,00	5,00
3,30	27,66	5,00
11,00	0,12	2,00
-999,00	207,00	1,00
4,70	85,00	1,00
-999,00	36,33	1,00
10,40	0,10	3,00
7,40	1,04	4,00
2,10	521,00	5,00
-999,00	100,00	1,00
-999,00	35,00	4,00
7,70	0,01	4,00
17,90	0,01	1,00
6,10	62,00	1,00
8,20	0,12	1,00
8,40	1,35	3,00
11,90	0,02	3,00
10,80	0,05	3,00
13,80	1,70	1,00
14,30	3,50	1,00
-999,00	250,00	5,00
15,20	0,48	2,00
10,00	10,00	4,00
11,90	1,62	2,00
6,50	192,00	4,00
7,50	2,50	5,00
-999,00	4,29	2,00
10,60	0,28	3,00
7,40	4,24	1,00
8,40	6,80	2,00
5,70	0,75	2,00
4,90	3,60	3,00
-999,00	14,83	5,00
3,20	55,50	5,00
-999,00	1,40	2,00
8,10	0,06	2,00
11,00	0,90	2,00
4,90	2,00	3,00
13,20	0,10	2,00
9,70	4,19	4,00
12,80	3,50	1,00
-999,00	4,05	1,00

Summary of computational transaction
Raw Input	view raw input (R code)
Raw Output	view raw output of R engine
Computing time	1 seconds
R Server	'Gwilym Jenkins' @ 72.249.127.135
R Engine error message	Error in eval(expr, envir, enclos) : y values must be 0 <= y <= 1 Calls: brglm -> fit.proc -> eval -> eval Execution halted

\begin{tabular}{lllllllll}
\hline
Summary of computational transaction \tabularnewline
Raw Input & view raw input (R code)  \tabularnewline
Raw Output & view raw output of R engine  \tabularnewline
Computing time & 1 seconds \tabularnewline
R Server & 'Gwilym Jenkins' @ 72.249.127.135 \tabularnewline
R Engine error message & Error in eval(expr, envir, enclos) : y values must be 0 <= y <= 1
Calls: brglm -> fit.proc -> eval -> eval
Execution halted
 \tabularnewline
\hline
\end{tabular}
%Source: https://freestatistics.org/blog/index.php?pk=113827&T=0

[TABLE]
[ROW][C]Summary of computational transaction[/C][/ROW]
[ROW][C]Raw Input[/C][C]view raw input (R code) [/C][/ROW]
[ROW][C]Raw Output[/C][C]view raw output of R engine [/C][/ROW]
[ROW][C]Computing time[/C][C]1 seconds[/C][/ROW]
[ROW][C]R Server[/C][C]'Gwilym Jenkins' @ 72.249.127.135[/C][/ROW]
[ROW][C]R Engine error message[/C][C]Error in eval(expr, envir, enclos) : y values must be 0 <= y <= 1
Calls: brglm -> fit.proc -> eval -> eval
Execution halted
[/C][/ROW]
[/TABLE]
Source: https://freestatistics.org/blog/index.php?pk=113827&T=0

Globally Unique Identifier (entire table): ba.freestatistics.org/blog/index.php?pk=113827&T=0

As an alternative you can also use a QR Code:

The GUIDs for individual cells are displayed in the table below:

Summary of computational transaction
Raw Input	view raw input (R code)
Raw Output	view raw output of R engine
Computing time	1 seconds
R Server	'Gwilym Jenkins' @ 72.249.127.135
R Engine error message	Error in eval(expr, envir, enclos) : y values must be 0 <= y <= 1 Calls: brglm -> fit.proc -> eval -> eval Execution halted

Parameters (Session):

Parameters (R input):

R code (references can be found in the software module):

library(brglm)
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)
rc <- summary(r)$coeff
try(hm <- hosmerlem(y[1,],r$fitted.values),silent=T)
try(hm,silent=T)
bitmap(file='test0.png')
ra <- roc.analysis(r)
dev.off()
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)
}
bitmap(file='test1.png')
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()
load(file='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='mytable.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='mytable1.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='mytable2.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='mytable3.tab')

Free Statistics

Description of Statistical Computation

Tree of Dependent Computations

Dataset

Tables (Output of Computation)

Figures (Output of Computation)

Input Parameters & R Code