################################################################################################
########## Linear models and testing of several coefficients simultaneously
################################################################################################
setwd("pathway ... /Users/tonje/stk3100/oppgaver2014")     # set work directory where you want R to work from


###Read the data set:
# either you save the data locally on you computer, load package and dataset using the pathway where you saved the data
birthweight = read.table("pathway ... /data/birthweight.txt", header=T)

# or you read the data directly from the webpage:
birthweight = read.table("http://www.uio.no/studier/emner/matnat/math/STK3100/h14/data/birthweight.txt", header=T)

birthweight
#  svlengde vekt sex
#        40 2968   1
#        38 2795   1
#        40 3163   1

#plot the data:
plot(x=birthweight$svlengde, y=birthweight$vekt, col=birthweight$sex)



######Start calculations:
# make the variable sex a factorial vector
birthweight$sex = as.factor(birthweight$sex)


#a
#model number 1:
model1 = lm(vekt~sex+svlengde-1,data=birthweight)
model.matrix(model1)


#add two columns (covariates) with length for boys and girls separately

#model number 2:
birthweight$x4 = birthweight$svlengde*(birthweight$sex==1)
birthweight$x5 = birthweight$svlengde*(birthweight$sex==2)

bw.full <- lm(vekt~sex+x4+x5-1,data=birthweight)
model.matrix(bw.full)

#b
bw.restricted <- lm(vekt~sex-1,data=birthweight)        #beta1=beta2=0 when H0 is true
model.matrix(bw.restricted) 

#Theory: 4.15 in J&H, page 57.
q = 2                                                   #the number of parameters restricted to zero
p = 4                                                   #the number of parameters in the unrestricted model
n = 24                                                  #the number of observations

RSS_restricted = sum(bw.restricted$residuals^2)
RSS_full = sum(bw.full$residuals^2)
est_var_full = summary(bw.full)$sigma^2

Fstat <- (RSS_restricted - RSS_full)/(q*est_var_full)   #is F(q,n-p)-distributed
pvalue <- pf(Fstat,q,n-p,lower.tail=F)                  #Reject H0