R Under development (unstable) (2023-10-22 r85388 ucrt) -- "Unsuffered Consequences"
Copyright (C) 2023 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library("hhh4contacts")
Loading required package: surveillance
Loading required package: sp
Loading required package: xtable
This is surveillance 1.21.1; see 'package?surveillance' or
https://surveillance.R-Forge.R-project.org/ for an overview.
> library("surveillance")
> data("measlesWeserEms")
> 
> ## basic power-law model
> WPL <- W_powerlaw(maxlag = 5, normalize = TRUE)
> measlesModel <- list(
+     end = list(f = addSeason2formula(~1), offset = population(measlesWeserEms)),
+     ar = list(f = ~1),
+     ne = list(f = ~1 + log(pop), weights = WPL),
+     family = "NegBin1", data = list(pop = population(measlesWeserEms))
+ )
> measlesFit <- hhh4(stsObj = measlesWeserEms, control = measlesModel)
> 
> ## matrix of adjacency orders
> nbmat <- neighbourhood(measlesWeserEms)
> 
> ## fake group-specific power law (single group)
> WPLgfake <- addGroups2WFUN(WPL, groups = factor(rep("d", ncol(nbmat))))
NOTE: 'groups' has only one level ...
> stopifnot(identical(WPL$w(0.5, nbmat), WPLgfake$w(0.5, nbmat)),
+           identical(WPL$dw(0.5, nbmat), WPLgfake$dw(0.5, nbmat)[[1]]),
+           identical(WPL$d2w(0.5, nbmat), WPLgfake$d2w(0.5, nbmat)[[1]]))
> stopifnot(all.equal(
+     measlesFit,
+     update(measlesFit, ne = list(weights = WPLgfake), use.estimates = FALSE),
+     ignore = "control"))
> 
> ### uncomment below to check derivatives with multiple groups
> ### (time-consuming and analytical derivatives already verified)
> ## WPLgroups <- addGroups2WFUN(WPL, factor(substr(colnames(nbmat), 1, 4) == "0345"))
> ## pars_groups <- c(0.5, 2)
> ## dwnum <- sapply(seq_along(nbmat), function (i)
> ##    numDeriv::grad(function (wpar) WPLgroups$w(wpar, nbmat)[i], x = pars_groups))
> ## stopifnot(all.equal(dwnum[1,], c(WPLgroups$dw(pars_groups, nbmat)[[1]])),
> ##          all.equal(dwnum[2,], c(WPLgroups$dw(pars_groups, nbmat)[[2]])))
> ## d2wnum <- sapply(seq_along(nbmat), function (i)
> ##    numDeriv::hessian(function (wpar) WPLgroups$w(wpar, nbmat)[i], x = pars_groups))
> ## stopifnot(all.equal(d2wnum[1,], c(WPLgroups$d2w(pars_groups, nbmat)[[1]])),
> ##          all.equal(d2wnum[4,], c(WPLgroups$d2w(pars_groups, nbmat)[[3]])),
> ##          abs(c(d2wnum[c(2,3),]) - 0) < .Machine$double.eps)
> 
> ### check score vector and Fisher info of all model parameters
> ## measlesModelGrouped <- modifyList(measlesModel, list(ne=list(weights=WPLgroups)))
> ## hhh4(measlesWeserEms, measlesModelGrouped, check.analyticals = TRUE)
> 
> proc.time()
   user  system elapsed 
   2.42    0.43    2.79