# Binäärinen data: Alhainen syntymäpaino vastasyntyneillä (Ripley, 1999, MAS with S-Plus)

library(mass)
options(contrasts=c("contr.treatment", "contr.poly"))
data(birthwt)
attach(birthwt)
race <- factor(race, labels=c("white", "black", "other"))
table(ptl) 	# aikaisemmat ennenaikaisuuteen liittyvät ongelmat 
#  0   1   2   3 
# 159  24   5   1 

ptd <- factor(ptl > 0)
table(ptd)
# ptd
# FALSE  TRUE 
#  159    30 

table(ftv)	# lääkärissäkäynnit raskauden 1. kolmanneksella
#ftv
#  0   1   2   3   4   6 
# 100  47  30   7   4   1 
ftv <- factor(ftv)
levels(ftv)[-(1:2)] <- "2+"
table(ftv)  # 
#ftv
#  0   1  2+ 
#100  47  42

# Esimerkki
#####
z <- c(0,1,2,3,4)
z <- factor(z)
levels(z)[1:2]
#[1] "0" "1"
levels(z)[-(1:2)]
#[1] "2" "3" "4"
###

bwt <- data.frame(low=factor(low), age, lwt, race,
    smoke=(smoke>0), ptd, ht=(ht>0), ui=(ui>0), ftv)
detach(); rm(race, ptd, ftv)

birthwt.glm <- glm(low ~ ., family=binomial, data=bwt)
summary(birthwt.glm, cor=F)
###########################
# Tehtävä 1.20
###########################
Call:
glm(formula = low ~ ., family = binomial, data = bwt)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.7038  -0.8068  -0.5009   0.8836   2.2151  

Coefficients:
             Estimate Std. Error z value Pr(>|z|)   
(Intercept)  0.822706   1.240174   0.663  0.50709   
age         -0.037220   0.038530  -0.966  0.33404   
lwt         -0.015651   0.007048  -2.221  0.02637 * 
raceblack    1.192231   0.534428   2.231  0.02569 * 
raceother    0.740513   0.459769   1.611  0.10726   
smokeTRUE    0.755374   0.423246   1.785  0.07431 . 
ptdTRUE      1.343654   0.479409   2.803  0.00507 **
htTRUE       1.912974   0.718586   2.662  0.00776 **
uiTRUE       0.680162   0.463464   1.468  0.14222   
ftv1        -0.436331   0.477792  -0.913  0.36112   
ftv2+        0.178939   0.455227   0.393  0.69426   
---
Signif. codes:  0  `***'  0.001  `**'  0.01  `*'  0.05  `.'  0.1  ` '  1 

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 234.67  on 188  degrees of freedom
Residual deviance: 195.48  on 178  degrees of freedom
AIC: 217.48
###############
# (a) Muuttujalla pdt (aikaisemmat raskauteen liittyvät ongelmat) on suurin z-arvo. 
#  z-arvon mielessä paras selittäjä. ht (kohonut verenpaine) lähes yhtä hyvä.
###############
# AIC = -2*(logaritmoidun uf:n maksimi) + 2* (mallin parametrien lkm)
# uf = uskottavuusfunktio
###############
# (b) update

###############
birthwt.next <- update(birthwt.glm, . ~ . - ptd) 
birthwt.next
birthwt.next <- update(birthwt.glm, . ~ . - lwt)
birthwt.next
birthwt.next <- update(birthwt.glm, . ~ . - race) 
birthwt.next
birthwt.next <- update(birthwt.glm, . ~ . - ht) 
birthwt.next
birthwt.next <- update(birthwt.glm, . ~ . - ftv) 
birthwt.next
birthwt.next <- update(birthwt.glm, . ~ . - age) 
birthwt.next
birthwt.next <- update(birthwt.glm, . ~ . - smoke) 
birthwt.next
birthwt.next <- update(birthwt.glm, . ~ . - ui) 
birthwt.next
##################
# age:n pudottaminen mallista huonontaa vähiten mallin sopivuutta aineistoon
# (res. dev. 196.4).
# Siis lähinnä mallia birthwt.glm oleva yhtä selitt. pienempi malli on sellainen,
# josta age pudotettu pois
# - ftv:n (lääkärissäkäynnit) pudottaminen johtaa lähes yhtä hyvään (res. dev. 196.8)
#   pienempään malliin.
##################
# (c) 
# - Suurin AIC (223.6) on mallilla, josta ptd on pois: ptd:n poistaminen mallista
#   kasvattaa eniten AIC:tä. Seuraavaksi ht (222.9) eli korkea verenpaine. 
# - Ainoa mallia birthwt.glm yhtä selittäjää pienempi malli, jolla on pienempi AIC (214.8) kuin 
# birthwt.glm, saadaan poistamalla ftv (lääkärissäkäyntien lkm).
##################
# Tehtävä 1.21
############# 
# (a) dropterm
#############
dropterm(birthwt.glm, test= "Chisq")
########
Single term deletions

Model:
low ~ age + lwt + race + smoke + ptd + ht + ui + ftv
       Df Deviance     AIC     LRT  Pr(Chi)   
<none>     195.476 217.476                    
age     1  196.417 216.417   0.942 0.331796   
lwt     1  200.949 220.949   5.474 0.019302 * 
race    2  201.227 219.227   5.751 0.056380 . 
smoke   1  198.674 218.674   3.198 0.073718 . 
ptd     1  203.584 223.584   8.109 0.004406 **
ht      1  202.934 222.934   7.458 0.006314 **
ui      1  197.585 217.585   2.110 0.146342   
ftv     2  196.834 214.834   1.358 0.507077   
##########---
# ftv huonoin selittäjä (p-arvo), sitten age
# - jätetään ftv pois
####
# (b) addterm
####
# alkumalli
#
birthwt.start <- glm(low ~ptd,family=binomial, data=bwt) 
birthwt.start
birthwt.start <- glm(low ~lwt,family=binomial, data=bwt) 
birthwt.start
birthwt.start <- glm(low ~race,family=binomial, data=bwt) 
birthwt.start
birthwt.start <- glm(low ~ht,family=binomial, data=bwt) 
birthwt.start
birthwt.start <- glm(low ~ftv,family=binomial, data=bwt) 
birthwt.start
birthwt.start <- glm(low ~age,family=binomial, data=bwt) 
birthwt.start
birthwt.start <- glm(low ~smoke,family=binomial, data=bwt) 
birthwt.start
birthwt.start <- glm(low ~ui,family=binomial, data=bwt) 
birthwt.start
######

birthwt.start <- glm(low ~ ptd, family=binomial, data=bwt)
birthwt.start
addterm(birthwt.start,  ~ . + age + lwt + race + smoke + ptd + ht + ui +ftv,test= "Chisq")
birthwt.start <- glm(low ~ ptd+race, family=binomial, data=bwt)
birthwt.start
addterm(birthwt.start,  ~ . + age + lwt + race + smoke + ptd + ht + ui +ftv,test= "Chisq")
birthwt.start <- glm(low ~ ptd+race+smoke, family=binomial, data=bwt)
birthwt.start
addterm(birthwt.start,  ~ . + age + lwt + race + smoke + ptd + ht + ui +ftv,test= "Chisq")
birthwt.start <- glm(low ~ ptd+race+smoke+lwt, family=binomial, data=bwt)
birthwt.start
addterm(birthwt.start,  ~ . + age + lwt + race + smoke + ptd + ht + ui +ftv,test= "Chisq")
birthwt.start <- glm(low ~ ptd+race+smoke+lwt+ht, family=binomial, data=bwt)
birthwt.start
addterm(birthwt.start,  ~ . + age + lwt + race + smoke + ptd + ht + ui +ftv,test= "Chisq")
birthwt.start <- glm(low ~ ptd+race+smoke+lwt+ht+ui, family=binomial, data=bwt)
birthwt.start
addterm(birthwt.start,  ~ . + age + lwt + race + smoke + ptd + ht + ui +ftv,test= "Chisq")
birthwt.start <- glm(low ~ ptd+race+smoke+lwt+ht+ui+ftv, family=binomial, data=bwt)
birthwt.start
addterm(birthwt.start,  ~ . + age + lwt + race + smoke + ptd + ht + ui +ftv)
############
# Tehtävä 1.22
############# 
# (a) dropterm
#############
dropterm(birthwt.glm)
########
Single term deletions

Model:
low ~ age + lwt + race + smoke + ptd + ht + ui + ftv
       Df Deviance    AIC
<none>      195.48 217.48
age     1   196.42 216.42
lwt     1   200.95 220.95
race    2   201.23 219.23
smoke   1   198.67 218.67
ptd     1   203.58 223.58
ht      1   202.93 222.93
ui      1   197.59 217.59
ftv     2   196.83 214.83
##########---
# - jätetään ftv pois
# 
# (b) addterm

birthwt.start <- glm(low ~ ptd, family=binomial, data=bwt)
birthwt.start
addterm(birthwt.start,  ~ . + age + lwt + race + smoke + ptd + ht + ui +ftv)
birthwt.start <- glm(low ~ ptd+age, family=binomial, data=bwt)
birthwt.start
addterm(birthwt.start,  ~ . + age + lwt + race + smoke + ptd + ht + ui +ftv)
birthwt.start <- glm(low ~ ptd+age+ht, family=binomial, data=bwt)
birthwt.start
addterm(birthwt.start,  ~ . + age + lwt + race + smoke + ptd + ht + ui +ftv)
birthwt.start <- glm(low ~ ptd+age+ht+lwt, family=binomial, data=bwt)
birthwt.start
addterm(birthwt.start,  ~ . + age + lwt + race + smoke + ptd + ht + ui +ftv)
birthwt.start <- glm(low ~ ptd+age+ht+lwt+ui, family=binomial, data=bwt)
birthwt.start
addterm(birthwt.start,  ~ . + age + lwt + race + smoke + ptd + ht + ui +ftv)
birthwt.start <- glm(low ~ ptd+age+ht+lwt+ui+smoke, family=binomial, data=bwt)
birthwt.start
addterm(birthwt.start,  ~ . + age + lwt + race + smoke + ptd + ht + ui +ftv)
birthwt.start <- glm(low ~ ptd+race+smoke+lwt+ht+ui+ftv, family=binomial, data=bwt)
birthwt.start
addterm(birthwt.start,  ~ . + age + lwt + race + smoke + ptd + ht + ui +ftv)
############
(b) # AIC = -2*(logaritmoidun uf:n maksimi) + 2* (mallin parametrien lkm)
# uf = uskottavuusfunktio

############
# Tehtävä 1.23
############
# (a)
#####
#stepAIC(object, scope, scale, direction=c("both", "backward", "forward"), 
#             trace=1, keep=NULL, steps=1000, use.start=FALSE, k=2, ...)
#####
birthwt.step <- stepAIC(birthwt.glm, trace=F)
birthwt.step$anova

Stepwise Model Path 
Analysis of Deviance Table

Initial Model:
low ~ age + lwt + race + smoke + ptd + ht + ui + ftv

Final Model:
low ~ lwt + race + smoke + ptd + ht + ui


   Step Df Deviance Resid. Df Resid. Dev      AIC
1       NA       NA       178   195.4755 217.4755
2 - ftv  2 1.358185       180   196.8337 214.8337
3 - age  1 1.017866       181   197.8516 213.8516

####
# (b)
####
birthwt.step2 <- stepAIC(birthwt.glm, ~ .^2 + I(scale(age)^2)
   + I(scale(lwt)^2), trace = F)
birthwt.step2$anova
###################
Stepwise Model Path 
Analysis of Deviance Table

Initial Model:
low ~ age + lwt + race + smoke + ptd + ht + ui + ftv

Final Model:
low ~ age + lwt + smoke + ptd + ht + ui + ftv + age:ftv + smoke:ui


        Step Df  Deviance Resid. Df Resid. Dev      AIC
1            NA        NA       178   195.4755 217.4755
2  + age:ftv  2 12.474897       176   183.0006 209.0006
3 + smoke:ui  1  3.056805       175   179.9438 207.9438
4     - race  2  3.129586       177   183.0734 207.0734
summary(birthwt.step2, cor=F)$coef
                    Estimate  Std. Error    z value     Pr(>|z|)
(Intercept)      -0.58238913 1.420834491 -0.4098923 0.6818849477
age               0.07553844 0.053944995  1.4002864 0.1614275844
lwt              -0.02037234 0.007487695 -2.7207766 0.0065128768
smokeTRUE         0.78004747 0.420043213  1.8570648 0.0633019332
ptdTRUE           1.56030401 0.496626115  3.1418082 0.0016790798
htTRUE            2.06567991 0.748329570  2.7603879 0.0057732762
uiTRUE            1.81849631 0.666669703  2.7277320 0.0063771389
ftv1              2.92106773 2.284093458  1.2788740 0.2009414315
ftv2+             9.24445985 2.650495048  3.4878239 0.0004869688
age.ftv1         -0.16182328 0.096736157 -1.6728314 0.0943604977
age.ftv2+        -0.41101103 0.118553331 -3.4668872 0.0005265227
smokeTRUE.uiTRUE -1.91664380 0.972365601 -1.9711144 0.0487107960
#################
# Tehtävä 1.24
#################

table(bwt$low, predict(birthwt.step2) > 0)
   
    FALSE TRUE
  0   116   14
  1    28   31

# Malli birthwt.step2 on ennusteen mielessä hieman parempi

table(bwt$low, predict(birthwt.step) > 0)
    FALSE TRUE
  0   117   13
  1    33   26

table(bwt$low, predict(birthwt.glm) > 0)
 
FALSE TRUE
  0   116   14
  1    37   22

# Onko vääriin ennustuksiin suhtauduttava symmetrisesti?

cbind(bwt$low,birthwt.step$fitted.values)

######################################
#1.25
###################
birthwt.step <- stepAIC(birthwt.glm, trace=F)
birthwt.step$anova

#Final Model:
#low ~ lwt + race + smoke + ptd + ht + ui

n <- 189
birthwt.step <- stepAIC(birthwt.glm, trace=F, k = log(n))
birthwt.step$anova

#Final Model:
#low ~ lwt + ptd + ht

birthwt.step2 <- stepAIC(birthwt.glm, ~ .^2 + I(scale(age)^2)
   + I(scale(lwt)^2), trace = F)
birthwt.step2$anova

Final Model:
low ~ age + lwt + smoke + ptd + ht + ui + ftv + age:ftv + smoke:ui

birthwt.step2 <- stepAIC(birthwt.glm, ~ .^2 + I(scale(age)^2)
   + I(scale(lwt)^2), trace = F, k = log(n))
birthwt.step2$anova

# Final Model:
# low ~ lwt + ptd + ht

#########
# Tehtävä 1.26
#########
table(bwt$low, predict(birthwt.step2) > 0)  # BIC: low ~ lwt + ptd + ht

#     FALSE TRUE
#  0   120   10
#  1    39   20

table(bwt$low, predict(birthwt.step2) > 0)  # AIC: 
# low ~ age + lwt + smoke + ptd + ht + ui + ftv + age:ftv + smoke:ui
  
#    FALSE TRUE
#  0   116   14
#  1    28   31
####################