#Crump 2014 Season of birth and other perinatal risk factors for melanoma
#Int. J. Epidemiol. (2014) doi: 10.1093/ije/dyt277
#Table 1 Data
dat<-structure(list(`Birth Month` = c("January", "February", "March",
"April", "May", "June", "July", "August", "September", "October",
"November", "December"), Births = c(297241L, 291749L, 334787L,
331520L, 324466L, 305267L, 310043L, 300578L, 291196L, 278533L,
252598L, 253596L), `Person-years (Millions)` = c(5.47, 5.36,
6.2, 6.09, 5.88, 5.45, 5.44, 5.26, 5.13, 4.84, 4.41, 4.42), `CMM cases` = c(144L,
139L, 162L, 170L, 166L, 127L, 131L, 126L, 116L, 115L, 99L, 100L
), `Cases per 100 000 person-years` = c(2.63, 2.59, 2.61, 2.79,
2.82, 2.33, 2.41, 2.4, 2.26, 2.37, 2.25, 2.26)), .Names = c("Birth Month",
"Births", "Person-years (Millions)", "CMM cases", "Cases per 100 000 person-years"
), row.names = c(NA, 12L), class = "data.frame")

#Table 2 Data
age.brackets<-seq(0,30,by=5)+2
age.cases<-c(4,6,32,181,395,505,472)



#Simulate researcher behaviour
bday<-sample(1:360, 100000, replace=T)
tday<-sample(120:210, 100000, replace=T)

age.sim<-cbind(bday,tday,0,NA)
colnames(age.sim)[3]<-"Age at Test"
age.sim[which(age.sim[,1]<=age.sim[,2]),3]<-1
age.sim<-age.sim[sort(age.sim[,1], index.return=T)$ix,]

mo<-1
for(i in seq(0,330,by=30)){
age.sim[which(age.sim[,1]>i & age.sim[,1]<=(i+30)),4]<-mo
mo=mo+1
}

pred.ages=NULL
for(i in 1:12){
pred.ages<-rbind(pred.ages,mean(age.sim[which(age.sim[,4]==i),3]))
}




#estimate cases by month
rate.per.yr<-(395-181)/5
rate.per.month<-rate.per.yr/12
agesize<-3595055/7
dec.cases<-rate.per.month*6 + 181
births.per.mo<-dat[1:12,2]
mo<-12*pred.ages
est<-((dec.cases+mo*rate.per.month)/agesize)*births.per.mo


#plots
plot(1000000*dat[,3]/dat[,2], type="b", xlab="", xaxt="n",
ylab= "person-years of follow up per birth",
main= "Average age?")
axis(side=1,at=1:12,labels=dat[,1])
lines(1:12,17.4+pred.ages, type="l", col="Red", lwd=3)

par(mfrow=c(2,1))
plot(age.brackets,age.cases, type="b",
xlab="Age Bracket (Middle Value)", ylab="# of Cases")
abline(v=c(17.5,18.5), lty=2, lwd=3)
abline(h=seq(0,500,by=100),lty=1)
abline(h=seq(0,500,by=50),lty=3)

plot(age.brackets,age.cases/agesize, type="b",
xlab="Age Bracket (Middle Value)", ylab="Estimated Prevalence")
abline(v=c(17.5,18.5), lty=2, lwd=3)
abline(h=seq(0,1,by=.2)/1000,lty=1)
abline(h=seq(0,1,by=.1)/1000,lty=3)


par(mfrow=c(2,1))
plot(dat[,4], type="b", xaxt="n",
ylim=c(min(dat[,4],est),max(dat[,4],est)),
ylab=colnames(dat)[4], xlab="Birth Month",
main="# of CMM cases")
axis(side=1,at=1:12,labels=dat[,1])
lines(est, type="b", col="Blue")


plot(est-dat[,4], type="b", col="Blue",
ylab="Model-data",  xaxt="n", xlab="Birth Month",
main="# of CMM cases")
axis(side=1,at=1:12,labels=dat[,1])
abline(h=0)