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'.
1 2 3 4 |
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
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 | 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()
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.