#  analyses of titanic data     jh 2009.01.28

library(Epi)

######  titanic data ###############
 
titanic = read.table("titanic.txt",header=T)
titanic[1:10,]
titanic[427:436,]

# if wish, drop passenger names;

titanic=titanic[,3:7]
summary(titanic)

# Define as Lexis object with timescales calendar time and age 

L=Lexis(entry=list(per=1912),exit=list( per=last, age=last-born),exit.status = dead,data=titanic)
names(L)
plot(L,cex.lab=2,xlab="year",cex.axis=1.5)
summary(L$lex.dur)

LL <- splitLexis(  L, breaks = seq(1911,2005,5), time.scale="per")
LL <- splitLexis( LL, breaks = c(0,1,seq(0,105,5)), time.scale="age" )
length(LL$lex.dur)
LL[1:30,]

LL.male=LL[LL$ifemale==0,];   LL.male[101:124,]
LL.female=LL[LL$ifemale==1,]; LL.female[101:124,]
                        
                        
########## USA death.rates #########################################

us.rates=read.table("USmx-5x5.txt",header=T) ; dim(us.rates)

# augment .. use 1991-1995 for 1996-2000 and for 2001-2005

us.rates=rbind(us.rates,us.rates[451:475,],us.rates[451:475,])
dim(us.rates)

# drop 1901-05 and 1906-10

us.rates=us.rates[51:525,] ; dim(us.rates)

# array them (dropping ages 105-)

us.males  =array( us.rates[,4], c(25,19) ); us.males  = us.males[1:22,];
us.females=array( us.rates[,3], c(25,19) ); us.females= us.females[1:22,];
rm(us.rates)
                          
## log death rates by age band

par(mfrow = c(1,2), # 2 x 2 pictures on one plot
          pty = "s")       # square plotting region,
                           # independent of device size
# 1911-1915

plot(1:22,log2(us.males[1:22,1]),
    col="blue",main="log2 rates, 1911-1915",xlab="age band",
    ylim=c(-13,-1),cex.lab=2,cex.main=2,cex.axis=2,cex=2,ylab="")
points(1:22,log2(us.females[1:22,1]),
    col="red",ylab="",cex.lab=2,cex.main=2,cex.axis=2,cex=2)

# 1991-1995

plot(1:22,log2(us.males[1:22,19]),col="blue",main="log2 rates, 1991-1995",xlab="age band",cex.lab=2,cex.main=2,cex.axis=2,cex=2,ylab="",ylim=c(-13,-1))
points(1:22,log2(us.females[1:22,19]),col="red",cex.lab=2,cex.main=2,cex.axis=2,cex=2)

############# death rates in 3 age-bands, from 1910 to 1995

par(mfrow = c(1,3),pty = "s")  

plot(1:17,log2(us.males[1,1:17]),col="blue",ylim=c(-11,-3),main="log2 rates, 0-1 age band",xlab="17 quinquennia 1911-1995",ylab="",cex.lab=2,cex.main=2,cex.axis=2,cex=2)
points(1:17,log2(us.females[1,1:17]),col="red",cex=2)

plot(1:17,log2(us.females[7,1:17]),col = "red",cex.lab=2,cex.main=2,cex.axis=2,cex=2,
ylab="",ylim=c(-11,-3),main="log2 rates, 25-29 age band",xlab="17 quinquennia 1911-1995")
points(1:17,log2(us.males[7,1:17]),col = "blue",cex=2)

plot(1:17,log2(us.females[15,1:17]),col = "red",cex.lab=2,cex.main=2,cex.axis=2,ylab="",cex=2,
ylim=c(-11,-3),main="log2 rates,  65-69 age band",xlab="17 quinquennia 1911-1995")
points(1:17,log2(us.males[15,1:17]),col = "blue",cex=2)
                           

######### Death rates in Titanic survivors vs US General Population ##########

# Overall no. of deaths and p-t

observed.deaths = sum(LL[,"lex.Xst"]); observed.deaths
person.time = sum(LL[,"lex.dur"]); person.time


# ************** aggregates within age-period cells, using aggregate ***********
	
# age categories (a.c for short) 

# deal with 01- and 1-4 .. ie convert to categories

LL$a.c = 1*(LL$age==0) + 2 * ( (LL$age > 0) & (LL$age < 6) ) + 
         (2+floor(LL$age/5)) * ( LL$age > 5) ; LL[1:30,c("age","a.c")]

# lex.Xst and lex.dur are the numbers of deaths and numebrs of years
# in the lexis rectangle

#  t  is short for total  

# p is short for (calendar) period ; m.cat for whether male

t.us=aggregate(LL[,c("lex.Xst","lex.dur")],      
           by=list(a.cat=LL$a.c,
                   p.cat=floor((LL$per-1907)/5),
                   m.cat=1-LL$ifemale,
                   c.cat=LL$class ),
           sum) ; 
        
       dim(t.us); t.us[1:10,]
 
  #  augment with a column of rates and one of null expectations
   
  t.us$usrate = NA ;
  n.rows = length(t.us$a.cat)
  
  for (k in 1:n.rows) {
  	i = t.us[k,"a.cat"]; j = t.us[k,"p.cat"]; s=t.us[k,"m.cat"];
  	if (s==1) t.us$usrate[k] = us.males[i,j] else t.us$usrate[k] = us.females[i,j]
  	} 
  	
  t.us[1:20,]  
  
  # titanic (index) vs USA   hr.mh = O/E (see footnote in assignment)
  
  mh.num = sum( t.us$lex.Xst ) # observed deaths
  mh.den = sum( t.us$usrate * t.us$lex.dur ) # expected deaths
  mh.ratio = mh.num/mh.den
  
  c( round(mh.num,2), round(mh.den,2), round(mh.ratio,2) )
  
  # males .. subset then go back to hratio  t.us=t.us[t.us$m.cat==0,]  
              
  null.chi = abs ( (mh.num - mh.den ) / sqrt(mh.den) )
   
  c("ChiSq=",round(null.chi^2,2), "Chi",round(null.chi,2), "P-value", 
    round(2*(1-pnorm(null.chi)),4) ) 
  
  test.based.ci.95 = round(mh.ratio^(1+1.96*c(1,-1)/null.chi),2) ; test.based.ci.95   
     
  # --