## com39a of october 2021,
#  old roman era egypt data again, now for bayes stk 4021,
#  with Weibull models, lazy bayes approximations, etc. 

#  two questions addressed here, via draws from posterior distributions: 
#  (a) gamma = Prob(live longer than 40.0 years), men, women, and ratio 
#  (b) mu = median = a (log 2)^(1/b), men, women, and difference 

library(MASS)
# to get mvrnorm on board 

egypt <- scan("egypt_data", skip=5)
# range: 1.5 to 96.0 

yym <- egypt[1:82] # 82 men 
nnm <- length(yym)
yyw <- egypt[83:141] # 59 women 
nnw <- length(yyw)

## looking at empirical survival curves: 
yval <- seq(0,97,by=0.10)
Shatm <- 0*yval 
Shatw <- 0*yval 

for (j in 1:length(yval))
{
Shatm[j] <- mean(1*(yym >= yval[j]))
Shatw[j] <- mean(1*(yyw >= yval[j]))
}

matplot(yval,cbind(Shatm,Shatw),type="l",
   lwd=2.5,xlab="age",ylab="survival curves")

# F = 1 - exp(-(y/a)^b) 
# f = exp(-(y/a)^b) b y^(b-1) / a^b

## loglik analysis for the men: 
logLm <- function(ab)
{
a <- ab[1]
b <- ab[2]
# 
auxm <- -(yym/a)^b + log(b) + (b-1)*log(yym) - b*log(a)
sum(auxm)
}

minuslogLm <- function(ab)
{-logLm(ab)}

nilsm <- nlm(minuslogLm,c(10,1),hessian=T)
MLm <- nilsm$estimate
Jhatm <- nilsm$hessian 
Sigmam <- solve(Jhatm) 
sem <- sqrt(diag(Sigmam))
showm <- cbind(MLm,sem)
print(round(showm,4))

#         MLm    sem
#[1,] 36.9910 3.1393 a 
#[2,]  1.3617 0.1261 b 

## loglik analysis for the women: 
logLw <- function(ab)
{
a <- ab[1]
b <- ab[2]
# 
auxw <- -(yyw/a)^b + log(b) + (b-1)*log(yyw) - b*log(a)
sum(auxw)
}

minuslogLw <- function(ab)
{-logLw(ab)}

nilsw <- nlm(minuslogLw,c(10,1),hessian=T)
MLw <- nilsw$estimate
Jhatw <- nilsw$hessian 
Sigmaw <- solve(Jhatw) 
sew <- sqrt(diag(Sigmaw))
showw <- cbind(MLw,sew)
print(round(showw,4))
 
#         MLw    sew
#[1,] 28.8095 2.5549 a 
#[2,]  1.5457 0.1530 b 

## we're now prepared for (lazy) Bayesian analysis
#  of the relevant interest parameters 

#  prob = Prob(live longer than 40.0 year) 
y0 <- 40.0 
sim <- 10^4

thetam <- mvrnorm(sim,MLm,Sigmam)
plot(thetam) 
am <- thetam[ ,1]
bm <- thetam[ ,2]
probm <- exp(-(y0/am)^bm)

thetaw <- mvrnorm(sim,MLw,Sigmaw)
plot(thetaw) 
aw <- thetaw[ ,1]
bw <- thetaw[ ,2]
probw <- exp(-(y0/aw)^bw)

print( round(quantile(probm,c(0.05,0.50,0.95)),4) ) 
print( round(quantile(probw,c(0.05,0.50,0.95)),4) ) 

#    5%    50%    95% 
# 0.2594 0.3291 0.3955 men   : much more likely to have long lives 
# 0.1239 0.1913 0.2579 women : bad life prospects 

## looking at the medians:   
medm <- am*(log(2))^(1/bm)
medw <- aw*(log(2))^(1/bw)

plot(medm,medw)

delta <- medm - medw 

print( round(quantile(medm,c(0.05,0.50,0.95)),4) ) 
print( round(quantile(medw,c(0.05,0.50,0.95)),4) ) 
print( round(quantile(delta,c(0.05,0.50,0.95)),4) ) 

plot(density(delta))

#     5%     50%     95% 
# 23.8205 28.2126 32.6861 men
# 19.0247 22.6906 26.3848 women 
# -0.1954  5.5090 11.3031 median men - median women 

## looking at empirical survival curves,
## now together with the ML fitted curves: 

yval <- seq(0,97,by=0.10)
Shatm <- 0*yval 
Shatw <- 0*yval 

Shatmpara <- exp(-(yval/MLm[1])^MLm[2])
Shatwpara <- exp(-(yval/MLw[1])^MLw[2])

for (j in 1:length(yval))
{
Shatm[j] <- mean(1*(yym >= yval[j]))
Shatw[j] <- mean(1*(yyw >= yval[j]))
}

par(bty="l")
par(cex.lab=1.4142) 

matplot(yval,cbind(Shatm,Shatw, Shatmpara,Shatwpara),
   type="l",col=c(1,1,2,2),lty=c(1,2,1,2),
   lwd=2.5,xlab="age",ylab="survival curves, nonpara and Weibull")

text(60,0.40,"men")
text(20,0.25,"women")