# 5. harjoitus, 8. vko 2002

# Aineisto budworm: Collett, Modellin Binary Data, Chapman & Hall, 1991 & 1999, s. 75


#######
#5.1
#######
ldose <- rep(0:5, 2)
numdead <- c(1, 4, 9, 13, 18, 20, 0, 2, 6, 10, 12, 16)
number <-rep(20,12) # painomuuttuja
y <- numdead/number # vastemuuttuja

# (a)
############

budworm.lg <- glm(y ~ ldose,weight=number, family=binomial)
budworm.lg$coefficients
# (Intercept)       ldose 
#  -2.766087    1.006807 

# (b) Malli M
#############

budworm.lg$coefficients
#  -2.766087    1.006807 
lk<- -2.766087+1.006807*ldose
e<-exp(lk)
pihat<- e/(1+e)

dbinom(numdead, 20, pihat)
# [1] 0.37138363 0.17754498 0.08523514 0.13493922 0.10388903 0.13943161
# [7] 0.29518147 0.23488723 0.18801011 0.14995309 0.03563229 0.07759358
prod(dbinom(numdead, 20, pihat))
# [1] 5.93709e-11 ; uf:n maksimi

log(prod(dbinom(numdead, 20, pihat)))
# [1] -23.54722  ; luf:n maksimi

sum(log(dbinom(numdead, 20, pihat)))
# -23.54722 = L_M ; luf:n maksimi
LM<- -23.54722

# (c) Täysi malli M_f
######################

log(prod(dbinom(numdead, 20, y)))
# [1] -15.0552  ; luf:n maksimi

sum(log(dbinom(numdead, 20, y)))
# -15.0552 = L_M ; luf:n maksimi

Lf<- -15.0552

#######################################

#5.2
#########
# (a)

# Devianssi: D(M) = 2*(L_f - L_M)
DM<- 2*(Lf-LM)
DM
# [1] 16.98404

# (b)
#####
# H_0: beta=0

budworm0 <- glm(y ~ 1,weight=number, family=binomial)
D0<-budworm0$deviance
D0
# [1] 124.8756
DM<-budworm.lg$deviance
DM
# [1] 16.98403
# D0-DM = 2*(LF-L0)-2*(Lf-LM) = 2*(LM-L0)
D0-DM
107.8916
pchisq(107.892, 1,lower.tail = FALSE) # D0-DM noudattaa Khii2-jakaumaa (df=1), kun H_0 tosi
# 2.838307e-25; testisuureen arvoon liittyvä p-arvo
# p-arvo pieni; H0 hylätään; dose siis merkitsevä selittäjä
 
anova(budworm.lg,test="Chisq")
#      Df Deviance Resid. Df Resid. Dev P(>|Chi|)
# NULL                     11    124.876          
# ldose  1  107.892        10     16.984 2.839e-25

#################################################

#5.3
#####
y1<-numdead; y2<- 20-numdead
y1;y2
# [1]  1  4  9 13 18 20  0  2  6 10 12 16
# [1] 19 16 11  7  2  0 20 18 14 10  8  4
y1hat<-20*pihat;y2hat<-20-y1hat
y1hat; y2hat
# [1]  1.183690  2.937611  6.405654 11.264859 15.584449 18.123750  1.183690
# [8]  2.937611  6.405654 11.264859 15.584449 18.123750
# [1] 18.816310 17.062389 13.594346  8.735141  4.415551  1.876250 18.816310
# [8] 17.062389 13.594346  8.735141  4.415551  1.876250

y12<-c(y1,y2)
y12hat<-c(y1hat,y2hat)
ye<- y12/y12hat

yepos <- ye[ye>0]
ypos <- y12[y12>0]
DM<- 2*sum(ypos*log(yepos))
DM
# [1] 16.98403
budworm.lg$deviance
# [1] 16.98403

####################

# 5.4
# Sukupuolen ja annoksen yhdysvaikutus mukaan

sex <- factor(rep(c("M", "F"), c(6, 6)))
budworm.sli <- glm(y ~ sex*ldose,weight=number, family=binomial)
anova(budworm.sli,test="Chisq")

# sex:ldose  1    1.763         8      4.994     0.184
# Käyrien muoto ei eroa kovin merkitsevästi

budworm.sl <- glm(y ~ sex + ldose,weight=number, family=binomial)
budworm.sl$deviance
# [1] 6.757064

budworm.sli$deviance
# [1] 4.993727

# Valinta mallien välillä perustuu devianssien erotukseen
budworm.sl$deviance-budworm.sli$deviance
# [1] 1.763337; ei merkitsevä
# Valitaan yksinkertaisempi malli, jos devianssi ei merkitsevä

# Valitaan malli "budworm.sl"

###################
# Sukupuolten käyrät, yhdysvaikutus

plot(c(1,32), c(0,1), type="n", xlab="dose",
   ylab="prob", log="x")
text(2^ldose, numdead/20, labels=as.character(sex))
ld <- seq(0, 5, 0.1)
lines(2^ld, predict(budworm.sl, data.frame(ldose=ld,
  sex=factor(rep("M", length(ld)), levels=levels(sex))),
  type="response"), col=3)
lines(2^ld, predict(budworm.sl, data.frame(ldose=ld,
  sex=factor(rep("F", length(ld)), levels=levels(sex))),
  type="response"), lty=2, col=2)

#Sukupuolten käyrät, ei yhdysvaikutusta


lines(2^ld, predict(budworm.sli, data.frame(ldose=ld,
  sex=factor(rep("M", length(ld)), levels=levels(sex))),
  type="response"), col=4)
lines(2^ld, predict(budworm.sli, data.frame(ldose=ld,
  sex=factor(rep("F", length(ld)), levels=levels(sex))),
  type="response"), lty=2, col=5)

############################
# 5.7
# Uudet AIDS-tapaukset

# (a)
y <- scan()
0  0  3  0  1  1  1  2  2
4  2  8  0  3  4  5  2  2
2  5  4  3 15 12  7 14 6
10 14  8 19 10  7 20 10 19
 
i <- 1:36
aids.reg <- glm(y ~ i,family=poisson)  # poisson pienellä vaikka nimi
aids.reg
# aids.reg$coefficients

# (b)
plot(i,y)
f <- aids.reg$fitted.values
lines(f)

# Käyrä toisella tavalla
plot(i,y,xlim=c(0,40),ylim=c(0,25))
plot(i,y,main="AIDSiin sairastuneet Englannissa 1983-1985",
xlab="Aika: 1=joulukuu 1982, 36=marraskuu 1985", ylab="AIDSiin sairastuneiden Lkm kuukaudessa")
#####
#plot(c(1,40), c(0,25), type="n", xlab="Aika",ylab="AIDSiin sairastuneet/kk")
# text(i,y, labels=as.character("x"))
lines(predict(aids.reg, data.frame(i), type="response"))

# (c)
#####
aids.reg
#Degrees of Freedom: 35 Total (i.e. Null);  34 Residual
#Null Deviance:      190.2 
#Residual Deviance: 62.36        AIC: 177.7 

anova(aids.reg,test="Chisq")

# Df Deviance Resid. Df Resid. Dev P(>|Chi|)
# NULL                    35     190.17          
# i     1   127.81        34      62.36 1.238e-29

pchisq(62.36, 34 , ncp=0, lower.tail = F)
#[1] 0.002131935

f # pieniä estimoituja frekvenssejä paljon

# (d) Yhdistetään luokkia
######

yy <- scan()
5 5 6 8 7 5 6 5 7 15 12  7 14 6
10 14  8 19 10  7 20 10 19


ii <- c(3.5,8,10.5,12,14,16,18,20,21.5,23:36)
aids1.reg <- glm(yy ~ ii,family=poisson)
aids1.reg  
pchisq(26.155, 21, ncp=0, lower.tail = F)
# [1] 0.2006005 ; yhteensopivuus O.K.