# Auswertung fuer Beispiel 13.12:


ds <- read.table("Buerk.txt",header=T)
str(ds)
Y <- ds$mortal
X <- ds[,1:21]


formel <- "Y ~ X[,1]"
for (j in 2:21)
{
		formel <- paste(formel," + X[,", encodeString(j), "]")
}
res <- glm(formula = formel, family=binomial())

summary(res)

th.hat <- res$coefficients
F.hat <- cbind(1,as.matrix(X)) %*% th.hat
# Log-Likelihood für spaetere Berechnungen:
LL <- (- res$deviance)


tmp <- sort(F.hat,index.return=TRUE)
F.hat <- tmp$x
Y <- Y[tmp$ix]
Y.hat.iso <- IsoMeans(Y)

# Ausduennen fuer Teppichfransenplot:
par(mai=c(0.5,0.5,0.1,0.1),cex=1.2)
F.hat.1 <- F.hat[Y == 1]
F.hat.1l <- F.hat.1[1:20]
F.hat.1r <- F.hat.1[(length(F.hat.1)-19):length(F.hat.1)]
F.hat.1m <- sample(F.hat.1[21:(length(F.hat.1)-20)], size=200)
F.hat.1 <- c(F.hat.1l,F.hat.1m,F.hat.1r)
F.hat.0 <- F.hat[Y == 0]
F.hat.0l <- F.hat.0[1:20]
F.hat.0r <- F.hat.0[(length(F.hat.0)-19):length(F.hat.0)]
F.hat.0m <- sample(F.hat.0[21:(length(F.hat.0)-20)], size=600)
F.hat.0 <- c(F.hat.0l,F.hat.0m,F.hat.0r)
plot(F.hat.1,rep(1.03,240),pch="|",xlab="",ylab="",xlim=c(F.hat[1],6),ylim=c(-0.005,1.005))
points(F.hat.0,rep(-0.03,640),pch="|")
abline(h=0)
abline(h=1)
zz <- seq(F.hat[1],6,length.out=500)
lines(zz,1/(exp(-zz)+1),col="blue",lwd=4)
lines(F.hat,Y.hat.iso,col="black",lwd=2)


N0 <- sum(1 - Y)
N1 <- sum(Y)
n <- N0 + N1

Sens <- rep(0,n+1)
Spez <- rep(0,n+1)
Sens[1] <- N1
Sens[2:(n+1)] <- N1 - cumsum(Y)
Sens <- Sens/N1
Spez[1] <- 0
Spez[2:(n+1)] <- cumsum(1 - Y)
Spez <- Spez/N0

par(mai=c(1,1,0.14,0.14),cex=1.2)
plot(1 - Spez, Sens, type="l",lwd=3, xlim=c(0.037,0.963),ylim=c(0.037,0.963),xlab="1 - Spezifität",ylab="Sensitivität")
abline(a=0,b=1,lty=3)


# P-Werte fuer die 21 Kovariablen
# nach der Profil-Likelihood-Methode:

PV <- rep(1,21)
for (j in 1:21)
{
	if (j == 1)
	{
		formel <- "Y ~ X[,2]"
		for (k in 3:21)
		{
			formel <- paste(formel," + X[,", encodeString(k), "]")
		}
	}
	if (j > 1)
	{
		formel <- "Y ~ X[,1]"
		for (k in 2:21)
		{
			if (k != j)
			{
				formel <- paste(formel," + X[,", encodeString(k), "]")
			}
		}
	}
	
	res0 <- glm(formula = formel, family=binomial())
	LL0 <- (- res0$deviance)
	
	PV[j] <- 1 - pchisq(2*LL - 2*LL0, 1)
}

PVadj <- p.adjust(PV,method="holm")
cbind(PV,PVadj)