Documents/R/nowcast.R

#####################################################################
# Test script to compare delay adjustment with the one made in
# Hoehle and an der Heiden (2014), Biometrics, See:
# http://people.su.se/~mhh/pubs/hoehle_anderheiden2014-preprint.pdf
# This functionality is implemented in the R package surveillance
# as function 'nowcast'.
####################################################################

#Load our package.
library("euromomo")

#Load surveillance package.
library("surveillance")

#Small helper function to take the global momoFile$dLastFullWeek
#and take nWeeks back.
nWeeksFromNow <- function(nWeeks) {
    sort(seq(momoFile$dLastFullWeek,by="-1 week",length.out=nWeeks))
}
#Small helper function to concatenate the last element of a vector on the vector
repLast <- function(x) { c(x,x[length(x)])}

#Read options using the example data.
parseDefaultsFile("defaults.txt")
checkOptions()

# Read in the raw data
momoFile <- readmomofile()

#===========================================================
# Modify Day of Aggregation for reading all data
optsAll <- optsBackup <- getOption("euromomo")
optsAll[["DayOfAggregation"]] <- Sys.Date()
optsAll[["BaselineSeasons"]] <- 7
options("euromomo"=optsAll)
#Read data with this date
momoFileAll <- readmomofile()
#Restore old options.
options("euromomo"=optsBackup)
#====================================================

### read holidays
holiday.file<-holiday(holiday.filename=getOption("euromomo")$HolidayFile)

#Create the groups (as stored in the option file)
momo <- makeGroups(momoFile$momo)
momoAll <- makeGroups(momoFileAll$momo)

### The groups are determined from the options
groups<-names(getOption("euromomo")[["groups"]])

### Make a PDF with all 4+1 default age groups
pdf(file="delayadjustments.pdf",width=8,height=5)
par(mfrow=c(2,3))

#Do comparsion for one group (Here it's the 0-4)
for (i in groups) {
  #i<-"momodefault1"
  groupOpts <- getOption("euromomo")[["groups"]][[i]]
  groupLabel <- groupOpts["label"]
  cat("Group",groupOpts["label"],"\n")

  back<-as.numeric(groupOpts["back"])
  groupIndicator <- momo[, paste0("group_",i)]
  momoGroup <- momo[groupIndicator,]
  momoGroupAll <- momoAll[momoAll[,paste0("group_",i)],]

  #Get reporting triangle
  rTList <- df2ReportingTriangle(momo, groupIndicator, back, dWeeks=momoFile$dWeeks, dLastFullWeek=momoFile$dLastFullWeek) # something about the group
  rTDF <- rT2DataFrame(rTList$cumRT)

  #Make sts object with the same information.
  sts <- linelist2sts( momoGroup, dateCol="DoD", aggregate.by="1 week",dRange=momoFile$dWeeks)
  stsAll <- linelist2sts( momoGroupAll, dateCol="DoD", aggregate.by="1 week",dRange=momoFileAll$dWeeks)

  #Can plot it, if we want to.
  plotPeriod <- nWeeksFromNow(10)
  whichPlot <- epoch(sts) %in% plotPeriod

  if (FALSE) {
    #Show data as weekly data. To understand the plotting symbols check strptime and
    plot(sts[whichPlot,],xaxis.tickFreq=list("%m"=atChange,"%Y"=atChange),
         xaxis.labelFreq=list("%V"=atChange),
         xaxis.labelFormat="%G-W%V",legend.opts=NULL,
         xlab="",las=1,cex.axis=0.8,main=groupLabel,
         ylim=c(0,max(observed(sts[whichPlot,]),observed(stsAll[whichPlot,]))),
         lwd=1)
    lines(seq_len(sum(whichPlot)+1)-0.5,repLast(observed(stsAll[whichPlot,])),type="s",lwd=3,col="darkgreen")
    legend(x="topright",legend=paste0("DoA=",
                                      c(options("euromomo")$euromomo$DayOfAggregation,as.character(optsAll[["DayOfAggregation"]]))),lwd=c(1,3),col=c("black","darkgreen"))
  }

  #Prior with a large variance
  prior <- structure("poisgamma",mean.lambda=mean(observed(sts)),var.lambda=5*var(observed(sts)))
  #noWeeks between dLastFullWeek and StartDelayEst
  nWeeksDelayEst <- as.numeric((momoFile$dLastFullWeek - ISOweek2date(paste(getOption("euromomo")[["StartDelayEst"]],"-1",sep=""))) / 7)

  #Call nowcast function from R package surveillance.
  nc <- nowcast(now=momoFile$dLastFullWeek,when=nWeeksFromNow(back+1),
                data=momoGroup,dEventCol="DoD",dReportCol="DoR",
                method=c("bayes.trunc"), #,"bayes.notrunc.bnb","lawless","bayes.trunc"),
                aggregate.by="1 week",
                D=back,
                m=nWeeksDelayEst, #Moving window for the delay estimation. NULL=take all ATM it takes ALL.
                control=list(dRange=plotPeriod,N.tInf.max=max(observed(sts))*3,N.tInf.prior=prior))

  # Delay adjustment using the momo algorithm.
  drTDF <- delay.nb(rTDF,holiday=holiday.file)

  #Show time series and posterior median forecast/nowcast
  plot(nc,xaxis.tickFreq=list("%d"=atChange,"%m"=atChange),
       xaxis.labelFreq=list("%d"=at2ndChange),
       xaxis.labelFormat="%G-W%V",
       xlab="Time (days)",ylab="No.",
       col=c(NA,"gray","blue"),lty=c(1,1,1,1),lwd=c(1,1,3),legend.opts=NULL,
       main=groupLabel, ylim=c(0,max(observed(stsAll[whichPlot,]),nc@pi,na.rm=TRUE)))

  #Show truth also
  idx <- which( epoch(stsAll) %in% epoch(nc))
  lines(seq_len(length(idx)+1)-0.5,repLast(observed(stsAll[idx,])),col="magenta",lwd=3,type="s")

  # Which index to show of drTDF
  idxShow <- which(drTDF$ISOweek %in% ISOweek(plotPeriod))
  lines(seq_len(length(idxShow)), drTDF[idxShow,"cnb"],lwd=3,type="b")

  lines(seq_len(length(idxShow)), drTDF[idxShow,"cnb"] - 1.96*sqrt(drTDF[idxShow,"v.cnb"]),lwd=1,lty=2)
  lines(seq_len(length(idxShow)), drTDF[idxShow,"cnb"] + 1.96*sqrt(drTDF[idxShow,"v.cnb"]),lwd=1,lty=2)
}

plot(c(0,0),type="n",axes=FALSE,xlab="",ylab="")
#legend(x="center",c("bayes.trunc","delay.nb"),col=c("blue","black"),lwd=3,bg="white",lty=c(1,1))
legend(x="center",c("bayes.trunc","delay.nb","truth"),col=c("blue","black","magenta"),lwd=3,bg="white",lty=c(1,1,1))

#Close graphics device.
dev.off()
thl-mjv/euromomo documentation built on May 31, 2019, 10:43 a.m.