serf.fit: Fit Serfing's regression model
In ensoesie/epidemics: A Method for Defining Influenza Epidemic Onset

Description Usage Arguments Value Examples

Takes a time series and fits a regression model using the method described by Serfling (1963) in 'Methods for Current Statistical Analysis of Excess Pneumonia-Influenza Deaths'. Helper function 'serf.design' constructs the design matrix according to parameters 'trend.poly', 'harmonic', and 'period'.

serf.fit(ts = g.ts, t1 = 1, t2 = length(ts), t3 = min(t2 + 1,
  length(ts)), t4 = min(t2 + 1, length(ts)), trend.poly = 1,
  harmonic = c(1, 2), period = 52.2, ylim = range(ts[t1:t4]),
  plot = TRUE, output = TRUE, thresh = 2.5)

`ts`	response vector of numeric values representing the time series to be modeled. By default this is set to 'g.ts' the global time series which MUST be specified by the user before calling this function
`t1`	is the integer index of the parameter 'ts' where model fitting begins; defaults to 1
`t2`	is the integer index of the parameter 'ts' where model fitting ends; defaults to length(ts)
`t3`	is the integer index of the time series where fitted values begin; defaults to min(t2+1,length(ts))
`t4`	is the integer index of the time series where fitted values end; defaults to min(t2+1,length(ts))
`trend.poly`	is specified as an integer vector representing the degrees of the polynomial terms to be included; defaults to 1 i.e. single linear trend term
`harmonic`	is also specified as an integer vector representing fractional multiples of the period e.g. 'harmonic=c(1,2)' specifies two harmonic pairs (two sine terms and two cosine terms) of one period and one half-period respectively
`period`	numeric value; defaults to 52.2
`ylim`	range of y-axis for plotting
`plot`	logical, generates plot of times series with threshold and epidemic observations; defaults to TRUE
`output`	logical, outputs fitted values and epidemic observations; default set to TRUE
`thresh`	is a numeric value specifying the epidemic threshold as a multiple of the residual standard deviation; defaults to 2.5

epidemic numeric 0/1, 1 indicates epidemic period and 0 otherwise

library(splines)
          library(MASS)
          library(date)

          options(warn=2)

          data(wa_who) #load WHO ILI data for Washington state
          names(wa_who) <- c('id','ili','total','iliperc') #rename columns

          head(wa_who)

          wa_who$year <- floor(wa_who/100)[1] #extract year
          wa_who$week <- (wa_who %% 100)[1]   #extract week

          names(wa_who) <- c('id','ili','total','iliperc', 'year', 'week') #add column names
          datelist <- sort(unique(wa_who$id)) #save sorted dates
          year <- floor(datelist/100)
          week <- datelist %% 100

          ili <- rep(NA,length(datelist))
          tot <- rep(NA,length(datelist))

          for (i in 1:length(datelist)) {
           ili[i] <- sum(wa_who$ili[wa_who$id==datelist[i]],na.rm=TRUE)
           tot[i] <- sum(wa_who$total[wa_who$id==datelist[i]],na.rm=TRUE)
           }

           ili2.dat <- data.frame(id=(1:length(datelist)), year, week, ili=ili,total=tot)
           ili2.dat$perc <- 100*round(ili2.dat$ili/ili2.dat$total,3)


           # Global data and plotting parameters
           # Change as needed

           g.ts <- ili2.dat$perc
           g.ts[is.nan(g.ts)] <- 0
           g.period <- 52.2
           g.mon <- 'Oct'
           g.yr <- 1997 #data start year
           g.xlab <- 'Calendar year'
           g.ylab <- 'Percentage of samples positive'
           g.main <- 'WHO'
           g.dates <- TRUE

           resvec <- rep(0,length(g.ts))
           stoptime <- length(g.ts) - 259
           resmat <- matrix(NA,stoptime,length(g.ts))

           for (i in 1:stoptime) resmat[i,i:(i+259)] <- serf.fit(g.ts, t1=i,t2=i+259,harmonic=1,thresh=1.96,output=T)$epi.thresh

           for (i in 1:length(g.ts)) resvec[i] <- max(resmat[,i],na.rm=TRUE)

           ilires.dat <- data.frame(year=ili2.dat$year,week=ili2.dat$week,perc=ili2.dat$perc,epidemic=resvec)
           #write.csv(ilires.dat, ' ', row.names=FALSE) #write output to file


           #pdf('who-timeseries.pdf') #save plot to desired location
           #tsplot()
           #points(which(resvec==1),g.ts[resvec==1],pch=19,col='red',cex=0.8)
           #dev.off()