###################### ESIMERKKI 6.1 ##############################

beetle <- read.table("http://users.jyu.fi/~junyblom/beetle.dat", header=TRUE) 
attach(beetle)

beetle.lgt <- glm(killed/total~dose, family=binomial("logit"), weights=total)
coef(beetle.lgt)

# Logitin käänteisfunktio

 invlogit <- function(x,b){
 # x = annos, b = kerroinvektori
 z <- 1/(1 + exp(-b[1] - b[2]*x))
 names(z) <- "Todennäköisyys"
 z
                     }


 oldpar <- par(no.readonly=TRUE)
 par(mfrow=c(1,2), oma=c(1,1,1,1))

 par(mar = c(3,4,0,0)) 

 plot(dose,killed/total, xlab="", ylab=" ")
 curve(invlogit(x,coef(beetle.lgt)), from=1.69,to=1.89, 
 add=TRUE)

 mtext("Annos (10-kantainen logaritmi)", font=1, side=1, line=2.5) 
 mtext("Todennäköisyys", font=1, side=2, line=2.5)

 par(mar=c(3,4,0,0))

 plot(10^dose,killed/total, xlab="", ylab=" ")
 curve(invlogit(log10(x),coef(beetle.lgt)), from=49,to=77, 
 add=TRUE)

 mtext("Annos (mg/l)", font=1, side=1, line=2.5) 
 mtext("Todennäköisyys", font=1, side=2, line=2.5)

 par(oldpar)

###################### ESIMERKKI 6.2 ##################

# Komplementaarisessa log-log -tapauksessa kertoimet on etsitty pienimmän 
# neliösumman avulla:

 x <- seq(.01,.99,len=99)
 out.ll <- lm(log(x/(1-x))~log(-log(1-x)))

 beetle <- read.table("http://users.jyu.fi/~junyblom/beetle.dat", 
 header=T)

 attach(beetle)

 beetle.lgt <- glm(killed/total~dose, family=binomial, weights=total)
 b.lgt <- coef(beetle.lgt)
 summary(beetle.lgt)
1 - pchisq(11.232, df=6)

 beetle.pbt <- glm(killed/total~dose, weights=total, 
 family=binomial(link="probit"))
 b.pbt <- coef(beetle.pbt)
 summary(beetle.pbt)
1 - pchisq(10.12, df=6)

 beetle.cll <- glm(killed/total~dose, weights=total, 
 family=binomial(link="cloglog"))
 b.cll <- coef(beetle.cll)
 summary(beetle.cll)
1 - pchisq(3.4464, df=6)

oldpar <- par(no.readonly=TRUE)
 par(mfrow=c(1,2), oma=c(1,1,1,1))

 par(mar = c(3,3,0,0)) 

 curve(log(x/(1-x)), from=0, to=1, xlab=" ", ylab=" ")
 curve(qnorm(x, sd=pi/sqrt(3)), from=0, to=1, add=T, lty=2)
 curve(coef(out.ll)[2]*log(-log(1-x))+coef(out.ll)[1], from=0, to=1, add=T, lty=3)

 mtext("Todennäköisyys", font=1, side=1, line=2.5) 
 mtext("Arvo", font=1, side=2, line=2.5)

 par(mar=c(3,4,0,0))

 plot(10^dose,100*killed/total, ylim=c(0,100), xlab="", 
 ylab=" ")
 
 curve(100*exp(b.lgt[1]+b.lgt[2]*log10(x))/(1+exp(b.lgt[1]+b.lgt[2]*log10(x))), 
 from=49, to=77, add=T)

 curve(100*pnorm(b.pbt[1]+b.pbt[2]*log10(x)), from=49, to=77, 
 add=T, lty=2, col=2)
 
 curve(100*(1-exp(-exp(b.cll[1]+b.cll[2]*log10(x)))), from=49, to=77, 
 add=T, lty=3, col=3)

 mtext("Annos (mg/l)", font=1, side=1, line=2.5) 
 mtext("Kuolleita (%)", font=1, side=2, line=2.5)

 par(oldpar)

### TAULUKKO 6.1 ###############

M.l <- cbind(killed,total*beetle.lgt$fit,total-killed,
total*(1-beetle.lgt$fit))

M.p <- cbind(killed,total*beetle.pbt$fit,total-killed,
total*(1-beetle.pbt$fit))

M.c <- cbind(killed,total*beetle.cll$fit,total-killed,
total*(1-beetle.cll$fit))

M.k <- cbind(M.l[,1:2],M.p[,2],M.c[,2])
colnames(M.k) <- c("Havaitut","Logit","Probit","Kompl.log-log")
round(M.k,1)

M.e <- cbind(M.l[,3:4],M.p[,4],M.c[,4])
colnames(M.e) <- c("Havaitut","Logit","Probit","Kompl.log-log")
round(M.e,1)


######################## IWLS estimointi ########################

eps <- 1
n <- total
y <- killed
X <- cbind(1,dose)
b <- rep(0,ncol(X))
names(b) <- c("Intercept","dose")
loglik <- c(y%*%X%*%b) - n%*%log(1+exp(X%*%b))
print(loglik)
while(eps>1e-10){
prob <- exp(X%*%b)/(1 + exp(X%*%b))
w <- n*prob*(1-prob)
u <- (y - n*prob)/w
lmw <- lm(u ~ X -1, weights=w)
eps <- max(abs(d <- coef(lmw)))
b <- b + d
loglik <- c(y%*%X%*%b) - n%*%log(1+exp(X%*%b))
a <- c(b, loglik)
names(a)[3] <- "loglik"
print(a)
}
summary(lmw)$cov.unscaled
display(beetle.lgt)
cov.matrix(beetle.lgt)

################ ESIMERKKI 6.4 ###################################

## Funktioita
## Luottamusväli-funktio lineaariselle prediktorille

raja <- function(x,obj,level=.95){
a <- .5*(1-level)
X <- cbind(1,x)
cov <- summary(obj)$cov.scaled
s <- sqrt(diag(X%*%cov%*%t(X)))
b <- qnorm(a)*s
b
                                  }

## Funktio: Bootstrap-rajat lineaariselle prediktorille

raja.boot <- function(x,b.boot,level=.95){
a <- .5*(1-level)
X <- cbind(1,x)
fits <- X%*%b.boot
limits <- apply(fits, 1, quantile, probs= c(a,1-a))
                                  }

## Logitin käänteisfunktio = logistisen jakauman kertymäfunktio
plogis <- function(x){
z <- exp(x)/(1+exp(x))
z
                     }

logit <- function(p){
z <- log(p/(1-p))
z
                    }


## Luottamusväli todennäköisyydelle, teoria

plot(10^dose,killed/total, xlab="Annos (mg/l)", ylab="Todennäköisyys")
curve(plogis(b.lgt[1]+b.lgt[2]*log10(x)), from=49, to=77, add=T)
curve(plogis(b.lgt[1]+b.lgt[2]*log10(x) + raja(log10(x),beetle.lgt)),  
from=49, to=77, lty=2, add=T)
curve(plogis(b.lgt[1]+b.lgt[2]*log10(x) - raja(log10(x),beetle.lgt)),  
from=49, to=77, lty=2, add=T)

# bootstrapilla

N <- 9999
b <- matrix(0,2,N)
for(i in 1:N){
y <- rbinom(8, prob=fitted(beetle.lgt), size=total)
b[,i] <- coef(glm(y/total ~ dose, weights=total, family=binomial))
}

confint(beetle.lgt)
Q <- t(apply(b,1,quantile, c(.025,.975)))
rownames(Q) <- names(coef(beetle.lgt))
Q
curve(plogis(raja.boot(log10(x), b.boot=b))[1,], from=49, to=77, 
add=TRUE, lty=3)


############## ESIMERKKI 6.5 #################################

heart <- read.table("http://users.jyu.fi/~junyblom/heart.dat", header=T)

attach(heart)
sc <- factor(sc)
smoke <- factor(smoke)
alc <- factor(alc)
hist <- factor(hist)
age <- factor(age)


heart1.glm <- glm(trouble/total ~ smoke + alc + age + hist + sc,
weights=total, family=binomial)

summary(heart1.glm)

# Ristitulosuhteet l. OR:t

 round(exp(coef(heart1.glm)[-1]),1)

# Ristitulosuhteiden luottamusvälit

  t(round(exp(confint(heart1.glm)[-1,]),1))

# Entisten ja nykyisten tupakoijien vertailu vertailu

 d <- coef(heart1.glm)["smoke3"]-coef(heart1.glm)["smoke2"]
 C <- summary(heart1.glm)$cov.unscaled[c("smoke2","smoke3"),
 c("smoke2","smoke3")]
 s.e. <- sqrt(C[1,1]+C[2,2]- 2*C[1,2])
 c(d,s.e.)
 exp(d + c(-1.96,1.96)*s.e.)

################### ESIMERKKI 6.6 #########################

B <- 9999 # bootstrap-toistoja
dev.boot <- numeric(B)              # bootstrap devianssit tähän
coef.boot <- matrix(0,B,8)          # kertoimet tähän
colnames(coef.boot) <- 
names(coef(heart1.glm))
k <- 72                             # luokkien lukumäärä saturoidussa mallissa
phat <- fitted(heart1.glm)
for(i in 1:B){
y <- rbinom(k,prob=phat,size=total)                  # bootstrap-otos
out <- glm(y/total ~ smoke + alc + age + hist + sc,
weights=total, family=binomial)                      # mallin sovitus
dev.boot[i] <- out$dev
coef.boot[i,] <- coef(out)
}

#Bootstrap p-arvo, vertailu asymptoottisiin tuloksiin

mean(dev.boot >= heart1.glm$dev) #devianssin bootstrap p-arvo

hist(dev.boot, breaks=100, prob=TRUE)

############### JÄÄNNÖKSET JA DIAGNOSTIIKKA ####################

## Kuva 6.4

oldpar <- par(no.readonly=TRUE)
 par(mfrow=c(1,2), oma=c(1,1,1,1))

 par(mar = c(3,3,0,0)) 

 hav.tn <- trouble/total
 plot(heart1.glm$lin, log(hav.tn/(1-hav.tn)), xlab=" ",
 ylab=" ")
 abline(0,1)


 mtext("Lineaarinen prediktori", font=1, side=1, line=2.5) 
 mtext("Empiirinen logit", font=1, side=2, line=2.5)

 par(mar=c(3,4,0,0))

 plot(heart1.glm$lin,(trouble/total),
 xlab=" ", ylab=" ")
 lines(sort(heart1.glm$lin),heart1.glm$fit[order(heart1.glm$lin)])


 mtext("Lineaarinen prediktori", font=1, side=1, line=2.5) 
 mtext("Todennäköisyys", font=1, side=2, line=2.5)

 par(oldpar)

## Kuva 6.5

oldpar <- par(no.readonly=TRUE)
 par(mfrow=c(1,2), oma=c(1,1,1,1))

 par(mar = c(3,3,0,0)) 

 phat <- heart1.glm$fit
 odfr <- total*phat
 r <- (trouble - odfr)/sqrt(odfr*(1-phat))
 qqnorm(r, main=" ", ylab=" ", xlab=" "); abline(0,1)

 mtext("Teoreettiset kvantiilit", font=1, side=1, line=2.5) 
 mtext("Pearsonin jäännökset", font=1, side=2, line=2.5)

 par(mar=c(3,4,0,0))

 d <- resid(heart1.glm)
 qqnorm(d, main=" ", ylab=" ", xlab=" "); abline(0,1)

 mtext("Teoreettiset kvantiilit", font=1, side=1, line=2.5) 
 mtext("Devianssijäännökset", font=1, side=2, line=2.5)

 par(oldpar)