The models we fitted so far, accounted for autoregressive and spatial spread through neighbourhood weights. Meyer and Held (2017), incorporated social contact matrices, to add additional information about contact patterns of different age groups in hhh4contacts. This is what we are going to implement and compare in this chapter.
6.1 Mathematical Formulation
By stratifying the counts also by age groups, a third dimension is incorporated. The local transmission autoregressive part of the model, is now being included in the spatiotemporal part by giving it unit weight in the power-law formulation (\(w_{ji} = (o_{ji} + 1)^{-\rho}\)). So the 2D model formula becomes 3D with endemic and epidemic spatiotemporal-only components:
So the product of the contact matrix \(C\) and the spatial weights \(W\) (normalized), determines how previous counts affect the mean \(\mu_{git}\). By parameterising the contact matrix with a power parameter \(\kappa\), the contact matrix can be adjusted. Then \(\kappa\) measures the amount of transmission between subgroups, where
\[ C^{\kappa} := E \Lambda^{\kappa} E^{-1}, \]
and \(E\) is the matrix of eigenvectors and \(\Lambda\) is the diagonal matrix of eigenvalues. More information on the parameterisation and inference can be found in the correspondind publication from Meyer and Held (2017).
6.2 Using surveillance
The addition of age groups makes the model structure slightly more complex. We follow the same approach as in the demo of hhh4contacts, since we also use the example data provided by this package. For further details, we refer to demo("hhh4contacts").
First, we create an sts object stratified by age group \(g = 1,\dots,M\) and district \(i = 1,\dots,N\). Using hhh4contacts::noroBE(), we request all stratifications and a flattened representation. Thereby, the original dimension \(T \times M \times N\) become reshaped into a matrix of dimension \(T \times (MN)\), where rows correspond to time points and columns to age group-district combinations. Following hhh4contacts, districts vary faster than age groups.
## age and district stratified GROUPING <-c(1, 2, 2, 4, 4, 2) # 6 age groupsnoroBEall <-noroBE(by ="all", flatten =TRUE, agegroups = GROUPING)head(observed(noroBEall), 1)
Using hhh4contacts::stratum() we can extract names of the specific strata.
## names of groups and districtsDISTRICTS <-unique(stratum(noroBEall, 1))NDISTRICTS <-length(DISTRICTS)GROUPS <-unique(stratum(noroBEall, 2))NGROUPS <-length(GROUPS)
For the contact matrix, we use the function contactmatrix() from hhh4contacts. It retrieves various social contact matrices for Germany from the POLYMOD survey (Mossong et al. (2008)).
## aggregate to the same six age groups as the above countsCgrouped <-contactmatrix(which ="reciprocal", # estimated by the Wallinga et al (2006) methodtype ="all",grouping = GROUPING # age group specification)Cgrouped
We then construct indicator matrices for age groups and districts, which are used as fixed effects in the model.
DATAt <-list(t =epoch(noroBEall) -1)## setup a model matrix with group indicatorsMMG <-sapply(GROUPS, function (g) { index <-which(stratum(noroBEall, which =2) == g) res <-col(noroBEall) res[] <- res %in% index res}, simplify =FALSE, USE.NAMES =TRUE)str(MMG)
## setup model matrix with district indicatorsMMR <-sapply(DISTRICTS, function (r) { index <-which(stratum(noroBEall, which =1) == r) res <-col(noroBEall) res[] <- res %in% index res}, simplify =FALSE, USE.NAMES =TRUE)str(MMR)
To allow seasonal patterns to differ between age groups, we create an age-group-specific seasonal model matrix by multiplying each age-group indicator with sine and cosine terms of period 52 weeks.
## setup model matrix of group-specific seasonal termsMMgS <-with(c(MMG, DATAt), unlist(lapply(X = GROUPS,FUN =function (g) { gIndicator <-get(g) res <-list(gIndicator *sin(2* pi * t/52), gIndicator *cos(2* pi * t/52))names(res) <-paste0(c("sin", "cos"), "(2 * pi * t/52).", g) res }), recursive =FALSE, use.names =TRUE))str(MMgS)
List of 12
$ sin(2 * pi * t/52).00-04: num [1:208, 1:72] 0 0.121 0.239 0.355 0.465 ...
$ cos(2 * pi * t/52).00-04: num [1:208, 1:72] 1 0.993 0.971 0.935 0.885 ...
$ sin(2 * pi * t/52).05-14: num [1:208, 1:72] 0 0 0 0 0 0 0 0 0 0 ...
$ cos(2 * pi * t/52).05-14: num [1:208, 1:72] 0 0 0 0 0 0 0 0 0 0 ...
$ sin(2 * pi * t/52).15-24: num [1:208, 1:72] 0 0 0 0 0 0 0 0 0 0 ...
$ cos(2 * pi * t/52).15-24: num [1:208, 1:72] 0 0 0 0 0 0 0 0 0 0 ...
$ sin(2 * pi * t/52).25-44: num [1:208, 1:72] 0 0 0 0 0 0 0 0 0 0 ...
$ cos(2 * pi * t/52).25-44: num [1:208, 1:72] 0 0 0 0 0 0 0 0 0 0 ...
$ sin(2 * pi * t/52).45-64: num [1:208, 1:72] 0 0 0 0 0 0 0 0 0 0 ...
$ cos(2 * pi * t/52).45-64: num [1:208, 1:72] 0 0 0 0 0 0 0 0 0 0 ...
$ sin(2 * pi * t/52).65+ : num [1:208, 1:72] 0 0 0 0 0 0 0 0 0 0 ...
$ cos(2 * pi * t/52).65+ : num [1:208, 1:72] 0 0 0 0 0 0 0 0 0 0 ...
As we include separate fixed effects for age groups and districts (rather than their interactions), we construct the corresponding formulas manually for the control list.
Next, we fit a model with a power-adjusted contact matrix. Since the power parameter is estimated via profile likelihood, model fitting takes a bit longer.
### fit power-adjusted contact matrix via profile likelihoodptm <-proc.time()noroEEcpowerFit <-fitC(noroEEcFit, Cgrouped, normalize =TRUE, truncate =TRUE)
initial value 14420.425960
iter 2 value 14418.693846
iter 3 value 14414.397722
iter 4 value 14408.695992
iter 5 value 14405.588920
iter 6 value 14405.438849
iter 7 value 14405.326531
iter 8 value 14405.312660
iter 8 value 14405.312513
iter 8 value 14405.312513
final value 14405.312513
converged
Comparing their AIC values shows that parameterising the contact matrix improves the model fit.
6.3 Using RTMB
In RTMB we first use the same flat matrix strucutre and make use of the same indicator matrices, as well as seasonal terms matrix.
First, we create a data list with the flat norovirus counts, add all indicator matrices and the contact matrix.
noro_all <-observed(noroBEall)## Make data suitable for objective function n_time <-nrow(noro_all)n_unit <-ncol(noro_all)data <-list(Y = noro_all, # noro case matrix n_time = n_time,n_unit = n_unit, # all district x groupsdistricts = DISTRICTS,n_districts = NDISTRICTS,groups = GROUPS,n_groups = NGROUPS,W_n =neighbourhood(noroBEall)[1:NDISTRICTS, 1:NDISTRICTS], # neighborhood matrixCgrouped = Cgrouped, # non normalized unexpanded contact matrix EV =eigen(Cgrouped /rowSums(Cgrouped), symmetric =FALSE), MMR = MMR, # model matrix with district indicatorsMMG= MMG, # model matrix with group indicatorsMMgS = MMgS # model matrix with age specific season)
We build a helper function for the mean, where we use the Kronecker product to get a big weight matrix (N * M x N * M).
Then we build the endemic- and neighbourhood component matrices out of the parameters by looping through either districts or groups, to add the corresponding parameters (also age specific season) to the matrix of the global intercept of each component (therefore we add a 0 to the district and group parameter vectors at the beginning).
meanEEc <-function(parms, data) { n_unit <- data$n_unit n_time <- data$n_time groups <- data$groups districts <- data$districts## transform parameters neweights.d <-exp(parms$neweights.logd) power <-exp(parms$logpower)## power-law including own region with unit weight W <- (data$W_n +1)^(-neweights.d)## age group contact matrix EV <- data$EV C <- EV$vectors %*%diag(EV$values^power) %*%solve(EV$vectors)## Combine them WC <-kronecker(C, W)## normalize rows WCnorm <- WC/rowSums(WC)REPORT(WCnorm)## create full parameter vectors ## 0 added global intercept without district- or group-specific effect end.group <-group_par(parms, "end", "", which = groups[-1], threeparts =FALSE) end.group_full <-c(0, end.group) end.district <-group_par(parms, "end", "", which = districts[-1], threeparts =FALSE) end.district_full <-c(0, end.district) ne.group <-group_par(parms, "ne", "", which = groups[-1], threeparts =FALSE) ne.group_full <-c(0, ne.group) ne.district <-group_par(parms, "ne", "", which = districts[-1], threeparts =FALSE) ne.district_full <-c(0, ne.district)## ENDEMIC## sesonality sinS <-group_par(parms, "end", "sin", which = groups) cosS <-group_par(parms, "end", "cos", which = groups) end.mat <-matrix(parms$end.1, n_time, n_unit)for (g inseq_along(groups)) { end.mat <- end.mat + end.group_full[[g]] * MMG[[g]] }for (d inseq_along(districts)) { end.mat <- end.mat + end.district_full[[d]] * MMR[[d]] }for (g inseq_along(groups)) { end.mat <- end.mat + sinS[[g]] * MMgS[[2*g-1]] end.mat <- end.mat + cosS[[g]] * MMgS[[2*g]] } end.exppred <-exp(end.mat) end <- end.exppred[-1, ]## EPIDEMIC ne.mat <-matrix(parms$ne.1, n_time, n_unit)for (g inseq_along(groups)) { ne.mat <- ne.mat + ne.group_full[[g]] * MMG[[g]] }for (d inseq_along(districts)) { ne.mat <- ne.mat + ne.district_full[[d]] * MMR[[d]] }## data Y_lag <- data$Y[1:(n_time -1), ] ne.exppred <-exp(ne.mat) ne <- ne.exppred[-1, ] * (Y_lag %*% WCnorm)return(list(mu = end + ne,end = end,epi = ne,epi.ne = ne,end.exppred = end.exppred,ne.exppred = ne.exppred,neW = WCnorm))}objectiveEEcpower <-function(parms, data) { Y <-OBS(data$Y) fit <-meanEEc(parms, data)## negative log likelihood nll <--sum(dnbinom2(Y[-1, ], mu = fit$mu, size =exp(parms$log_overdispersion), log =TRUE)) nll}
Since we want to compare the parameter estimates we name the parameters similar to the hhh4 model and order them accordingly.
## create group specific seasonality parameters sinS <-rep(0, data$n_groups)names(sinS) <-paste0("end.sin.", GROUPS)cosS <-rep(0, data$n_groups)names(cosS) <-paste0("end.cos.", GROUPS)## order them as in hhh4contacts S <-c(sinS, cosS)[order(c(seq_along(sinS),seq_along(cosS)))]## create group and district specific (named) fixed effectsne.group <-rep(0, data$n_groups -1)names(ne.group) <-paste0("ne.", GROUPS[-1])ne.district =rep(0, data$n_districts -1)names(ne.district) <-paste0("ne.", DISTRICTS[-1])end.group <-rep(0.1, data$n_groups -1)names(end.group) <-paste0("end.", GROUPS[-1])end.district =rep(0, data$n_districts -1)names(end.district) <-paste0("end.", DISTRICTS[-1])## Parameter objects are gathered in a list that also serves as initial guess when fitting the model:parms <-c(list(ne.1 =0), # intercept phi ne.group, ne.district,list(end.1 =0), # intercept nu end.group, end.district, S,list(neweights.logd =log(2), # distance decay parameterlog_overdispersion =0, # overdispersionlogpower =0# power contact matrix ))
We fit the model as usual and calculate uncertainties.
## initate runtimeptm <-proc.time()## The objective function is processed by RTMB using the callobjEEcpower <-MakeADFun(cmb(objectiveEEcpower, data), parms)# Fitting the model## runtime without processing objectiveptm_wo <-proc.time()## Optimize the model using nlminboptEEcpower <-nlminb(objEEcpower$par, objEEcpower$fn, objEEcpower$gr)optEEcpower$objective## Uncertainties are now calculated using:sdrEEcpower <-sdreport(objEEcpower)runtimeEEcpower <-proc.time() - ptmruntime_woEEcpower <-proc.time() - ptm_wosdrEEcpower
6.4 Comparison Table
Parameter estimates and runtime for <code>hhh4</code> and <code>RTMB</code> implementations with standard errors (SE) in brackets and differences between both implementations.
hhh4
RTMB
diff.
Parameters
ne.1
-0.641909877 (0.141895268)
-0.642009479 (0.173119403)
9.96e-05 (-0.0312)
ne.05-14
-1.696712794 (0.27865492)
-1.696564373 (0.314039136)
-0.000148 (-0.0354)
ne.15-24
-1.437598286 (0.274175641)
-1.437430317 (0.278818989)
-0.000168 (-0.00464)
ne.25-44
-2.051832422 (0.274475407)
-2.051515251 (0.293541579)
-0.000317 (-0.0191)
ne.45-64
-0.759469413 (0.128301562)
-0.759300885 (0.16039412)
-0.000169 (-0.0321)
ne.65+
0.424232248 (0.107243456)
0.424323871 (0.107370148)
-9.16e-05 (-0.000127)
ne.frkr
-1.148536351 (0.24596229)
-1.148434916 (0.245905568)
-0.000101 (5.67e-05)
ne.lich
0.214200753 (0.159126561)
0.214218508 (0.161468963)
-1.78e-05 (-0.00234)
ne.mahe
0.108530217 (0.189265185)
0.108384927 (0.192469371)
0.000145 (-0.0032)
ne.mitt
0.01382475 (0.148209959)
0.013787247 (0.149489667)
3.75e-05 (-0.00128)
ne.neuk
0.373250527 (0.139316411)
0.373201762 (0.141019048)
4.88e-05 (-0.0017)
ne.pank
0.71509765 (0.135106967)
0.71506207 (0.138147689)
3.56e-05 (-0.00304)
ne.rein
0.670839816 (0.130456726)
0.670801334 (0.132777451)
3.85e-05 (-0.00232)
ne.span
0.207200784 (0.14372353)
0.207111067 (0.144611796)
8.97e-05 (-0.000888)
ne.zehl
0.938878342 (0.128337523)
0.938817403 (0.130784139)
6.09e-05 (-0.00245)
ne.scho
0.433492551 (0.131180989)
0.433502916 (0.132882988)
-1.04e-05 (-0.0017)
ne.trko
0.746046993 (0.132009801)
0.74606413 (0.134342501)
-1.71e-05 (-0.00233)
end.1
-0.885741801 (0.11271468)
-0.885713676 (0.115787848)
-2.81e-05 (-0.00307)
end.05-14
-1.722361281 (0.17137752)
-1.722388366 (0.171617357)
2.71e-05 (-0.00024)
end.15-24
-1.149537138 (0.14970184)
-1.149597615 (0.151328713)
6.05e-05 (-0.00163)
end.25-44
0.210326615 (0.11731632)
0.210229101 (0.117667976)
9.75e-05 (-0.000352)
end.45-64
-0.113851153 (0.128140949)
-0.113915339 (0.128126898)
6.42e-05 (1.41e-05)
end.65+
-0.130972249 (0.182180873)
-0.13105175 (0.185457493)
7.95e-05 (-0.00328)
end.frkr
-0.073203685 (0.115676445)
-0.073172195 (0.117151009)
-3.15e-05 (-0.00147)
end.lich
-0.042038833 (0.136533017)
-0.042025109 (0.138028064)
-1.37e-05 (-0.0015)
end.mahe
-0.022641648 (0.14025298)
-0.022514572 (0.142999517)
-0.000127 (-0.00275)
end.mitt
-0.030395757 (0.121063314)
-0.030339901 (0.12195768)
-5.59e-05 (-0.000894)
end.neuk
-0.330771933 (0.135842862)
-0.330697636 (0.136137268)
-7.43e-05 (-0.000294)
end.pank
0.372286245 (0.122408053)
0.372351157 (0.123492019)
-6.49e-05 (-0.00108)
end.rein
-0.906120335 (0.205133998)
-0.906041995 (0.205763093)
-7.83e-05 (-0.000629)
end.span
-0.44074959 (0.139420987)
-0.440674628 (0.139425862)
-7.5e-05 (-4.87e-06)
end.zehl
-0.033357158 (0.139956804)
-0.033222968 (0.140116816)
-0.000134 (-0.00016)
end.scho
0.07972706 (0.119481083)
0.079719347 (0.119846356)
7.71e-06 (-0.000365)
end.trko
-0.939677801 (0.216505885)
-0.939711985 (0.216870748)
3.42e-05 (-0.000365)
end.sin(2 * pi * t/52).00-04
0.719887911 (0.075644204)
0.719835427 (0.0827625)
5.25e-05 (-0.00712)
end.cos(2 * pi * t/52).00-04
-0.875776299 (0.069228606)
-0.875790888 (0.070424068)
1.46e-05 (-0.0012)
end.sin(2 * pi * t/52).05-14
0.456080055 (0.140552978)
0.456021921 (0.141845219)
5.81e-05 (-0.00129)
end.cos(2 * pi * t/52).05-14
-0.073455731 (0.145525225)
-0.073518521 (0.14593679)
6.28e-05 (-0.000412)
end.sin(2 * pi * t/52).15-24
0.089673925 (0.099669375)
0.089654211 (0.100140342)
1.97e-05 (-0.000471)
end.cos(2 * pi * t/52).15-24
-0.612169348 (0.105351559)
-0.612164512 (0.106259368)
-4.84e-06 (-0.000908)
end.sin(2 * pi * t/52).25-44
0.066929973 (0.05401224)
0.066924674 (0.054299443)
5.3e-06 (-0.000287)
end.cos(2 * pi * t/52).25-44
-0.518159038 (0.054961138)
-0.518148667 (0.055168543)
-1.04e-05 (-0.000207)
end.sin(2 * pi * t/52).45-64
0.049550469 (0.083096145)
0.049570528 (0.083589977)
-2.01e-05 (-0.000494)
end.cos(2 * pi * t/52).45-64
-0.765506844 (0.077147929)
-0.765497175 (0.077333814)
-9.67e-06 (-0.000186)
end.sin(2 * pi * t/52).65+
-0.141382806 (0.129709515)
-0.141342433 (0.130560298)
-4.04e-05 (-0.000851)
end.cos(2 * pi * t/52).65+
-1.350585868 (0.126488017)
-1.350547056 (0.128515774)
-3.88e-05 (-0.00203)
neweights.logd
0.752830336 (0.067450869)
0.752835373 (0.067613271)
-5.04e-06 (-0.000162)
-log(overdisp)
1.127533356 (0.053653197)
1.127554936 (0.053654161)
-2.16e-05 (-9.63e-07)
logpower
-0.703612013 (0.165639272)
-0.703801931 (0.165712013)
0.00019 (-7.27e-05)
Metrics
Log-Likelihood
-14405.312513112
-14405.31251669
3.578e-06
Runtime (seconds)
149.693
2.408
147.285
The estimates are again quite similar. The use of profile likelihood to estimate the power parameter for the contact matrix, slows the fitting with hhh4contacts down. It takes 62 times longer. This makes quite a difference, especially when the dimensions are larger, as this can create a bottleneck in rolling forecasts later on.
6.5 Test Array Structure Togeteher with Design Matrix Structure
Since building the indicator matrices is complex and not necessary for RTMB models, we make use of a design matrix approach using the meta data, fromulas and model.matrix().
Therefore, we first need to build a data frame of the meta data (dimensions and covariates). We use the district and age group vectors from the beginning.
dimens <-list(time =1:n_time, districts = DISTRICTS, groups = GROUPS)meta_data <-expand.grid(dimens)head(meta_data)
Now, using formulas, we create a model.matrix for each component. The seasonal terms are given in the formula directly.
Xend <-model.matrix(~ districts + groups +sin(2* pi * (time -1) /52):groups +cos(2* pi * (time -1) /52):groups, data = meta_data)Xne <-model.matrix(~ districts + groups, data = meta_data)colnames(Xend) <-paste0("end.", colnames(Xend))colnames(Xne) <-paste0("ne.", colnames(Xne))
We tweak the names of the columns (parameters) to be the same as in the first RTMB model to make them comparable.
## remove covariable group name from namessub_co_group <-function(MM) {colnames(MM) <-gsub("districts", "", colnames(MM))colnames(MM) <-gsub("groups", "", colnames(MM))colnames(MM) <-gsub("\\(Intercept\\)", "1", colnames(MM))return(MM)}Xend <-sub_co_group(Xend)Xne <-sub_co_group(Xne)## remove sinus and cosinus formulationscolnames(Xend) <-gsub("^end\\.([^:]+):(sin|cos)\\(.*\\)$", "end.\\2.\\1", colnames(Xend))## order the columns as in first rtmb modelXend <- Xend[, grep("end*", names(parms), value =TRUE)]Xne <- Xne[, grep("^ne\\.", names(parms), value =TRUE)]
For the observations \(Y\), we use a 3-dimensional array structure instead of the flattened matrix from before.
# List of sts objects per age groupnoroBElist <-noroBE(by ="all", flatten =FALSE, agegroups = GROUPING)# Get observed matricesobs_matrices <-lapply(noroBElist, observed)# Combine matrices into a 3-dimensional array [time, district, agegroup]Y <-simplify2array(obs_matrices)
We can see directly that the function to calculate \(\mu_{git}\) is less complex and shorter. The model matrices (\(X_{end}, X_{ne}\)) are multiplied with the corresponding parameter vectors and reshaped into the multidemnsional array structure of the obeservations \(Y\).
The weight matrices \(w_{ij}\) and \(c_{hg}\) are normalized sperarately. To get the weighted lagged \(Y_{g,i,t-1}\), we first multiply \(w_{ij}\) to the lagged \(Y_{i,t-1}\) while the third dimension (groups) is fixed and second multiply the resulting spatially weighted \(Y_{g,t-1}\), while the second dimension (districts) is fixed, with \(c_{hg}\). This is the same as taking the Kronecker product of both weight matrices, normalize the resulting ((gi) x (gi)) matrix and multiplying it with the flattened \(Y_{(g*i),t-1}\).
This way we could also scale up the dimensions as needed.
meanEEcMM <-function(parms,data) { Y <-OBS(data$Y)## EPIDEMIC## coefficients end.par <-do.call(c, parms[colnames(data$Xend)]) end.pre <- data$Xend %*% end.par end.exppred <-exp(end.pre)dim(end.exppred) <-dim(Y) end <- end.exppred[-1, , ]## ENDEMIC## power-law including own region with unit weight W <- (data$W_n +1)^(-exp(parms$neweights.logd)) Wnorm <- W /rowSums(W)## age group contact matrix EV <- data$EV C <- EV$vectors %*%diag(EV$values^exp(parms$logpower)) %*%solve(EV$vectors) Cnorm <- C /rowSums(C)## coefficients ne.par <-do.call(c, parms[colnames(data$Xne)]) ne.pre <- data$Xne %*% ne.par ne.exppred <-exp(ne.pre)dim(ne.exppred) <-dim(Y)## lagged Y Y_lag <- Y[1:(data$n_time -1), , ]## empty Y weighted Y_weigh <- Y_lag## sequential dimension wise weighting## spatialfor (g inseq_len(data$n_groups)) { Y_weigh[, , g] <- Y_lag[, , g] %*% Wnorm # another helper (scale up to n dimensions) }## groupsfor (d inseq_len(data$n_districts)) { Y_weigh[, d, ] <- Y_weigh[, d, ] %*% Cnorm } ne <- ne.exppred[-1, , ] * Y_weighreturn(list(mu = end + ne,end = end,epi = ne,epi.ne = ne,end.exppred = end.exppred,ne.exppred = ne.exppred,neW = Wnorm,neC = Cnorm))}objectiveEEcpowerMM <-function(parms, data) { Y <-OBS(data$Y)## get the fit fit <-meanEEcMM(parms, data)## negative log likelihood nll <--sum(dnbinom2(Y[-1, , ], mu = fit$mu, size =exp(parms$log_overdispersion), log =TRUE)) nll}
## initate runtimeptm <-proc.time()## The objective function is processed by RTMB using the callobjEEcpowerMM <-MakeADFun(cmb(objectiveEEcpowerMM, dataS), parms)# Fitting the model## runtime without processing objectiveptm_wo <-proc.time()## Optimize the model using nlminboptEEcpowerMM <-nlminb(objEEcpowerMM$par, objEEcpowerMM$fn, objEEcpowerMM$gr)## Uncertainties are now calculated using:sdrEEcpowerMM <-sdreport(objEEcpowerMM)runtimeEEcpowerMM <-proc.time() - ptmruntime_woEEcpowerMM <-proc.time() - ptm_wo
The difference between the parameter estimates of both RTMB models and their negative log-likelihood is quite small.
Parameter estimates and runtime for <code>hhh4</code> and <code>RTMB</code> implementations with standard errors (SE) in brackets and differences between both implementations.
hhh4
RTMB
diff.
Parameters
ne.1
-0.641909877 (0.141895268)
-0.641957086 (0.173110111)
4.72e-05 (-0.0312)
ne.05-14
-1.696712794 (0.27865492)
-1.696187836 (0.313975277)
-0.000525 (-0.0353)
ne.15-24
-1.437598286 (0.274175641)
-1.437211943 (0.278828189)
-0.000386 (-0.00465)
ne.25-44
-2.051832422 (0.274475407)
-2.051142576 (0.293508605)
-0.00069 (-0.019)
ne.45-64
-0.759469413 (0.128301562)
-0.759122187 (0.16038972)
-0.000347 (-0.0321)
ne.65+
0.424232248 (0.107243456)
0.42436964 (0.107382446)
-0.000137 (-0.000139)
ne.frkr
-1.148536351 (0.24596229)
-1.148553619 (0.245910904)
1.73e-05 (5.14e-05)
ne.lich
0.214200753 (0.159126561)
0.213885572 (0.161478491)
0.000315 (-0.00235)
ne.mahe
0.108530217 (0.189265185)
0.108193734 (0.192479615)
0.000336 (-0.00321)
ne.mitt
0.01382475 (0.148209959)
0.013643994 (0.149485065)
0.000181 (-0.00128)
ne.neuk
0.373250527 (0.139316411)
0.373127037 (0.141010777)
0.000123 (-0.00169)
ne.pank
0.71509765 (0.135106967)
0.714954594 (0.138140153)
0.000143 (-0.00303)
ne.rein
0.670839816 (0.130456726)
0.670684959 (0.132768099)
0.000155 (-0.00231)
ne.span
0.207200784 (0.14372353)
0.206961656 (0.144604803)
0.000239 (-0.000881)
ne.zehl
0.938878342 (0.128337523)
0.938686378 (0.130775331)
0.000192 (-0.00244)
ne.scho
0.433492551 (0.131180989)
0.433337127 (0.132874319)
0.000155 (-0.00169)
ne.trko
0.746046993 (0.132009801)
0.745847151 (0.134333554)
2e-04 (-0.00232)
end.1
-0.885741801 (0.11271468)
-0.885715194 (0.115788987)
-2.66e-05 (-0.00307)
end.05-14
-1.722361281 (0.17137752)
-1.722560848 (0.171622501)
2e-04 (-0.000245)
end.15-24
-1.149537138 (0.14970184)
-1.149739576 (0.15134887)
0.000202 (-0.00165)
end.25-44
0.210326615 (0.11731632)
0.210087857 (0.117674962)
0.000239 (-0.000359)
end.45-64
-0.113851153 (0.128140949)
-0.114067211 (0.128133939)
0.000216 (7.01e-06)
end.65+
-0.130972249 (0.182180873)
-0.131192302 (0.185490919)
0.00022 (-0.00331)
end.frkr
-0.073203685 (0.115676445)
-0.073089473 (0.117164223)
-0.000114 (-0.00149)
end.lich
-0.042038833 (0.136533017)
-0.041813842 (0.138040315)
-0.000225 (-0.00151)
end.mahe
-0.022641648 (0.14025298)
-0.022406018 (0.14301862)
-0.000236 (-0.00277)
end.mitt
-0.030395757 (0.121063314)
-0.030269544 (0.121970964)
-0.000126 (-0.000908)
end.neuk
-0.330771933 (0.135842862)
-0.330749236 (0.136160208)
-2.27e-05 (-0.000317)
end.pank
0.372286245 (0.122408053)
0.372404999 (0.123506033)
-0.000119 (-0.0011)
end.rein
-0.906120335 (0.205133998)
-0.906181882 (0.205832337)
6.15e-05 (-0.000698)
end.span
-0.44074959 (0.139420987)
-0.440562674 (0.139434964)
-0.000187 (-1.4e-05)
end.zehl
-0.033357158 (0.139956804)
-0.033219415 (0.140137157)
-0.000138 (-0.00018)
end.scho
0.07972706 (0.119481083)
0.079785017 (0.119857785)
-5.8e-05 (-0.000377)
end.trko
-0.939677801 (0.216505885)
-0.939512936 (0.216857465)
-0.000165 (-0.000352)
end.sin(2 * pi * t/52).00-04
0.719887911 (0.075644204)
0.719776096 (0.082752194)
0.000112 (-0.00711)
end.cos(2 * pi * t/52).00-04
-0.875776299 (0.069228606)
-0.875711894 (0.070421811)
-6.44e-05 (-0.00119)
end.sin(2 * pi * t/52).05-14
0.456080055 (0.140552978)
0.455938352 (0.14185353)
0.000142 (-0.0013)
end.cos(2 * pi * t/52).05-14
-0.073455731 (0.145525225)
-0.073459259 (0.145953006)
3.53e-06 (-0.000428)
end.sin(2 * pi * t/52).15-24
0.089673925 (0.099669375)
0.089679413 (0.100147108)
-5.49e-06 (-0.000478)
end.cos(2 * pi * t/52).15-24
-0.612169348 (0.105351559)
-0.612185529 (0.106266327)
1.62e-05 (-0.000915)
end.sin(2 * pi * t/52).25-44
0.066929973 (0.05401224)
0.066954723 (0.054304638)
-2.48e-05 (-0.000292)
end.cos(2 * pi * t/52).25-44
-0.518159038 (0.054961138)
-0.518151337 (0.055172264)
-7.7e-06 (-0.000211)
end.sin(2 * pi * t/52).45-64
0.049550469 (0.083096145)
0.04958936 (0.08359812)
-3.89e-05 (-0.000502)
end.cos(2 * pi * t/52).45-64
-0.765506844 (0.077147929)
-0.765492694 (0.07733986)
-1.42e-05 (-0.000192)
end.sin(2 * pi * t/52).65+
-0.141382806 (0.129709515)
-0.141287087 (0.130581347)
-9.57e-05 (-0.000872)
end.cos(2 * pi * t/52).65+
-1.350585868 (0.126488017)
-1.350624782 (0.128530254)
3.89e-05 (-0.00204)
neweights.logd
0.752830336 (0.067450869)
0.752807519 (0.067618895)
2.28e-05 (-0.000168)
-log(overdisp)
1.127533356 (0.053653197)
1.127537861 (0.053653661)
-4.5e-06 (-4.63e-07)
logpower
-0.703612013 (0.165639272)
-0.7038947 (0.165697968)
0.000283 (-5.87e-05)
Metrics
Log-Likelihood
-14405.312513112
-14405.312520802
7.69e-06
Runtime (seconds)
149.693
2.408
147.285
Although the comparison of the model matrix RTMB model with the hhh4 one looks similar to our first comparison table, using the model.matrix approach makes the code much easier to follow as it is more intuitive. Formulating the model is also easier and less error-prone using this approach, since one does not have to construct the flattened indicator matrices as before, and it is known by many R users as other modelling packages such as lme4 use model matrices too. Moreover, we don’t have to construct the parameter names.
Meyer, Sebastian, and Leonhard Held. 2017. “Incorporating Social Contact Data in Spatio-Temporal Models for Infectious Disease Spread.”Biostatistics 18 (2): 338–51. https://doi.org/10.1093/biostatistics/kxw051.
Mossong, Joël, Niel Hens, Mark Jit, Philippe Beutels, Kari Auranen, Rafael Mikolajczyk, Marco Massari, et al. 2008. “Social Contacts and Mixing Patterns Relevant to the Spread of Infectious Diseases.”PLOS Medicine 5 (3): e74. https://doi.org/10.1371/journal.pmed.0050074.