# ***** Problem S8 *****
birth<-read.table("http://www.uio.no/studier/emner/matnat/math/STK3100/data/Birth.txt",header=T,row.names=NULL)
#models with different numer of covariates:
mod10 <- glm(children~age+I(age^2) ,family=poisson,data=birth) # canonical link: log
mod11 <- glm(children~age ,family=poisson,data=birth) # canonical link: log
mod12 <- glm(children~1 ,family=poisson,data=birth) # canonical link: log
#In function formula. There it is used to inhibit the interpretation of operators such as "+", "-", "*" and "^"
#as formula operators, so they are used as arithmetical operators. This is interpreted as a symbol by terms.formula.
summary(mod10) #beta_0, beta_1, beta_2
summary(mod11) #beta_0, beta_1
summary(mod12) #beta_0
#a) test H: beta_2=0
# likelihood ratio test
anova10_11<-anova(mod11,mod10) # look under the variable called the "Deviance"
q = 1
Chi2 <- anova10_11[[2,4]]
1-pchisq(Chi2,q) # p-value is large thus fail to reject H0 and beta_2 = 0 is not rejected
# Wald test
b<-mod10$coeff # beta_hat
C1<-cbind(0,0,1) # testing for the last beta only
r = 0
W1<-t(C1%*%b - r)%*%solve((C1)%*%vcov(mod10)%*%t(C1))%*%(C1%*%b - r) # Wald statistic
1-pchisq(W1,q) # p-value is large thus fail to reject H0 and beta_2 = 0 is not rejected
#b) test H: beta_1=beta_2=0
# likelihood ratio test
anova10_12<-anova(mod12,mod10) # look under the variable called the "Deviance"
q = 2
Chi2 <- anova10_12[[2,4]] # look under the variable called the "Deviance"
1-pchisq(Chi2,q) # p-value is low thus we have strong evidence to reject Ho: beta_1=beta_2=0
C2<-cbind(c(0,0),diag(c(1,1))) # testing for the last two betas
r = c(0,0)
W2<-t(C2%*%b-r)%*%solve(C2%*%vcov(mod10)%*%t(C2))%*%(C2%*%b-r) # Wald statistic
1-pchisq(W2,2) # p-value is low thus we have strong evidence to reject Ho: beta_1=beta_2=0