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

Sun, 19 Dec 2010 10:12:57 +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/19/t12927565245kkaa0qda54tzj7.htm/, Retrieved Sun, 05 May 2024 04:51:37 +0000

Statistical Computations at FreeStatistics.org, Office for Research Development and Education, URL https://freestatistics.org/blog/index.php?pk=112280, Retrieved Sun, 05 May 2024 04:51:37 +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

147

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] [b7765ad69c3ab250b1ef04c2ab1247ec] [Current]
-    D      [Bias-Reduced Logistic Regression] [] [2010-12-21 18:56:29] [20c5a34fea7ed3b9b27ff444f2eb4dfe] 
-   P       [Bias-Reduced Logistic Regression] [] [2010-12-21 21:31:33] [07fa8844ca5618cd0482008937d9acea] 

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	'Sir Ronald Aylmer Fisher' @ 193.190.124.24
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 & 'Sir Ronald Aylmer Fisher' @ 193.190.124.24 \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=112280&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]'Sir Ronald Aylmer Fisher' @ 193.190.124.24[/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=112280&T=0

Globally Unique Identifier (entire table): ba.freestatistics.org/blog/index.php?pk=112280&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	'Sir Ronald Aylmer Fisher' @ 193.190.124.24
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