# Birthweight example (cf. lecture chapter 1):
#---------------------

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

# Plot data
plot(firstborn$weeks,firstborn$weight,col=c(rep('blue',28),rep('red',28)),
     pch=c(rep(15,28),rep(16,28)),xlab="Weeks",ylab="Weight")
legend(36,4500,c("Boys","Girls"),col=c("blue","red"),pch=c(15,16),bty="n")

# Fit linear regression model (cf. lecture 1) and plot regression lines 
fit=lm(weight~weeks+factor(sex),data=firstborn)
summary(fit)
abline(fit$coef[1],fit$coef[2],col='blue',lwd=2)
abline(fit$coef[1]+fit$coef[3],fit$coef[2],col='red',lwd=2)

# We can consider different models
# M0: mui=b0
# M1: mui=b0+b1*weeks
# M2; mui=b0+b1*weeks+b2*sex 
# M3: mui=b0+b1*weeks+b2*sex+b3*weeks*sex
# MNULL: mui=0

# The result of M2 is in "fit"
fit.M2 <- fit
#The other models are obtained as
fit.M1=lm(weight~weeks,data=firstborn) #M1
fit.M3=lm(weight~weeks*factor(sex),data=firstborn) #M3
fit.M0=lm(weight~1,data=firstborn) #M0

# ANOVA-tables can then be generated by
anova(fit.M0,fit.M1,fit.M2,fit.M3)

#But to get the F-values for comparing nested models
#from the lecture you have to call anova pairwise: 
#((SSE0-SSE1)/(p1-p0))/(SSE1/(n-p1))
anova(fit.M0,fit.M1)
#Check F-value:
(anova(fit.M0,fit.M1)$Sum[2]/anova(fit.M0,fit.M1)$Df[2])/(anova(fit.M0,fit.M1)$RSS[2]/anova(fit.M0,fit.M1)$Res.Df[2])
#More than one way to Rome
((anova(fit.M0,fit.M1)$RSS[1]-anova(fit.M0,fit.M1)$RSS[2])/((anova(fit.M0,fit.M1)$Res.Df[1])-anova(fit.M0,fit.M1)$Res.Df[2]))/(anova(fit.M0,fit.M1)$RSS[2]/anova(fit.M0,fit.M1)$Res.Df[2])
#And corresponding t-value squared:
7.632^2
#Also interesting, see p-value for 'weeks' in summary for M1:
summary(fit.M1)
  
anova(fit.M1,fit.M2)
anova(fit.M2,fit.M3)

fit.MNULL=lm(weight~-1,data=firstborn) #MNULL
anova(fit.MNULL,fit.M0)
# Just to check
sum(firstborn$weight^2)

# We could also model weeks very flexible by a 1-way ANOVA 
#(weeks are considered as factors/groups)
#M4: mui=b0+sum_(j=37)^42 b_j x_ij, with 
#x_ij=I(baby i born in week j) and b_36=0
factor(firstborn$weeks)
fit1way=lm(weight~factor(weeks),data=firstborn)
summary(fit1way)

# and extend this to a 2-way anova by
fit2way=lm(weight~factor(weeks)+factor(sex),data=firstborn)
fit2wayint=lm(weight~factor(weeks)*factor(sex),data=firstborn)
anova(fit.M0,fit1way)
anova(fit1way,fit2way)
anova(fit2way,fit2wayint)
# It could be worthwhile to check the linearity assumption 
#in M2: mui=b0+b1*weeks+b2*sex
anova(fit.M2,fit2way)