#### Exercise 6.3 Female mortality

######### Models
mortality <- read.table("mortality.csv",header=T,sep=",")
mortality <- mortality[,-c(4,6,8,10,11)]
mortality <- na.omit(mortality)
attach(mortality)
library(MASS)


model2 <- glm.nb(Female_death ~ factor(Age)+factor(Year)+
	offset(L_female_exp))
summary(model2,correlation=F)


######### choose m for age and year jointly
AIC.calc.both<-function(age1,age2,year1,year2)
{
AIC <- matrix(0,nrow=age2-age1+1,ncol=year2-year1+1)
for(i in age1:age2)
	{for(j in year1:year2)
		{
		print(c(i,j))
		model.ij <- glm.nb(Female_death ~ poly(Age,i)+poly(Year,j)+offset(L_female_exp))
		AIC[i-age1+1,j-year1+1] <- -model.ij$twologlik + 2*(i+j)

		}
	
	}
AIC
}

AIC.both <- AIC.calc.both(20,30,2,8)

> AIC.both
          [,1]     [,2]     [,3]     [,4]     [,5]     [,6]     [,7]
 [1,] 50219.93 50220.64 50221.64 50215.08 50212.83 50214.24 50214.31
 [2,] 50215.99 50216.70 50217.69 50211.13 50208.87 50210.28 50210.35
 [3,] 50208.86 50209.57 50210.57 50204.01 50201.75 50203.15 50203.23
 [4,] 50194.84 50195.57 50196.57 50190.01 50187.75 50189.14 50189.22
 [5,] 50190.13 50190.86 50191.87 50185.31 50183.05 50184.45 50184.53
 [6,] 50185.67 50186.40 50187.41 50180.86 50178.60 50180.00 50180.08
 [7,] 50184.46 50185.20 50186.20 50179.64 50177.39 50178.78 50178.86
 [8,] 50185.93 50186.69 50187.73 50181.25 50179.05 50180.46 50180.51
 [9,] 50183.41 50184.28 50185.40 50179.20 50177.19 50178.66 50178.60
[10,] 50184.66 50185.55 50186.69 50180.54 50178.56 50180.05 50179.96
[11,] 50187.25 50188.13 50189.28 50183.13 50181.15 50182.64 50182.55

> min(AIC.both)
[1] 50177.19   #### occurs at i=9 and j=5 i.e. p=28 and q=6

model.final <- glm.nb(Female_death ~ poly(Age,28)+poly(Year,6)+
	offset(L_female_exp))
summary(model.final,correlation=F)


-model2$twologlik+2*(109+54)      ### AIC for fully parametrized model
[1] 50275.45
 

## check final model:
-model.final$twologlik+2*(28+6)      ### AIC for final model
[1] 50177.19


##### 3d plots
mortality <- read.table("mortality.csv",header=T,sep=",")
mortality <- mortality[,-c(3,5,7,9,11)]
mortality <- na.omit(mortality)
attach(mortality)
library(MASS)
library(lattice)
model.final <- glm.nb(Female_death ~ poly(Age,28)+poly(Year,6)+
	offset(L_female_exp))


postscript("3dplot.female.1.ps")
wireframe(log(q_female)~Year*Age,data=mortality,  scales = list(arrows = FALSE), drape = TRUE, 
colorkey = FALSE, screen = list(z = 30, x = -60),xlab="Year",ylab="Age",
	zlab=list("Log death rate",rot=90)) 
dev.off()

postscript("3dplot.female.2.ps")
wireframe((model.final$linear.predictors-L_female_exp)[Age<101]~Year[Age<101]*Age[Age<101],
	scales = list(arrows = FALSE), drape = TRUE, 
	colorkey = FALSE, screen = list(z = 30, x = -60),xlab="Year",ylab="Age",
	zlab=list("Log death rate",rot=90),zlim=c(-11,1),ylim=c(0,109)) 
dev.off()