# Low birth weight data
# =============================================
#
# Summary: 
# The data on 189 births were collected at Baystate Medical Center,
# Springfield, Massachusetts, during 1986. The dataset contains an 
# indicator of low infant birth weight as a response and several 
# risk factors associated with low birth weight. The actual birth 
# weight is also included in the dataset. 
#
# Data description: 
# The dataset consists of the following 10 variables:
#Response vaiable:
#  low: indicator of birth weight less than 2.5kg 
#Potential covariates:
#  age: mother's age in years 
#  lwt: mother's weight in pounds at last menstrual period 
#  race: mothers 'race' ("white", "black", "other") 
#  smoke: smoking status during pregnancy 
#  ht: history of hypertension 
#  ui: presence of uterine irritability 
#  ftv: number of physician visits during the first trimester 
#  ptl: number of previous premature labours 
#  bwt: birth weight in grams
#
#======================================================


# Read data:
data="http://www.uio.no/studier/emner/matnat/math/STK3100/data/lowbirthweight.txt"
lowbirthweight=read.table(data,header=T)

# Inspect first lines of data file:
head(lowbirthweight)

# Fit logistic regression:
fit.all=glm(low~age+lwt+smoke+factor(race)+ht,family=binomial,data=lowbirthweight)
summary(fit.all)

# Estimated covariance matrix:
covhat=summary(fit.all)$cov.scaled
print(covhat,3)

# Estimates and standard errors 
# (cf. output from summary-command):
cbind(fit.all$coef, sqrt(diag(covhat)))

# Estimated correlation matrix:
corhat=cov2cor(covhat)
print(corhat,3)

# Test for no effect of age 

# Wald test directly from the summary(fit.all); cf above

# Likelihood ratio statistic from scratch
# First fit model without age:
fit.noage=glm(low~lwt+smoke+factor(race)+ht,family=binomial,data=lowbirthweight)
# Then extract loglikelihoods and compute P-value
L1=logLik(fit.all)
L0=logLik(fit.noage)
DiffL=as.numeric(-2*(L0-L1))
DiffL
Pval=1-pchisq(DiffL,1)
Pval

# Likelihood ratio test using deviance (details about deviance later)
anova(fit.noage,fit.all,test="LRT")

# Score test (Rao's efficient score test):
anova(fit.noage,fit.all,test="Rao")

# Test of all covariates is easier done with the drop1-command:
drop1(fit.all,test="LRT")

# Test of remaining covariates in a model without age
drop1(fit.noage,test="LRT")

# All remaining covariates are significant
# Summary of this model:
summary(fit.noage)