surveillance (version 1.8-0)

hhh4: Random effects HHH model fit as described in Paul and Held (2011)

Description

Fits a Poisson or negative binomial model with conditional mean $$\mu_{it} = \lambda_{it} y_{i,t-1} + \phi_{it} \sum_{j\neq i} w_{ji} y_{j,t-1} + e_{it} \nu_{it}$$ containing epidemic and endemic components to a multivariate time series of counts. In the case of a negative binomial model, the conditional variance is $\mu_{it}(1+\psi\mu_{it})$ with overdispersion parameter $\psi$. The three unknown quantities of the mean $\mu_{it}$

• $\lambda_{it}$in the autoregressive (ar) component,
• $\phi_{it}$in the neighbor-driven (ne) component, and
• $\nu_{it}$in the endemic (end) component,

are decomposed additively on the log scale, $w_{ji}$ are known weights, and $e_{it}$ is a (multiplicative) offset. The linear predictors may contain random effects. It is also possible to include multiplicative offsets in the autoregressive and neighbor-driven components.

Usage

hhh4(stsObj, 
     control = list(
               ar = list(f = ~ -1,
                         offset = 1,
                         lag = 1,
                         weights = NULL,
                         initial = NULL
                         ),
               ne = list(f = ~ -1,
                         offset = 1,
                         lag = 1,
                         weights = neighbourhood(stsObj),
                         initial = NULL
                         ),
               end = list(f = ~ 1,
                          offset = 1,
                          initial = NULL
                          ),
               family = c("Poisson","NegBin1","NegBinM"),
               subset = 2:nrow(stsObj),
               optimizer = list(stop = list(tol=1e-5, niter=100),
                                regression = list(method="nlminb"),
                                variance = list(method="nlminb")),
               verbose = FALSE,
               start = list(fixed=NULL, random=NULL, sd.corr=NULL),
               data = list(t = epoch(stsObj)-1),
               keep.terms = FALSE
               ),
     check.analyticals = FALSE
   )

Arguments

stsObj

object of class "sts" containing the multivariate count data time series

control

control object, which is a list containing several components: [object Object],[object Object],[object Object],[object Object],[object Object],[object Object],[object Object],[object Object],[object Object],[object Object]

check.analyticals

logical (or a subset of c("numDeriv", "maxLik")), indicating if (how) the implemented analytical score vector and Fisher information matrix should be checked against numerical derivatives at the parameter starting values, usin

Value

• Returns an object of class "hhh4" with elements
• coefficients
  named vector with estimated (regression) parameters of the model
• se
  estimated standard errors (for regression parameters)
• cov
  covariance matrix (for regression parameters)
• Sigma
  estimated variance-covariance matrix of random effects
• Sigma.orig
  estimated variance parameters on internal scale used for optimization
• Sigma.cov
  inverse of marginal Fisher information (on internal scale), i.e., the asymptotic covariance matrix of Sigma.orig
• call
  the matched call
• dim
  vector with number of fixed and random effects in the model
• loglikelihood
  (penalized) loglikelihood evaluated at the MLE
• margll
  (approximate) log marginal likelihood should the model contain random effects
• convergence
  logical. Did optimizer converge?
• fitted.values
  fitted mean values $\mu_{i,t}$
• control
  control object of the fit
• terms
  the terms object used in the fit if keep.terms = TRUE and NULL otherwise
• stsObj
  the supplied stsObj
• lags
  named integer vector of length two containing the lags used for the epidemic components "ar" and "ne", respectively. The corresponding lag is NA if the component was not included in the model.
• nObs
  number of observations used for fitting the model
• nTime
  number of time points used for fitting the model
• nUnit
  number of units (e.g. areas) used for fitting the model
• runtime
  the proc.time-queried time taken to fit the model, i.e., a named numeric vector of length 5 of class "proc_time"

encoding

latin1

Details

For further details see vignette("hhh4") and the references.

References

Held, L., H�hle{Hoehle}, M., Hofmann, M. (2005): A statistical framework for the analysis of multivariate infectious disease surveillance counts. Statistical Modelling, 5, 187--199. Paul, M., Held, L. and Toschke, A. M. (2008): Multivariate modelling of infectious disease surveillance data. Statistics in Medicine, 27, 6250--6267.

Paul, M. and Held, L. (2011): Predictive assessment of a non-linear random effects model for multivariate time series of infectious disease counts. Statistics in Medicine, 30, 1118--1136

Held, L. and Paul, M. (2012): Modeling seasonality in space-time infectious disease surveillance data. Biometrical Journal, 54, 824--843 Meyer, S. and Held, L. (2014): Power-law models for infectious disease spread. Annals of Applied Statistics. In press. arXiv:1308.5115 arXiv-Link: http://arxiv.org/abs/1308.5115

See Also

algo.hhh, fe, ri

Examples

#####################################################################
# Fit some models from ?algo.hhh
#####################################################################

## univariate salmonella agona data
data(salmonella.agona)
# convert to sts class
salmonella <- disProg2sts(salmonella.agona)

# generate formula for temporal and seasonal trends
f.end <- addSeason2formula(f = ~ 1 + t, S=1, period=52)
model1 <- list(ar = list(f = ~ 1), end = list(f =f.end),
               family = "NegBin1")
# run model
res <- hhh4(salmonella, model1)
summary(res, idx2Exp=1, amplitudeShift=TRUE)


## multivariate time series: 
# measles cases in Lower Saxony, Germany
data(measles.weser)
measles <- disProg2sts(measles.weser)

# same model as above
summary(hhh4(measles, control=model1))

# now use region-specific intercepts in endemic component
f.end2 <- addSeason2formula(f = ~ -1 + fe(1, which=rep(TRUE, ncol(measles))) + t,
                            S = 1, period = 52)
model2 <- list(ar = list(f = ~ 1), 
               end = list(f = f.end2, offset = population(measles)),
               family = "NegBin1")
# run model
summary(hhh4(measles, control=model2), idx2Exp=1, amplitudeShift=TRUE)

# include autoregressive parameter phi for adjacent "Kreise"
# no linear trend in endemic component
f.end3 <- addSeason2formula(f = ~ -1 + fe(1, which=rep(TRUE, ncol(measles))), 
                            S = 1, period = 52)
model3 <- list(ar = list(f = ~ 1),
               ne = list(f = ~1),
               end = list(f = f.end3, offset= population(measles)),
               family = "NegBin1")
# run model
res3 <- hhh4(measles, control=model3)
summary(res3, idx2Exp=1:2, amplitudeShift=TRUE)

## check that neighbourhood weights array yields same result
.neweights <- array(neighbourhood(measles),
                   dim = c(rep(ncol(measles),2),nrow(measles)),
                   dimnames = c(dimnames(neighbourhood(measles)), list(NULL)))
res3_tv <- hhh4(measles, control = modifyList(model3,
                list(ne=list(weights=.neweights))))
.disregardComponents <- function(obj) {
    obj$call <- obj$control$ne$weights <- obj$runtime <- NULL
    obj
}
stopifnot(all.equal(.disregardComponents(res3), .disregardComponents(res3_tv)))

######################################################################
# Fit the models from the Paul & Held (2011) paper for the influenza data
# from Bavaria and Baden-Wuerttemberg (this takes some time!)
# For further documentation see also the vignette.
######################################################################

data("fluBYBW") 

###############################################################
## generate formula for temporal and seasonal trends
f.end <- addSeason2formula(f = ~ -1 + ri(type="iid", corr="all") + 
                               I((t-208)/100), S=3, period=52)

## details for optimizer
opt <- list(stop = list(tol=1e-5, niter=200),
            regression = list(method="nlminb"),
            variance = list(method="nlminb"))

##########################
## models 
# A0
cntrl_A0 <- list(ar = list(f = ~ -1),
                 end = list(f =f.end, offset = population(fluBYBW)),
                 family = "NegBin1", optimizer = opt, verbose = 1)
summary(res_A0 <- hhh4(fluBYBW,cntrl_A0))

# B0
cntrl_B0 <- list(ar = list(f = ~ 1),
                 end = list(f =f.end, offset = population(fluBYBW)),
                 family = "NegBin1", optimizer = opt, verbose=1)
res_B0 <- hhh4(fluBYBW,cntrl_B0)               
 

# C0
cntrl_C0 <- list(ar = list(f = ~ -1 + ri(type="iid", corr="all")),
                 end = list(f =f.end, offset = population(fluBYBW)),
                 family = "NegBin1", optimizer = opt, verbose=1)
res_C0 <- hhh4(fluBYBW,cntrl_C0)               


#A1

# weight matrix w_ji = 1/(No. neighbors of j) if j ~ i, and 0 otherwise
wji <- neighbourhood(fluBYBW)/rowSums(neighbourhood(fluBYBW))

cntrl_A1 <- list(ar = list(f = ~ -1),
                 ne = list(f = ~ 1, weights = wji),
                 end = list(f =f.end, offset = population(fluBYBW)),
                 family = "NegBin1", optimizer = opt, verbose=1)
res_A1 <- hhh4(fluBYBW,cntrl_A1)               


# B1
cntrl_B1 <- list(ar = list(f = ~ 1),
                 ne = list(f = ~ 1, weights = wji),
                 end = list(f =f.end, offset = population(fluBYBW)),
                 family = "NegBin1", optimizer = opt, verbose=1)
res_B1 <- hhh4(fluBYBW,cntrl_B1)               


# C1
cntrl_C1 <- list(ar = list(f = ~ -1 + ri(type="iid", corr="all")),
                 ne = list(f = ~ 1, weights = wji),
                 end = list(f =f.end, offset = population(fluBYBW)),
                 family = "NegBin1", optimizer = opt, verbose=1)
res_C1 <- hhh4(fluBYBW,cntrl_C1)               


#A2
cntrl_A2 <- list(ar = list(f = ~ -1),
                 ne = list(f = ~ -1 + ri(type="iid",corr="all"), weights=wji),
                 end = list(f =f.end, offset = population(fluBYBW)),
                 family = "NegBin1", optimizer = opt, verbose=1)
res_A2 <- hhh4(fluBYBW,cntrl_A2)               


# B2
cntrl_B2 <- list(ar = list(f = ~ 1),
                 ne = list(f = ~ -1 + ri(type="iid",corr="all"), weights =wji),
                 end = list(f =f.end, offset = population(fluBYBW)),
                 family = "NegBin1", optimizer = opt, verbose=1)
res_B2 <- hhh4(fluBYBW,cntrl_B2)               

# C2
cntrl_C2 <- list(ar = list(f = ~ -1 + ri(type="iid", corr="all")),
                 ne = list(f = ~ -1 + ri(type="iid",corr="all"), weights =wji),
                 end = list(f =f.end, offset = population(fluBYBW)),
                 family = "NegBin1", optimizer = opt, verbose=1,
                 start=list(fixed=fixef(res_B0),random=c(rep(0,140),
                         ranef(res_B0)), sd.corr=c(-.5,res_B0$Sigma.orig,0)))
res_C2 <- hhh4(fluBYBW,cntrl_C2)               


# D
cntrl_D <- list(ar = list(f = ~ 1),
                ne = list(f = ~ -1 + ri(type="iid"), weights = wji),
                end = list(f =addSeason2formula(f = ~ -1 + ri(type="car") + 
                                             I((t-208)/100), S=3, period=52), 
                          offset = population(fluBYBW)),
                family = "NegBin1", optimizer = opt, verbose=1)
res_D <- hhh4(fluBYBW,cntrl_D)