surveillance (version 1.5-4)

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.

Usage

hhh4(stsObj, 
     control = list(
               ar = list(f = ~ -1,
                         lag = 1,
                         weights = NULL,
                         initial = NULL
                         ),
               ne = list(f = ~ -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 = data.frame(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 ah4 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 components for random effects
• Sigma.orig
  estimated variance components on internal scale used for optimization
• Sigma.cov
  inverse of (minus Hessian matrix) for variance components
• 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
• lag
  specified lag of the autoregression, ATM always $= 1$
• 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

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

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
summary(hhh4(measles, control=model3), idx2Exp=1:2, amplitudeShift=TRUE)

######################################################################
# 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, data = data.frame(t = epoch(fluBYBW)-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, data=data.frame(t=epoch(fluBYBW)-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, data=data.frame(t=epoch(fluBYBW)-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, data=data.frame(t=epoch(fluBYBW)-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, data=data.frame(t=epoch(fluBYBW)-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, data=data.frame(t=epoch(fluBYBW)-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, data=data.frame(t=epoch(fluBYBW)-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, data=data.frame(t=epoch(fluBYBW)-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, data=data.frame(t=epoch(fluBYBW)-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, data=data.frame(t=epoch(fluBYBW)-1))
res_D <- hhh4(fluBYBW,cntrl_D)