## R-kode for logistisk regresjon, Avsnitt 12.10

x=seq(-3,3,0.1)
plot(x,exp(x)/(1+exp(x)),type="l",ylab="sannsynlighet")
title("Logistisk modell")

# Biller
biller=read.table("https://www.uio.no/studier/emner/matnat/math/STK3100/data/beetle.dat")
names(biller)=c("dose","antall","dode")
biller

plot(biller$dose,biller$dode/biller$antall,xlab="dose",ylab="andel døde",ylim=c(0,1))

billerfit=glm(cbind(dode,antall-dode)~dose,family=binomial,data=biller)
billerfit

arg=seq(1.65,1.9,0.01)
linpr=billerfit$coef[1]+billerfit$coef[2]*arg
prob=exp(linpr)/(1+exp(linpr))
lines(arg,prob)

## Eksempel 12.34 (Challenger-ulykken)

challenger = read.table("https://www.uio.no/studier/emner/matnat/math/STK1110/data/exmp12-34.txt",header=TRUE)
plot(challenger$temp,challenger$fail,xlab="Temperatur (F)",ylab="Feil")

# I praksis hvordan vi gjør logistisk regresjon:
fit=glm(fail~temp,family=binomial,data=challenger)
fit
summary(fit)

# plot(challenger$temp,challenger$fail,xlab="Temperatur (F)",ylab="Feil")
# Legger inn tilpasset modell i plottet
temparg=seq(50,85,1)
linpred=fit$coef[1]+fit$coef[2]*temparg
prob=exp(linpred)/(1+exp(linpred))
lines(temparg,prob)


# Definerer respons og forklaringsvariabel (som temperatur-70 grader F)
y = challenger$fail
x = challenger$temp-70


# Bestemmer maximum likelihood-estimatene - "the hard way" 

# Funksjon som beregner log-likelihooden:
loglik = function(b0,b1)
{
  sum(b0*y+b1*x*y-log(1+exp(b0+b1*x)))
}

# Funksjon som beregner den negatve log-likelihooden:
# (siden kommandoen "optim" minimerer, ikke maksimerer)
negloglik = function(beta,x,y)
{
  beta0 = beta[1]
  beta1 = beta[2]
  
  -sum(beta0*y+beta1*x*y-log(1+exp(beta0+beta1*x)))
}

# Maksimerer log-likelihooden (ved aa minimere den negative log-likelihooden numerisk):
ml.est = optim(c(0,0),negloglik,x=x,y=y)  
beta0.hat = ml.est$par[1]
beta1.hat = ml.est$par[2]
c(beta0.hat,beta1.hat)

# Funksjon som beregner (estimert) sannsynlighet for feil som funksjon av x
px = function(x)
{
  exp(beta0.hat+beta1.hat*x)/(1+exp(beta0.hat+beta1.hat*x))
}

# Plotter stimert sannsynlighet for feil som funksjon av temperatur (x+70) 
xx = seq(-15,10,1)
plot(xx+70,px(xx),type='l',ylim=c(0,1),xlab="Grader Farenheit",ylab="Sannsynlighet for feil")

# Estimerer sannsynligheten for feil ved temperaturen 31 grader F, som var temperaturen da
# Challenger ble skutt opp (her boer en vaere obs paa at dette punktet ligger langt utafor
# omraadet med data brukt for aa estimere modellen)
px(31-70)

## Eksempel 12.34 (forts.)

# Informasjonsmatrise 
inf.matrise = matrix(0,2,2)
lin = beta0.hat+beta1.hat*x
inf.matrise[1,1] = sum(exp(lin)/(1+exp(lin))^2)
inf.matrise[2,2] = sum(x^2*exp(lin)/(1+exp(lin))^2)
inf.matrise[1,2] = sum(x*exp(lin)/(1+exp(lin))^2)
inf.matrise[2,1] = inf.matrise[1,2]

# Estimert kovariansmatrise og standrardfeil
cov.matrise = solve(inf.matrise)
var.est = diag(cov.matrise)
se.est = sqrt(var.est)
se.est

## Eksempel 12.34 (forts.)

# Vi kan bruke likelihood ratio-testen til aa teste 
# H_0: beta1=0 mot alternativet H_a: beta1 != 0

# Under H_0 avhenger ikke sammsynligheten av forklarinsvariabelen x
# Maksimum likelihood-estimat for den sannsynligheten under H_0 er
p0.hat = sum(y)/length(y)

# Maksimum likelihood-estimatet for den tilsvarende verdi av beta_0 er
beta0.null = log(p0.hat/(1-p0.hat))

# Beregner -2*ln(LR) og P-verdien for likelihood ratio-testen
minus.2lnLR = 2*(loglik(beta0.hat,beta1.hat)-loglik(beta0.null,0))
1-pchisq(minus.2lnLR,1)