library(serocalculator)
library(Hmisc)
#>
#> Attaching package: 'Hmisc'
#> The following objects are masked from 'package:base':
#>
#> format.pval, units
library(dplyr)
#>
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:Hmisc':
#>
#> src, summarize
#> The following objects are masked from 'package:stats':
#>
#> filter, lag
#> The following objects are masked from 'package:base':
#>
#> intersect, setdiff, setequal, union
# load longitudinal parameter MC sample
fs::path_package(
"extdata",
"pertussis-iggpt2.rda",
package = "serocalculator") |>
load() # output: longpar.mc
lp = longpar.mc$IgG
day2yr <- 365.25;
ln.pars <- data.frame(
y1 = lp$y1, # peak levels
alpha = day2yr * lp$alpha, # decay rate (per year)
d = lp$r -1); # shape (r-1 offset from exp)
library(readr)
#read simulated data
raw.data = fs::path_package(
"extdata",
"SimulDat.txt",
package = "serocalculator") |>
readr::read_table() |>
rename(AGE = '"AGE"', IgG = '"IgG"')
#>
#> ── Column specification ────────────────────────────────────────────────────────
#> cols(
#> `"AGE"` = col_double(),
#> `"IgG"` = col_double()
#> )
cat.bnd = seq(0,80,by=5);
n.cats = length(cat.bnd)-1; # define age categories
age.cat = cut(raw.data$AGE, breaks = cat.bnd, right = FALSE) |> as.numeric()
clean.data = data.frame(
age = floor(raw.data$AGE) + 0.5,
IgG = raw.data$IgG,
cat = age.cat)
# define global parameters
cond <- data.frame(
nu=2.5, # B noise
eps=0.3, # M noise
y.low=2.5, # low cutoff
y.high=5e4); # high cutoff
lambda = 0.015 # initial estimate: starting value
log.lambda = log(lambda)
log.lmin = log(lambda/10);
log.lmax=log(10*lambda) # seroincidence rate interval
log.lambda.est = matrix(0, nrow=n.cats, ncol=6,
dimnames = list(
rep("",n.cats),
c("","ml.est", "lwr0.05", "upr0.95",
"LLF","N iter.")))
# loop over age categories
system.time(
for (k.cat in 1:n.cats) # per age category
{
cat(k.cat,"\n")
y.cat = as.vector(subset(clean.data,subset=clean.data$cat==as.character(k.cat),
select="IgG"))$IgG;
a.cat = as.vector(subset(clean.data,subset=clean.data$cat==as.character(k.cat),
select="age"))$age;
cs.data <- data.frame(y=y.cat, a=a.cat);
objfunc <- function(llam){
# add terms, e.g. for other antibodies
# res < fdev(llam,cs1,ln1,cond1) + fdev(llam,cs2,ln2,cond2) + ...
# return(res)
return(fdev(llam,cs.data,ln.pars,cond));
}
# seroincidence estimation
fit = nlm(objfunc,log.lambda,
hessian=TRUE,print.level=0,stepmax=(log.lmax-log.lmin)/4);
log.lambda.est[k.cat,] = c(k.cat,
fit$estimate,
fit$estimate + qnorm(c(0.05))*sqrt(1/fit$hessian),
fit$estimate + qnorm(c(0.95))*sqrt(1/fit$hessian),
fit$minimum,
fit$iterations)
log.lambda=log.lambda.est[k.cat,2]
})
#> 1
#> 2
#> 3
#> 4
#> 5
#> 6
#> 7
#> 8
#> 9
#> 10
#> 11
#> 12
#> 13
#> 14
#> 15
#> 16
#> user system elapsed
#> 33.838 0.011 33.850
# check results in log.lambda.est:
pander::pander(log.lambda.est)
| 1 |
-4.111 |
-4.782 |
-3.439 |
89.69 |
3 |
| 2 |
-4.495 |
-4.97 |
-4.02 |
123.6 |
5 |
| 3 |
-4.358 |
-4.71 |
-4.007 |
202.7 |
4 |
| 4 |
-4.116 |
-4.386 |
-3.845 |
262.9 |
4 |
| 5 |
-4.194 |
-4.451 |
-3.937 |
332.3 |
4 |
| 6 |
-3.983 |
-4.205 |
-3.761 |
395.4 |
4 |
| 7 |
-4.082 |
-4.3 |
-3.865 |
413.7 |
4 |
| 8 |
-4.449 |
-4.683 |
-4.216 |
330.1 |
5 |
| 9 |
-4.347 |
-4.566 |
-4.129 |
379.4 |
4 |
| 10 |
-4.376 |
-4.592 |
-4.161 |
406.5 |
3 |
| 11 |
-4.148 |
-4.34 |
-3.956 |
462.8 |
4 |
| 12 |
-4.273 |
-4.472 |
-4.074 |
455.5 |
4 |
| 13 |
-4.092 |
-4.278 |
-3.906 |
508.8 |
4 |
| 14 |
-4.128 |
-4.314 |
-3.942 |
508.6 |
3 |
| 15 |
-4.342 |
-4.531 |
-4.152 |
440.1 |
4 |
| 16 |
-4.27 |
-4.458 |
-4.082 |
518.1 |
4 |